nc_beta.qbk 7.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205
  1. [section:nc_beta_dist Noncentral Beta Distribution]
  2. ``#include <boost/math/distributions/non_central_beta.hpp>``
  3. namespace boost{ namespace math{
  4. template <class RealType = double,
  5. class ``__Policy`` = ``__policy_class`` >
  6. class non_central_beta_distribution;
  7. typedef non_central_beta_distribution<> non_central_beta;
  8. template <class RealType, class ``__Policy``>
  9. class non_central_beta_distribution
  10. {
  11. public:
  12. typedef RealType value_type;
  13. typedef Policy policy_type;
  14. // Constructor:
  15. non_central_beta_distribution(RealType alpha, RealType beta, RealType lambda);
  16. // Accessor to shape parameters:
  17. RealType alpha()const;
  18. RealType beta()const;
  19. // Accessor to non-centrality parameter lambda:
  20. RealType non_centrality()const;
  21. };
  22. }} // namespaces
  23. The noncentral beta distribution is a generalization of the __beta_distrib.
  24. It is defined as the ratio
  25. [expression X = [chi][sub m][super 2]([lambda]) \/ ([chi][sub m][super 2]([lambda])
  26. + [chi][sub n][super 2])]
  27. where [role serif_italic [chi][sub m][super 2]([lambda])]
  28. is a noncentral [role serif_italic [chi][super 2]]
  29. random variable with /m/ degrees of freedom, and [chi][sub n][super 2]
  30. is a central [role serif_italic [chi][super 2] ] random variable with /n/ degrees of freedom.
  31. This gives a PDF that can be expressed as a Poisson mixture
  32. of beta distribution PDFs:
  33. [equation nc_beta_ref1]
  34. where P(i;[lambda]\/2) is the discrete Poisson probablity at /i/, with mean
  35. [lambda]\/2, and I[sub x][super ']([alpha], [beta]) is the derivative of
  36. the incomplete beta function. This leads to the usual form of the CDF
  37. as:
  38. [equation nc_beta_ref2]
  39. The following graph illustrates how the distribution changes
  40. for different values of [lambda]:
  41. [graph nc_beta_pdf]
  42. [h4 Member Functions]
  43. non_central_beta_distribution(RealType a, RealType b, RealType lambda);
  44. Constructs a noncentral beta distribution with shape parameters /a/ and /b/
  45. and non-centrality parameter /lambda/.
  46. Requires a > 0, b > 0 and lambda >= 0, otherwise calls __domain_error.
  47. RealType alpha()const;
  48. Returns the parameter /a/ from which this object was constructed.
  49. RealType beta()const;
  50. Returns the parameter /b/ from which this object was constructed.
  51. RealType non_centrality()const;
  52. Returns the parameter /lambda/ from which this object was constructed.
  53. [h4 Non-member Accessors]
  54. Most of the [link math_toolkit.dist_ref.nmp usual non-member accessor functions]
  55. are supported: __cdf, __pdf, __quantile, __mean, __variance, __sd,
  56. __median, __mode, __hazard, __chf, __range and __support.
  57. Mean and variance are implemented using hypergeometric pfq functions and relations given in
  58. [@http://reference.wolfram.com/mathematica/ref/NoncentralBetaDistribution.html Wolfram Noncentral Beta Distribution].
  59. However, the following are not currently implemented:
  60. __skewness, __kurtosis and __kurtosis_excess.
  61. The domain of the random variable is \[0, 1\].
  62. [h4 Accuracy]
  63. The following table shows the peak errors
  64. (in units of [@http://en.wikipedia.org/wiki/Machine_epsilon epsilon])
  65. found on various platforms with various floating point types.
  66. The failures in the comparison to the [@http://www.r-project.org/ R Math library],
  67. seem to be mostly in the corner cases when the probablity would be very small.
  68. Unless otherwise specified any floating-point type that is narrower
  69. than the one shown will have __zero_error.
  70. [table_non_central_beta_CDF]
  71. [table_non_central_beta_CDF_complement]
  72. Error rates for the PDF, the complement of the CDF and for the quantile
  73. functions are broadly similar.
  74. [h4 Tests]
  75. There are two sets of test data used to verify this implementation:
  76. firstly we can compare with a few sample values generated by the
  77. [@http://www.r-project.org/ R library].
  78. Secondly, we have tables of test data, computed with this
  79. implementation and using interval arithmetic - this data should
  80. be accurate to at least 50 decimal digits - and is the used for
  81. our accuracy tests.
  82. [h4 Implementation]
  83. The CDF and its complement are evaluated as follows:
  84. First we determine which of the two values (the CDF or its
  85. complement) is likely to be the smaller, the crossover point
  86. is taken to be the mean of the distribution: for this we use the
  87. approximation due to: R. Chattamvelli and R. Shanmugam,
  88. "Algorithm AS 310: Computing the Non-Central Beta Distribution Function",
  89. Applied Statistics, Vol. 46, No. 1. (1997), pp. 146-156.
  90. [equation nc_beta_ref3]
  91. Then either the CDF or its complement is computed using the
  92. relations:
  93. [equation nc_beta_ref4]
  94. The summation is performed by starting at i = [lambda]/2, and then recursing
  95. in both directions, using the usual recurrence relations for the Poisson
  96. PDF and incomplete beta functions. This is the "Method 2" described
  97. by:
  98. Denise Benton and K. Krishnamoorthy,
  99. "Computing discrete mixtures of continuous
  100. distributions: noncentral chisquare, noncentral t
  101. and the distribution of the square of the sample
  102. multiple correlation coefficient",
  103. Computational Statistics & Data Analysis 43 (2003) 249-267.
  104. Specific applications of the above formulae to the noncentral
  105. beta distribution can be found in:
  106. Russell V. Lenth,
  107. "Algorithm AS 226: Computing Noncentral Beta Probabilities",
  108. Applied Statistics, Vol. 36, No. 2. (1987), pp. 241-244.
  109. H. Frick,
  110. "Algorithm AS R84: A Remark on Algorithm AS 226: Computing Non-Central Beta
  111. Probabilities", Applied Statistics, Vol. 39, No. 2. (1990), pp. 311-312.
  112. Ming Long Lam,
  113. "Remark AS R95: A Remark on Algorithm AS 226: Computing Non-Central Beta
  114. Probabilities", Applied Statistics, Vol. 44, No. 4. (1995), pp. 551-552.
  115. Harry O. Posten,
  116. "An Effective Algorithm for the Noncentral Beta Distribution Function",
  117. The American Statistician, Vol. 47, No. 2. (May, 1993), pp. 129-131.
  118. R. Chattamvelli,
  119. "A Note on the Noncentral Beta Distribution Function",
  120. The American Statistician, Vol. 49, No. 2. (May, 1995), pp. 231-234.
  121. Of these, the Posten reference provides the most complete overview,
  122. and includes the modification starting iteration at [lambda]/2.
  123. The main difference between this implementation and the above
  124. references is the direct computation of the complement when most
  125. efficient to do so, and the accumulation of the sum to -1 rather
  126. than subtracting the result from 1 at the end: this can substantially
  127. reduce the number of iterations required when the result is near 1.
  128. The PDF is computed using the methodology of Benton and Krishnamoorthy
  129. and the relation:
  130. [equation nc_beta_ref1]
  131. Quantiles are computed using a specially modified version of
  132. __bracket_solve,
  133. starting the search for the root at the mean of the distribution.
  134. (A Cornish-Fisher type expansion was also tried, but while this gets
  135. quite close to the root in many cases, when it is wrong it tends to
  136. introduce quite pathological behaviour: more investigation in this
  137. area is probably warranted).
  138. [endsect] [/section:nc_beta_dist]
  139. [/ nc_beta.qbk
  140. Copyright 2008 John Maddock and Paul A. Bristow.
  141. Distributed under the Boost Software License, Version 1.0.
  142. (See accompanying file LICENSE_1_0.txt or copy at
  143. http://www.boost.org/LICENSE_1_0.txt).
  144. ]