]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | [section:nc_chi_squared_dist Noncentral Chi-Squared Distribution] |
2 | ||
3 | ``#include <boost/math/distributions/non_central_chi_squared.hpp>`` | |
4 | ||
5 | namespace boost{ namespace math{ | |
6 | ||
7 | template <class RealType = double, | |
8 | class ``__Policy`` = ``__policy_class`` > | |
9 | class non_central_chi_squared_distribution; | |
10 | ||
11 | typedef non_central_chi_squared_distribution<> non_central_chi_squared; | |
12 | ||
13 | template <class RealType, class ``__Policy``> | |
14 | class non_central_chi_squared_distribution | |
15 | { | |
16 | public: | |
17 | typedef RealType value_type; | |
18 | typedef Policy policy_type; | |
19 | ||
20 | // Constructor: | |
21 | non_central_chi_squared_distribution(RealType v, RealType lambda); | |
22 | ||
23 | // Accessor to degrees of freedom parameter v: | |
24 | RealType degrees_of_freedom()const; | |
25 | ||
26 | // Accessor to non centrality parameter lambda: | |
27 | RealType non_centrality()const; | |
28 | ||
29 | // Parameter finders: | |
30 | static RealType find_degrees_of_freedom(RealType lambda, RealType x, RealType p); | |
31 | template <class A, class B, class C> | |
32 | static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c); | |
33 | ||
34 | static RealType find_non_centrality(RealType v, RealType x, RealType p); | |
35 | template <class A, class B, class C> | |
36 | static RealType find_non_centrality(const complemented3_type<A,B,C>& c); | |
37 | }; | |
38 | ||
39 | }} // namespaces | |
40 | ||
41 | The noncentral chi-squared distribution is a generalization of the | |
42 | __chi_squared_distrib. If X[sub i] are [nu] independent, normally | |
43 | distributed random variables with means [mu][sub i] and variances | |
44 | [sigma][sub i][super 2], then the random variable | |
45 | ||
46 | [equation nc_chi_squ_ref1] | |
47 | ||
48 | is distributed according to the noncentral chi-squared distribution. | |
49 | ||
50 | The noncentral chi-squared distribution has two parameters: | |
51 | [nu] which specifies the number of degrees of freedom | |
52 | (i.e. the number of X[sub i]), and [lambda] which is related to the | |
53 | mean of the random variables X[sub i] by: | |
54 | ||
55 | [equation nc_chi_squ_ref2] | |
56 | ||
57 | (Note that some references define [lambda] as one half of the above sum). | |
58 | ||
59 | This leads to a PDF of: | |
60 | ||
61 | [equation nc_chi_squ_ref3] | |
62 | ||
63 | where ['f(x;k)] is the central chi-squared distribution PDF, and | |
64 | ['I[sub v](x)] is a modified Bessel function of the first kind. | |
65 | ||
66 | The following graph illustrates how the distribution changes | |
67 | for different values of [lambda]: | |
68 | ||
69 | [graph nccs_pdf] | |
70 | ||
71 | [h4 Member Functions] | |
72 | ||
73 | non_central_chi_squared_distribution(RealType v, RealType lambda); | |
74 | ||
75 | Constructs a Chi-Squared distribution with /v/ degrees of freedom | |
76 | and non-centrality parameter /lambda/. | |
77 | ||
78 | Requires v > 0 and lambda >= 0, otherwise calls __domain_error. | |
79 | ||
80 | RealType degrees_of_freedom()const; | |
81 | ||
82 | Returns the parameter /v/ from which this object was constructed. | |
83 | ||
84 | RealType non_centrality()const; | |
85 | ||
86 | Returns the parameter /lambda/ from which this object was constructed. | |
87 | ||
88 | ||
89 | static RealType find_degrees_of_freedom(RealType lambda, RealType x, RealType p); | |
90 | ||
91 | This function returns the number of degrees of freedom /v/ such that: | |
92 | `cdf(non_central_chi_squared<RealType, Policy>(v, lambda), x) == p` | |
93 | ||
94 | template <class A, class B, class C> | |
95 | static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c); | |
96 | ||
97 | When called with argument `boost::math::complement(lambda, x, q)` | |
98 | this function returns the number of degrees of freedom /v/ such that: | |
99 | ||
100 | `cdf(complement(non_central_chi_squared<RealType, Policy>(v, lambda), x)) == q`. | |
101 | ||
102 | static RealType find_non_centrality(RealType v, RealType x, RealType p); | |
103 | ||
104 | This function returns the non centrality parameter /lambda/ such that: | |
105 | ||
106 | `cdf(non_central_chi_squared<RealType, Policy>(v, lambda), x) == p` | |
107 | ||
108 | template <class A, class B, class C> | |
109 | static RealType find_non_centrality(const complemented3_type<A,B,C>& c); | |
110 | ||
111 | When called with argument `boost::math::complement(v, x, q)` | |
112 | this function returns the non centrality parameter /lambda/ such that: | |
113 | ||
114 | `cdf(complement(non_central_chi_squared<RealType, Policy>(v, lambda), x)) == q`. | |
115 | ||
116 | [h4 Non-member Accessors] | |
117 | ||
118 | All the [link math_toolkit.dist_ref.nmp usual non-member accessor functions] | |
119 | that are generic to all distributions are supported: __usual_accessors. | |
120 | ||
121 | The domain of the random variable is \[0, +[infin]\]. | |
122 | ||
123 | [h4 Examples] | |
124 | ||
125 | There is a | |
126 | [link math_toolkit.stat_tut.weg.nccs_eg worked example] | |
127 | for the noncentral chi-squared distribution. | |
128 | ||
129 | [h4 Accuracy] | |
130 | ||
131 | The following table shows the peak errors | |
132 | (in units of [@http://en.wikipedia.org/wiki/Machine_epsilon epsilon]) | |
133 | found on various platforms with various floating point types. | |
134 | The failures in the comparison to the [@http://www.r-project.org/ R Math library], | |
135 | seem to be mostly in the corner cases when the probablity would be very small. | |
136 | Unless otherwise specified any floating-point type that is narrower | |
137 | than the one shown will have __zero_error. | |
138 | ||
139 | [table_non_central_chi_squared_CDF] | |
140 | ||
141 | [table_non_central_chi_squared_CDF_complement] | |
142 | ||
143 | Error rates for the quantile | |
144 | functions are broadly similar. Special mention should go to | |
145 | the `mode` function: there is no closed form for this function, | |
146 | so it is evaluated numerically by finding the maxima of the PDF: | |
147 | in principal this can not produce an accuracy greater than the | |
148 | square root of the machine epsilon. | |
149 | ||
150 | [h4 Tests] | |
151 | ||
152 | There are two sets of test data used to verify this implementation: | |
153 | firstly we can compare with published data, for example with | |
154 | Table 6 of "Self-Validating Computations of Probabilities for | |
155 | Selected Central and Noncentral Univariate Probability Functions", | |
156 | Morgan C. Wang and William J. Kennedy, | |
157 | Journal of the American Statistical Association, | |
158 | Vol. 89, No. 427. (Sep., 1994), pp. 878-887. | |
159 | Secondly, we have tables of test data, computed with this | |
160 | implementation and using interval arithmetic - this data should | |
161 | be accurate to at least 50 decimal digits - and is the used for | |
162 | our accuracy tests. | |
163 | ||
164 | [h4 Implementation] | |
165 | ||
166 | The CDF and its complement are evaluated as follows: | |
167 | ||
168 | First we determine which of the two values (the CDF or its | |
169 | complement) is likely to be the smaller: for this we can use the | |
170 | relation due to Temme (see "Asymptotic and Numerical Aspects of the | |
171 | Noncentral Chi-Square Distribution", N. M. Temme, Computers Math. Applic. | |
172 | Vol 25, No. 5, 55-63, 1993) that: | |
173 | ||
174 | F([nu],[lambda];[nu]+[lambda]) [asymp] 0.5 | |
175 | ||
176 | and so compute the CDF when the random variable is less than | |
177 | [nu]+[lambda], and its complement when the random variable is | |
178 | greater than [nu]+[lambda]. If necessary the computed result | |
179 | is then subtracted from 1 to give the desired result (the CDF or its | |
180 | complement). | |
181 | ||
182 | For small values of the non centrality parameter, the CDF is computed | |
183 | using the method of Ding (see "Algorithm AS 275: Computing the Non-Central | |
184 | #2 Distribution Function", Cherng G. Ding, Applied Statistics, Vol. 41, | |
185 | No. 2. (1992), pp. 478-482). This uses the following series representation: | |
186 | ||
187 | [equation nc_chi_squ_ref4] | |
188 | ||
189 | which requires just one call to __gamma_p_derivative with the subsequent | |
190 | terms being computed by recursion as shown above. | |
191 | ||
192 | For larger values of the non-centrality parameter, Ding's method can take | |
193 | an unreasonable number of terms before convergence is achieved. Furthermore, | |
194 | the largest term is not the first term, so in extreme cases the first term may | |
195 | be zero, leading to a zero result, even though the true value may be non-zero. | |
196 | ||
197 | Therefore, when the non-centrality parameter is greater than 200, the method due | |
198 | to Krishnamoorthy (see "Computing discrete mixtures of continuous distributions: | |
199 | noncentral chisquare, noncentral t and the distribution of the | |
200 | square of the sample multiple correlation coefficient", | |
201 | Denise Benton and K. Krishnamoorthy, Computational Statistics & | |
202 | Data Analysis, 43, (2003), 249-267) is used. | |
203 | ||
204 | This method uses the well known sum: | |
205 | ||
206 | [equation nc_chi_squ_ref5] | |
207 | ||
208 | Where P[sub a](x) is the incomplete gamma function. | |
209 | ||
210 | The method starts at the [lambda]th term, which is where the Poisson weighting | |
211 | function achieves its maximum value, although this is not necessarily | |
212 | the largest overall term. Subsequent terms are calculated via the normal | |
213 | recurrence relations for the incomplete gamma function, and iteration proceeds | |
214 | both forwards and backwards until sufficient precision has been achieved. It | |
215 | should be noted that recurrence in the forwards direction of P[sub a](x) is | |
216 | numerically unstable. However, since we always start /after/ the largest | |
217 | term in the series, numeric instability is introduced more slowly than the | |
218 | series converges. | |
219 | ||
220 | Computation of the complement of the CDF uses an extension of Krishnamoorthy's | |
221 | method, given that: | |
222 | ||
223 | [equation nc_chi_squ_ref6] | |
224 | ||
225 | we can again start at the [lambda]'th term and proceed in both directions from | |
226 | there until the required precision is achieved. This time it is backwards | |
227 | recursion on the incomplete gamma function Q[sub a](x) which is unstable. | |
228 | However, as long as we start well /before/ the largest term, this is not an | |
229 | issue in practice. | |
230 | ||
231 | The PDF is computed directly using the relation: | |
232 | ||
233 | [equation nc_chi_squ_ref3] | |
234 | ||
235 | Where ['f(x; v)] is the PDF of the central __chi_squared_distrib and | |
236 | ['I[sub v](x)] is a modified Bessel function, see __cyl_bessel_i. | |
237 | For small values of the | |
238 | non-centrality parameter the relation in terms of __cyl_bessel_i | |
239 | is used. However, this method fails for large values of the | |
240 | non-centrality parameter, so in that case the infinite sum is | |
241 | evaluated using the method of Benton and Krishnamoorthy, and | |
242 | the usual recurrence relations for successive terms. | |
243 | ||
244 | The quantile functions are computed by numeric inversion of the CDF. | |
245 | An improve starting quess is from | |
246 | Thomas Luu, | |
247 | [@http://discovery.ucl.ac.uk/1482128/, Fast and accurate parallel computation of quantile functions for random number generation, Doctorial Thesis, 2016]. | |
248 | ||
249 | There is no [@http://en.wikipedia.org/wiki/Closed_form closed form] | |
250 | for the mode of the noncentral chi-squared | |
251 | distribution: it is computed numerically by finding the maximum | |
252 | of the PDF. Likewise, the median is computed numerically via | |
253 | the quantile. | |
254 | ||
255 | The remaining non-member functions use the following formulas: | |
256 | ||
257 | [equation nc_chi_squ_ref7] | |
258 | ||
259 | Some analytic properties of noncentral distributions | |
260 | (particularly unimodality, and monotonicity of their modes) | |
261 | are surveyed and summarized by: | |
262 | ||
263 | Andrea van Aubel & Wolfgang Gawronski, Applied Mathematics and Computation, 141 (2003) 3-12. | |
264 | ||
265 | [endsect] [/section:nc_chi_squared_dist] | |
266 | ||
267 | [/ nc_chi_squared.qbk | |
268 | Copyright 2008 John Maddock and Paul A. Bristow. | |
269 | Distributed under the Boost Software License, Version 1.0. | |
270 | (See accompanying file LICENSE_1_0.txt or copy at | |
271 | http://www.boost.org/LICENSE_1_0.txt). | |
272 | ] | |
273 |