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