]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // boost\math\distributions\non_central_chi_squared.hpp |
2 | ||
3 | // Copyright John Maddock 2008. | |
4 | ||
5 | // Use, modification and distribution are subject to the | |
6 | // Boost Software License, Version 1.0. | |
7 | // (See accompanying file LICENSE_1_0.txt | |
8 | // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
9 | ||
10 | #ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP | |
11 | #define BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP | |
12 | ||
13 | #include <boost/math/distributions/fwd.hpp> | |
14 | #include <boost/math/special_functions/gamma.hpp> // for incomplete gamma. gamma_q | |
15 | #include <boost/math/special_functions/bessel.hpp> // for cyl_bessel_i | |
16 | #include <boost/math/special_functions/round.hpp> // for iround | |
17 | #include <boost/math/distributions/complement.hpp> // complements | |
18 | #include <boost/math/distributions/chi_squared.hpp> // central distribution | |
19 | #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks | |
20 | #include <boost/math/special_functions/fpclassify.hpp> // isnan. | |
21 | #include <boost/math/tools/roots.hpp> // for root finding. | |
22 | #include <boost/math/distributions/detail/generic_mode.hpp> | |
23 | #include <boost/math/distributions/detail/generic_quantile.hpp> | |
24 | ||
25 | namespace boost | |
26 | { | |
27 | namespace math | |
28 | { | |
29 | ||
30 | template <class RealType, class Policy> | |
31 | class non_central_chi_squared_distribution; | |
32 | ||
33 | namespace detail{ | |
34 | ||
35 | template <class T, class Policy> | |
36 | T non_central_chi_square_q(T x, T f, T theta, const Policy& pol, T init_sum = 0) | |
37 | { | |
38 | // | |
39 | // Computes the complement of the Non-Central Chi-Square | |
40 | // Distribution CDF by summing a weighted sum of complements | |
41 | // of the central-distributions. The weighting factor is | |
42 | // a Poisson Distribution. | |
43 | // | |
44 | // This is an application of the technique described in: | |
45 | // | |
46 | // Computing discrete mixtures of continuous | |
47 | // distributions: noncentral chisquare, noncentral t | |
48 | // and the distribution of the square of the sample | |
49 | // multiple correlation coeficient. | |
50 | // D. Benton, K. Krishnamoorthy. | |
51 | // Computational Statistics & Data Analysis 43 (2003) 249 - 267 | |
52 | // | |
53 | BOOST_MATH_STD_USING | |
54 | ||
55 | // Special case: | |
56 | if(x == 0) | |
57 | return 1; | |
58 | ||
59 | // | |
60 | // Initialize the variables we'll be using: | |
61 | // | |
62 | T lambda = theta / 2; | |
63 | T del = f / 2; | |
64 | T y = x / 2; | |
65 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
66 | T errtol = boost::math::policies::get_epsilon<T, Policy>(); | |
67 | T sum = init_sum; | |
68 | // | |
69 | // k is the starting location for iteration, we'll | |
70 | // move both forwards and backwards from this point. | |
71 | // k is chosen as the peek of the Poisson weights, which | |
72 | // will occur *before* the largest term. | |
73 | // | |
74 | int k = iround(lambda, pol); | |
75 | // Forwards and backwards Poisson weights: | |
76 | T poisf = boost::math::gamma_p_derivative(static_cast<T>(1 + k), lambda, pol); | |
77 | T poisb = poisf * k / lambda; | |
78 | // Initial forwards central chi squared term: | |
79 | T gamf = boost::math::gamma_q(del + k, y, pol); | |
80 | // Forwards and backwards recursion terms on the central chi squared: | |
81 | T xtermf = boost::math::gamma_p_derivative(del + 1 + k, y, pol); | |
82 | T xtermb = xtermf * (del + k) / y; | |
83 | // Initial backwards central chi squared term: | |
84 | T gamb = gamf - xtermb; | |
85 | ||
86 | // | |
87 | // Forwards iteration first, this is the | |
88 | // stable direction for the gamma function | |
89 | // recurrences: | |
90 | // | |
91 | int i; | |
92 | for(i = k; static_cast<boost::uintmax_t>(i-k) < max_iter; ++i) | |
93 | { | |
94 | T term = poisf * gamf; | |
95 | sum += term; | |
96 | poisf *= lambda / (i + 1); | |
97 | gamf += xtermf; | |
98 | xtermf *= y / (del + i + 1); | |
99 | if(((sum == 0) || (fabs(term / sum) < errtol)) && (term >= poisf * gamf)) | |
100 | break; | |
101 | } | |
102 | //Error check: | |
103 | if(static_cast<boost::uintmax_t>(i-k) >= max_iter) | |
104 | return policies::raise_evaluation_error( | |
105 | "cdf(non_central_chi_squared_distribution<%1%>, %1%)", | |
106 | "Series did not converge, closest value was %1%", sum, pol); | |
107 | // | |
108 | // Now backwards iteration: the gamma | |
109 | // function recurrences are unstable in this | |
110 | // direction, we rely on the terms deminishing in size | |
111 | // faster than we introduce cancellation errors. | |
112 | // For this reason it's very important that we start | |
113 | // *before* the largest term so that backwards iteration | |
114 | // is strictly converging. | |
115 | // | |
116 | for(i = k - 1; i >= 0; --i) | |
117 | { | |
118 | T term = poisb * gamb; | |
119 | sum += term; | |
120 | poisb *= i / lambda; | |
121 | xtermb *= (del + i) / y; | |
122 | gamb -= xtermb; | |
123 | if((sum == 0) || (fabs(term / sum) < errtol)) | |
124 | break; | |
125 | } | |
126 | ||
127 | return sum; | |
128 | } | |
129 | ||
130 | template <class T, class Policy> | |
131 | T non_central_chi_square_p_ding(T x, T f, T theta, const Policy& pol, T init_sum = 0) | |
132 | { | |
133 | // | |
134 | // This is an implementation of: | |
135 | // | |
136 | // Algorithm AS 275: | |
137 | // Computing the Non-Central #2 Distribution Function | |
138 | // Cherng G. Ding | |
139 | // Applied Statistics, Vol. 41, No. 2. (1992), pp. 478-482. | |
140 | // | |
141 | // This uses a stable forward iteration to sum the | |
142 | // CDF, unfortunately this can not be used for large | |
143 | // values of the non-centrality parameter because: | |
144 | // * The first term may underfow to zero. | |
145 | // * We may need an extra-ordinary number of terms | |
146 | // before we reach the first *significant* term. | |
147 | // | |
148 | BOOST_MATH_STD_USING | |
149 | // Special case: | |
150 | if(x == 0) | |
151 | return 0; | |
152 | T tk = boost::math::gamma_p_derivative(f/2 + 1, x/2, pol); | |
153 | T lambda = theta / 2; | |
154 | T vk = exp(-lambda); | |
155 | T uk = vk; | |
156 | T sum = init_sum + tk * vk; | |
157 | if(sum == 0) | |
158 | return sum; | |
159 | ||
160 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
161 | T errtol = boost::math::policies::get_epsilon<T, Policy>(); | |
162 | ||
163 | int i; | |
164 | T lterm(0), term(0); | |
165 | for(i = 1; static_cast<boost::uintmax_t>(i) < max_iter; ++i) | |
166 | { | |
167 | tk = tk * x / (f + 2 * i); | |
168 | uk = uk * lambda / i; | |
169 | vk = vk + uk; | |
170 | lterm = term; | |
171 | term = vk * tk; | |
172 | sum += term; | |
173 | if((fabs(term / sum) < errtol) && (term <= lterm)) | |
174 | break; | |
175 | } | |
176 | //Error check: | |
177 | if(static_cast<boost::uintmax_t>(i) >= max_iter) | |
178 | return policies::raise_evaluation_error( | |
179 | "cdf(non_central_chi_squared_distribution<%1%>, %1%)", | |
180 | "Series did not converge, closest value was %1%", sum, pol); | |
181 | return sum; | |
182 | } | |
183 | ||
184 | ||
185 | template <class T, class Policy> | |
186 | T non_central_chi_square_p(T y, T n, T lambda, const Policy& pol, T init_sum) | |
187 | { | |
188 | // | |
189 | // This is taken more or less directly from: | |
190 | // | |
191 | // Computing discrete mixtures of continuous | |
192 | // distributions: noncentral chisquare, noncentral t | |
193 | // and the distribution of the square of the sample | |
194 | // multiple correlation coeficient. | |
195 | // D. Benton, K. Krishnamoorthy. | |
196 | // Computational Statistics & Data Analysis 43 (2003) 249 - 267 | |
197 | // | |
198 | // We're summing a Poisson weighting term multiplied by | |
199 | // a central chi squared distribution. | |
200 | // | |
201 | BOOST_MATH_STD_USING | |
202 | // Special case: | |
203 | if(y == 0) | |
204 | return 0; | |
205 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
206 | T errtol = boost::math::policies::get_epsilon<T, Policy>(); | |
207 | T errorf(0), errorb(0); | |
208 | ||
209 | T x = y / 2; | |
210 | T del = lambda / 2; | |
211 | // | |
212 | // Starting location for the iteration, we'll iterate | |
213 | // both forwards and backwards from this point. The | |
214 | // location chosen is the maximum of the Poisson weight | |
215 | // function, which ocurrs *after* the largest term in the | |
216 | // sum. | |
217 | // | |
218 | int k = iround(del, pol); | |
219 | T a = n / 2 + k; | |
220 | // Central chi squared term for forward iteration: | |
221 | T gamkf = boost::math::gamma_p(a, x, pol); | |
222 | ||
223 | if(lambda == 0) | |
224 | return gamkf; | |
225 | // Central chi squared term for backward iteration: | |
226 | T gamkb = gamkf; | |
227 | // Forwards Poisson weight: | |
228 | T poiskf = gamma_p_derivative(static_cast<T>(k+1), del, pol); | |
229 | // Backwards Poisson weight: | |
230 | T poiskb = poiskf; | |
231 | // Forwards gamma function recursion term: | |
232 | T xtermf = boost::math::gamma_p_derivative(a, x, pol); | |
233 | // Backwards gamma function recursion term: | |
234 | T xtermb = xtermf * x / a; | |
235 | T sum = init_sum + poiskf * gamkf; | |
236 | if(sum == 0) | |
237 | return sum; | |
238 | int i = 1; | |
239 | // | |
240 | // Backwards recursion first, this is the stable | |
241 | // direction for gamma function recurrences: | |
242 | // | |
243 | while(i <= k) | |
244 | { | |
245 | xtermb *= (a - i + 1) / x; | |
246 | gamkb += xtermb; | |
247 | poiskb = poiskb * (k - i + 1) / del; | |
248 | errorf = errorb; | |
249 | errorb = gamkb * poiskb; | |
250 | sum += errorb; | |
251 | if((fabs(errorb / sum) < errtol) && (errorb <= errorf)) | |
252 | break; | |
253 | ++i; | |
254 | } | |
255 | i = 1; | |
256 | // | |
257 | // Now forwards recursion, the gamma function | |
258 | // recurrence relation is unstable in this direction, | |
259 | // so we rely on the magnitude of successive terms | |
260 | // decreasing faster than we introduce cancellation error. | |
261 | // For this reason it's vital that k is chosen to be *after* | |
262 | // the largest term, so that successive forward iterations | |
263 | // are strictly (and rapidly) converging. | |
264 | // | |
265 | do | |
266 | { | |
267 | xtermf = xtermf * x / (a + i - 1); | |
268 | gamkf = gamkf - xtermf; | |
269 | poiskf = poiskf * del / (k + i); | |
270 | errorf = poiskf * gamkf; | |
271 | sum += errorf; | |
272 | ++i; | |
273 | }while((fabs(errorf / sum) > errtol) && (static_cast<boost::uintmax_t>(i) < max_iter)); | |
274 | ||
275 | //Error check: | |
276 | if(static_cast<boost::uintmax_t>(i) >= max_iter) | |
277 | return policies::raise_evaluation_error( | |
278 | "cdf(non_central_chi_squared_distribution<%1%>, %1%)", | |
279 | "Series did not converge, closest value was %1%", sum, pol); | |
280 | ||
281 | return sum; | |
282 | } | |
283 | ||
284 | template <class T, class Policy> | |
285 | T non_central_chi_square_pdf(T x, T n, T lambda, const Policy& pol) | |
286 | { | |
287 | // | |
288 | // As above but for the PDF: | |
289 | // | |
290 | BOOST_MATH_STD_USING | |
291 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
292 | T errtol = boost::math::policies::get_epsilon<T, Policy>(); | |
293 | T x2 = x / 2; | |
294 | T n2 = n / 2; | |
295 | T l2 = lambda / 2; | |
296 | T sum = 0; | |
297 | int k = itrunc(l2); | |
298 | T pois = gamma_p_derivative(static_cast<T>(k + 1), l2, pol) * gamma_p_derivative(static_cast<T>(n2 + k), x2); | |
299 | if(pois == 0) | |
300 | return 0; | |
301 | T poisb = pois; | |
302 | for(int i = k; ; ++i) | |
303 | { | |
304 | sum += pois; | |
305 | if(pois / sum < errtol) | |
306 | break; | |
307 | if(static_cast<boost::uintmax_t>(i - k) >= max_iter) | |
308 | return policies::raise_evaluation_error( | |
309 | "pdf(non_central_chi_squared_distribution<%1%>, %1%)", | |
310 | "Series did not converge, closest value was %1%", sum, pol); | |
311 | pois *= l2 * x2 / ((i + 1) * (n2 + i)); | |
312 | } | |
313 | for(int i = k - 1; i >= 0; --i) | |
314 | { | |
315 | poisb *= (i + 1) * (n2 + i) / (l2 * x2); | |
316 | sum += poisb; | |
317 | if(poisb / sum < errtol) | |
318 | break; | |
319 | } | |
320 | return sum / 2; | |
321 | } | |
322 | ||
323 | template <class RealType, class Policy> | |
324 | inline RealType non_central_chi_squared_cdf(RealType x, RealType k, RealType l, bool invert, const Policy&) | |
325 | { | |
326 | typedef typename policies::evaluation<RealType, Policy>::type value_type; | |
327 | typedef typename policies::normalise< | |
328 | Policy, | |
329 | policies::promote_float<false>, | |
330 | policies::promote_double<false>, | |
331 | policies::discrete_quantile<>, | |
332 | policies::assert_undefined<> >::type forwarding_policy; | |
333 | ||
334 | BOOST_MATH_STD_USING | |
335 | value_type result; | |
336 | if(l == 0) | |
337 | return invert == false ? cdf(boost::math::chi_squared_distribution<RealType, Policy>(k), x) : cdf(complement(boost::math::chi_squared_distribution<RealType, Policy>(k), x)); | |
338 | else if(x > k + l) | |
339 | { | |
340 | // Complement is the smaller of the two: | |
341 | result = detail::non_central_chi_square_q( | |
342 | static_cast<value_type>(x), | |
343 | static_cast<value_type>(k), | |
344 | static_cast<value_type>(l), | |
345 | forwarding_policy(), | |
346 | static_cast<value_type>(invert ? 0 : -1)); | |
347 | invert = !invert; | |
348 | } | |
349 | else if(l < 200) | |
350 | { | |
351 | // For small values of the non-centrality parameter | |
352 | // we can use Ding's method: | |
353 | result = detail::non_central_chi_square_p_ding( | |
354 | static_cast<value_type>(x), | |
355 | static_cast<value_type>(k), | |
356 | static_cast<value_type>(l), | |
357 | forwarding_policy(), | |
358 | static_cast<value_type>(invert ? -1 : 0)); | |
359 | } | |
360 | else | |
361 | { | |
362 | // For largers values of the non-centrality | |
363 | // parameter Ding's method will consume an | |
364 | // extra-ordinary number of terms, and worse | |
365 | // may return zero when the result is in fact | |
366 | // finite, use Krishnamoorthy's method instead: | |
367 | result = detail::non_central_chi_square_p( | |
368 | static_cast<value_type>(x), | |
369 | static_cast<value_type>(k), | |
370 | static_cast<value_type>(l), | |
371 | forwarding_policy(), | |
372 | static_cast<value_type>(invert ? -1 : 0)); | |
373 | } | |
374 | if(invert) | |
375 | result = -result; | |
376 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
377 | result, | |
378 | "boost::math::non_central_chi_squared_cdf<%1%>(%1%, %1%, %1%)"); | |
379 | } | |
380 | ||
381 | template <class T, class Policy> | |
382 | struct nccs_quantile_functor | |
383 | { | |
384 | nccs_quantile_functor(const non_central_chi_squared_distribution<T,Policy>& d, T t, bool c) | |
385 | : dist(d), target(t), comp(c) {} | |
386 | ||
387 | T operator()(const T& x) | |
388 | { | |
389 | return comp ? | |
390 | target - cdf(complement(dist, x)) | |
391 | : cdf(dist, x) - target; | |
392 | } | |
393 | ||
394 | private: | |
395 | non_central_chi_squared_distribution<T,Policy> dist; | |
396 | T target; | |
397 | bool comp; | |
398 | }; | |
399 | ||
400 | template <class RealType, class Policy> | |
401 | RealType nccs_quantile(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& p, bool comp) | |
402 | { | |
403 | BOOST_MATH_STD_USING | |
404 | static const char* function = "quantile(non_central_chi_squared_distribution<%1%>, %1%)"; | |
405 | typedef typename policies::evaluation<RealType, Policy>::type value_type; | |
406 | typedef typename policies::normalise< | |
407 | Policy, | |
408 | policies::promote_float<false>, | |
409 | policies::promote_double<false>, | |
410 | policies::discrete_quantile<>, | |
411 | policies::assert_undefined<> >::type forwarding_policy; | |
412 | ||
413 | value_type k = dist.degrees_of_freedom(); | |
414 | value_type l = dist.non_centrality(); | |
415 | value_type r; | |
416 | if(!detail::check_df( | |
417 | function, | |
418 | k, &r, Policy()) | |
419 | || | |
420 | !detail::check_non_centrality( | |
421 | function, | |
422 | l, | |
423 | &r, | |
424 | Policy()) | |
425 | || | |
426 | !detail::check_probability( | |
427 | function, | |
428 | static_cast<value_type>(p), | |
429 | &r, | |
430 | Policy())) | |
431 | return (RealType)r; | |
432 | // | |
433 | // Special cases get short-circuited first: | |
434 | // | |
435 | if(p == 0) | |
436 | return comp ? policies::raise_overflow_error<RealType>(function, 0, Policy()) : 0; | |
437 | if(p == 1) | |
438 | return comp ? 0 : policies::raise_overflow_error<RealType>(function, 0, Policy()); | |
439 | // | |
440 | // This is Pearson's approximation to the quantile, see | |
441 | // Pearson, E. S. (1959) "Note on an approximation to the distribution of | |
442 | // noncentral chi squared", Biometrika 46: 364. | |
443 | // See also: | |
444 | // "A comparison of approximations to percentiles of the noncentral chi2-distribution", | |
445 | // Hardeo Sahai and Mario Miguel Ojeda, Revista de Matematica: Teoria y Aplicaciones 2003 10(1-2) : 57-76. | |
446 | // Note that the latter reference refers to an approximation of the CDF, when they really mean the quantile. | |
447 | // | |
448 | value_type b = -(l * l) / (k + 3 * l); | |
449 | value_type c = (k + 3 * l) / (k + 2 * l); | |
450 | value_type ff = (k + 2 * l) / (c * c); | |
451 | value_type guess; | |
452 | if(comp) | |
453 | { | |
454 | guess = b + c * quantile(complement(chi_squared_distribution<value_type, forwarding_policy>(ff), p)); | |
455 | } | |
456 | else | |
457 | { | |
458 | guess = b + c * quantile(chi_squared_distribution<value_type, forwarding_policy>(ff), p); | |
459 | } | |
460 | // | |
461 | // Sometimes guess goes very small or negative, in that case we have | |
462 | // to do something else for the initial guess, this approximation | |
463 | // was provided in a private communication from Thomas Luu, PhD candidate, | |
464 | // University College London. It's an asymptotic expansion for the | |
465 | // quantile which usually gets us within an order of magnitude of the | |
466 | // correct answer. | |
467 | // Fast and accurate parallel computation of quantile functions for random number generation, | |
468 | // Thomas LuuDoctorial Thesis 2016 | |
469 | // http://discovery.ucl.ac.uk/1482128/ | |
470 | // | |
471 | if(guess < 0.005) | |
472 | { | |
473 | value_type pp = comp ? 1 - p : p; | |
474 | //guess = pow(pow(value_type(2), (k / 2 - 1)) * exp(l / 2) * pp * k, 2 / k); | |
475 | guess = pow(pow(value_type(2), (k / 2 - 1)) * exp(l / 2) * pp * k * boost::math::tgamma(k / 2, forwarding_policy()), (2 / k)); | |
476 | if(guess == 0) | |
477 | guess = tools::min_value<value_type>(); | |
478 | } | |
479 | value_type result = detail::generic_quantile( | |
480 | non_central_chi_squared_distribution<value_type, forwarding_policy>(k, l), | |
481 | p, | |
482 | guess, | |
483 | comp, | |
484 | function); | |
485 | ||
486 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
487 | result, | |
488 | function); | |
489 | } | |
490 | ||
491 | template <class RealType, class Policy> | |
492 | RealType nccs_pdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x) | |
493 | { | |
494 | BOOST_MATH_STD_USING | |
495 | static const char* function = "pdf(non_central_chi_squared_distribution<%1%>, %1%)"; | |
496 | typedef typename policies::evaluation<RealType, Policy>::type value_type; | |
497 | typedef typename policies::normalise< | |
498 | Policy, | |
499 | policies::promote_float<false>, | |
500 | policies::promote_double<false>, | |
501 | policies::discrete_quantile<>, | |
502 | policies::assert_undefined<> >::type forwarding_policy; | |
503 | ||
504 | value_type k = dist.degrees_of_freedom(); | |
505 | value_type l = dist.non_centrality(); | |
506 | value_type r; | |
507 | if(!detail::check_df( | |
508 | function, | |
509 | k, &r, Policy()) | |
510 | || | |
511 | !detail::check_non_centrality( | |
512 | function, | |
513 | l, | |
514 | &r, | |
515 | Policy()) | |
516 | || | |
517 | !detail::check_positive_x( | |
518 | function, | |
519 | (value_type)x, | |
520 | &r, | |
521 | Policy())) | |
522 | return (RealType)r; | |
523 | ||
524 | if(l == 0) | |
525 | return pdf(boost::math::chi_squared_distribution<RealType, forwarding_policy>(dist.degrees_of_freedom()), x); | |
526 | ||
527 | // Special case: | |
528 | if(x == 0) | |
529 | return 0; | |
530 | if(l > 50) | |
531 | { | |
532 | r = non_central_chi_square_pdf(static_cast<value_type>(x), k, l, forwarding_policy()); | |
533 | } | |
534 | else | |
535 | { | |
536 | r = log(x / l) * (k / 4 - 0.5f) - (x + l) / 2; | |
537 | if(fabs(r) >= tools::log_max_value<RealType>() / 4) | |
538 | { | |
539 | r = non_central_chi_square_pdf(static_cast<value_type>(x), k, l, forwarding_policy()); | |
540 | } | |
541 | else | |
542 | { | |
543 | r = exp(r); | |
544 | r = 0.5f * r | |
545 | * boost::math::cyl_bessel_i(k/2 - 1, sqrt(l * x), forwarding_policy()); | |
546 | } | |
547 | } | |
548 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
549 | r, | |
550 | function); | |
551 | } | |
552 | ||
553 | template <class RealType, class Policy> | |
554 | struct degrees_of_freedom_finder | |
555 | { | |
556 | degrees_of_freedom_finder( | |
557 | RealType lam_, RealType x_, RealType p_, bool c) | |
558 | : lam(lam_), x(x_), p(p_), comp(c) {} | |
559 | ||
560 | RealType operator()(const RealType& v) | |
561 | { | |
562 | non_central_chi_squared_distribution<RealType, Policy> d(v, lam); | |
563 | return comp ? | |
564 | RealType(p - cdf(complement(d, x))) | |
565 | : RealType(cdf(d, x) - p); | |
566 | } | |
567 | private: | |
568 | RealType lam; | |
569 | RealType x; | |
570 | RealType p; | |
571 | bool comp; | |
572 | }; | |
573 | ||
574 | template <class RealType, class Policy> | |
575 | inline RealType find_degrees_of_freedom( | |
576 | RealType lam, RealType x, RealType p, RealType q, const Policy& pol) | |
577 | { | |
578 | const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom"; | |
579 | if((p == 0) || (q == 0)) | |
580 | { | |
581 | // | |
582 | // Can't a thing if one of p and q is zero: | |
583 | // | |
584 | return policies::raise_evaluation_error<RealType>(function, | |
585 | "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%", | |
586 | RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); | |
587 | } | |
588 | degrees_of_freedom_finder<RealType, Policy> f(lam, x, p < q ? p : q, p < q ? false : true); | |
589 | tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>()); | |
590 | boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); | |
591 | // | |
592 | // Pick an initial guess that we know will give us a probability | |
593 | // right around 0.5. | |
594 | // | |
595 | RealType guess = x - lam; | |
596 | if(guess < 1) | |
597 | guess = 1; | |
598 | std::pair<RealType, RealType> ir = tools::bracket_and_solve_root( | |
599 | f, guess, RealType(2), false, tol, max_iter, pol); | |
600 | RealType result = ir.first + (ir.second - ir.first) / 2; | |
601 | if(max_iter >= policies::get_max_root_iterations<Policy>()) | |
602 | { | |
603 | return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" | |
604 | " or there is no answer to problem. Current best guess is %1%", result, Policy()); | |
605 | } | |
606 | return result; | |
607 | } | |
608 | ||
609 | template <class RealType, class Policy> | |
610 | struct non_centrality_finder | |
611 | { | |
612 | non_centrality_finder( | |
613 | RealType v_, RealType x_, RealType p_, bool c) | |
614 | : v(v_), x(x_), p(p_), comp(c) {} | |
615 | ||
616 | RealType operator()(const RealType& lam) | |
617 | { | |
618 | non_central_chi_squared_distribution<RealType, Policy> d(v, lam); | |
619 | return comp ? | |
620 | RealType(p - cdf(complement(d, x))) | |
621 | : RealType(cdf(d, x) - p); | |
622 | } | |
623 | private: | |
624 | RealType v; | |
625 | RealType x; | |
626 | RealType p; | |
627 | bool comp; | |
628 | }; | |
629 | ||
630 | template <class RealType, class Policy> | |
631 | inline RealType find_non_centrality( | |
632 | RealType v, RealType x, RealType p, RealType q, const Policy& pol) | |
633 | { | |
634 | const char* function = "non_central_chi_squared<%1%>::find_non_centrality"; | |
635 | if((p == 0) || (q == 0)) | |
636 | { | |
637 | // | |
638 | // Can't do a thing if one of p and q is zero: | |
639 | // | |
640 | return policies::raise_evaluation_error<RealType>(function, | |
641 | "Can't find non centrality parameter when the probability is 0 or 1, only possible answer is %1%", | |
642 | RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); | |
643 | } | |
644 | non_centrality_finder<RealType, Policy> f(v, x, p < q ? p : q, p < q ? false : true); | |
645 | tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>()); | |
646 | boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); | |
647 | // | |
648 | // Pick an initial guess that we know will give us a probability | |
649 | // right around 0.5. | |
650 | // | |
651 | RealType guess = x - v; | |
652 | if(guess < 1) | |
653 | guess = 1; | |
654 | std::pair<RealType, RealType> ir = tools::bracket_and_solve_root( | |
655 | f, guess, RealType(2), false, tol, max_iter, pol); | |
656 | RealType result = ir.first + (ir.second - ir.first) / 2; | |
657 | if(max_iter >= policies::get_max_root_iterations<Policy>()) | |
658 | { | |
659 | return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" | |
660 | " or there is no answer to problem. Current best guess is %1%", result, Policy()); | |
661 | } | |
662 | return result; | |
663 | } | |
664 | ||
665 | } | |
666 | ||
667 | template <class RealType = double, class Policy = policies::policy<> > | |
668 | class non_central_chi_squared_distribution | |
669 | { | |
670 | public: | |
671 | typedef RealType value_type; | |
672 | typedef Policy policy_type; | |
673 | ||
674 | non_central_chi_squared_distribution(RealType df_, RealType lambda) : df(df_), ncp(lambda) | |
675 | { | |
676 | const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::non_central_chi_squared_distribution(%1%,%1%)"; | |
677 | RealType r; | |
678 | detail::check_df( | |
679 | function, | |
680 | df, &r, Policy()); | |
681 | detail::check_non_centrality( | |
682 | function, | |
683 | ncp, | |
684 | &r, | |
685 | Policy()); | |
686 | } // non_central_chi_squared_distribution constructor. | |
687 | ||
688 | RealType degrees_of_freedom() const | |
689 | { // Private data getter function. | |
690 | return df; | |
691 | } | |
692 | RealType non_centrality() const | |
693 | { // Private data getter function. | |
694 | return ncp; | |
695 | } | |
696 | static RealType find_degrees_of_freedom(RealType lam, RealType x, RealType p) | |
697 | { | |
698 | const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom"; | |
699 | typedef typename policies::evaluation<RealType, Policy>::type eval_type; | |
700 | typedef typename policies::normalise< | |
701 | Policy, | |
702 | policies::promote_float<false>, | |
703 | policies::promote_double<false>, | |
704 | policies::discrete_quantile<>, | |
705 | policies::assert_undefined<> >::type forwarding_policy; | |
706 | eval_type result = detail::find_degrees_of_freedom( | |
707 | static_cast<eval_type>(lam), | |
708 | static_cast<eval_type>(x), | |
709 | static_cast<eval_type>(p), | |
710 | static_cast<eval_type>(1-p), | |
711 | forwarding_policy()); | |
712 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
713 | result, | |
714 | function); | |
715 | } | |
716 | template <class A, class B, class C> | |
717 | static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c) | |
718 | { | |
719 | const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom"; | |
720 | typedef typename policies::evaluation<RealType, Policy>::type eval_type; | |
721 | typedef typename policies::normalise< | |
722 | Policy, | |
723 | policies::promote_float<false>, | |
724 | policies::promote_double<false>, | |
725 | policies::discrete_quantile<>, | |
726 | policies::assert_undefined<> >::type forwarding_policy; | |
727 | eval_type result = detail::find_degrees_of_freedom( | |
728 | static_cast<eval_type>(c.dist), | |
729 | static_cast<eval_type>(c.param1), | |
730 | static_cast<eval_type>(1-c.param2), | |
731 | static_cast<eval_type>(c.param2), | |
732 | forwarding_policy()); | |
733 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
734 | result, | |
735 | function); | |
736 | } | |
737 | static RealType find_non_centrality(RealType v, RealType x, RealType p) | |
738 | { | |
739 | const char* function = "non_central_chi_squared<%1%>::find_non_centrality"; | |
740 | typedef typename policies::evaluation<RealType, Policy>::type eval_type; | |
741 | typedef typename policies::normalise< | |
742 | Policy, | |
743 | policies::promote_float<false>, | |
744 | policies::promote_double<false>, | |
745 | policies::discrete_quantile<>, | |
746 | policies::assert_undefined<> >::type forwarding_policy; | |
747 | eval_type result = detail::find_non_centrality( | |
748 | static_cast<eval_type>(v), | |
749 | static_cast<eval_type>(x), | |
750 | static_cast<eval_type>(p), | |
751 | static_cast<eval_type>(1-p), | |
752 | forwarding_policy()); | |
753 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
754 | result, | |
755 | function); | |
756 | } | |
757 | template <class A, class B, class C> | |
758 | static RealType find_non_centrality(const complemented3_type<A,B,C>& c) | |
759 | { | |
760 | const char* function = "non_central_chi_squared<%1%>::find_non_centrality"; | |
761 | typedef typename policies::evaluation<RealType, Policy>::type eval_type; | |
762 | typedef typename policies::normalise< | |
763 | Policy, | |
764 | policies::promote_float<false>, | |
765 | policies::promote_double<false>, | |
766 | policies::discrete_quantile<>, | |
767 | policies::assert_undefined<> >::type forwarding_policy; | |
768 | eval_type result = detail::find_non_centrality( | |
769 | static_cast<eval_type>(c.dist), | |
770 | static_cast<eval_type>(c.param1), | |
771 | static_cast<eval_type>(1-c.param2), | |
772 | static_cast<eval_type>(c.param2), | |
773 | forwarding_policy()); | |
774 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
775 | result, | |
776 | function); | |
777 | } | |
778 | private: | |
779 | // Data member, initialized by constructor. | |
780 | RealType df; // degrees of freedom. | |
781 | RealType ncp; // non-centrality parameter | |
782 | }; // template <class RealType, class Policy> class non_central_chi_squared_distribution | |
783 | ||
784 | typedef non_central_chi_squared_distribution<double> non_central_chi_squared; // Reserved name of type double. | |
785 | ||
786 | // Non-member functions to give properties of the distribution. | |
787 | ||
788 | template <class RealType, class Policy> | |
789 | inline const std::pair<RealType, RealType> range(const non_central_chi_squared_distribution<RealType, Policy>& /* dist */) | |
790 | { // Range of permissible values for random variable k. | |
791 | using boost::math::tools::max_value; | |
792 | return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // Max integer? | |
793 | } | |
794 | ||
795 | template <class RealType, class Policy> | |
796 | inline const std::pair<RealType, RealType> support(const non_central_chi_squared_distribution<RealType, Policy>& /* dist */) | |
797 | { // Range of supported values for random variable k. | |
798 | // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. | |
799 | using boost::math::tools::max_value; | |
800 | return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); | |
801 | } | |
802 | ||
803 | template <class RealType, class Policy> | |
804 | inline RealType mean(const non_central_chi_squared_distribution<RealType, Policy>& dist) | |
805 | { // Mean of poisson distribution = lambda. | |
806 | const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::mean()"; | |
807 | RealType k = dist.degrees_of_freedom(); | |
808 | RealType l = dist.non_centrality(); | |
809 | RealType r; | |
810 | if(!detail::check_df( | |
811 | function, | |
812 | k, &r, Policy()) | |
813 | || | |
814 | !detail::check_non_centrality( | |
815 | function, | |
816 | l, | |
817 | &r, | |
818 | Policy())) | |
819 | return r; | |
820 | return k + l; | |
821 | } // mean | |
822 | ||
823 | template <class RealType, class Policy> | |
824 | inline RealType mode(const non_central_chi_squared_distribution<RealType, Policy>& dist) | |
825 | { // mode. | |
826 | static const char* function = "mode(non_central_chi_squared_distribution<%1%> const&)"; | |
827 | ||
828 | RealType k = dist.degrees_of_freedom(); | |
829 | RealType l = dist.non_centrality(); | |
830 | RealType r; | |
831 | if(!detail::check_df( | |
832 | function, | |
833 | k, &r, Policy()) | |
834 | || | |
835 | !detail::check_non_centrality( | |
836 | function, | |
837 | l, | |
838 | &r, | |
839 | Policy())) | |
840 | return (RealType)r; | |
841 | return detail::generic_find_mode(dist, 1 + k, function); | |
842 | } | |
843 | ||
844 | template <class RealType, class Policy> | |
845 | inline RealType variance(const non_central_chi_squared_distribution<RealType, Policy>& dist) | |
846 | { // variance. | |
847 | const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::variance()"; | |
848 | RealType k = dist.degrees_of_freedom(); | |
849 | RealType l = dist.non_centrality(); | |
850 | RealType r; | |
851 | if(!detail::check_df( | |
852 | function, | |
853 | k, &r, Policy()) | |
854 | || | |
855 | !detail::check_non_centrality( | |
856 | function, | |
857 | l, | |
858 | &r, | |
859 | Policy())) | |
860 | return r; | |
861 | return 2 * (2 * l + k); | |
862 | } | |
863 | ||
864 | // RealType standard_deviation(const non_central_chi_squared_distribution<RealType, Policy>& dist) | |
865 | // standard_deviation provided by derived accessors. | |
866 | ||
867 | template <class RealType, class Policy> | |
868 | inline RealType skewness(const non_central_chi_squared_distribution<RealType, Policy>& dist) | |
869 | { // skewness = sqrt(l). | |
870 | const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::skewness()"; | |
871 | RealType k = dist.degrees_of_freedom(); | |
872 | RealType l = dist.non_centrality(); | |
873 | RealType r; | |
874 | if(!detail::check_df( | |
875 | function, | |
876 | k, &r, Policy()) | |
877 | || | |
878 | !detail::check_non_centrality( | |
879 | function, | |
880 | l, | |
881 | &r, | |
882 | Policy())) | |
883 | return r; | |
884 | BOOST_MATH_STD_USING | |
885 | return pow(2 / (k + 2 * l), RealType(3)/2) * (k + 3 * l); | |
886 | } | |
887 | ||
888 | template <class RealType, class Policy> | |
889 | inline RealType kurtosis_excess(const non_central_chi_squared_distribution<RealType, Policy>& dist) | |
890 | { | |
891 | const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::kurtosis_excess()"; | |
892 | RealType k = dist.degrees_of_freedom(); | |
893 | RealType l = dist.non_centrality(); | |
894 | RealType r; | |
895 | if(!detail::check_df( | |
896 | function, | |
897 | k, &r, Policy()) | |
898 | || | |
899 | !detail::check_non_centrality( | |
900 | function, | |
901 | l, | |
902 | &r, | |
903 | Policy())) | |
904 | return r; | |
905 | return 12 * (k + 4 * l) / ((k + 2 * l) * (k + 2 * l)); | |
906 | } // kurtosis_excess | |
907 | ||
908 | template <class RealType, class Policy> | |
909 | inline RealType kurtosis(const non_central_chi_squared_distribution<RealType, Policy>& dist) | |
910 | { | |
911 | return kurtosis_excess(dist) + 3; | |
912 | } | |
913 | ||
914 | template <class RealType, class Policy> | |
915 | inline RealType pdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x) | |
916 | { // Probability Density/Mass Function. | |
917 | return detail::nccs_pdf(dist, x); | |
918 | ||
919 | ||
920 | template <class RealType, class Policy> | |
921 | RealType cdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x) | |
922 | { | |
923 | const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::cdf(%1%)"; | |
924 | RealType k = dist.degrees_of_freedom(); | |
925 | RealType l = dist.non_centrality(); | |
926 | RealType r; | |
927 | if(!detail::check_df( | |
928 | function, | |
929 | k, &r, Policy()) | |
930 | || | |
931 | !detail::check_non_centrality( | |
932 | function, | |
933 | l, | |
934 | &r, | |
935 | Policy()) | |
936 | || | |
937 | !detail::check_positive_x( | |
938 | function, | |
939 | x, | |
940 | &r, | |
941 | Policy())) | |
942 | return r; | |
943 | ||
944 | return detail::non_central_chi_squared_cdf(x, k, l, false, Policy()); | |
945 | } // cdf | |
946 | ||
947 | template <class RealType, class Policy> | |
948 | RealType cdf(const complemented2_type<non_central_chi_squared_distribution<RealType, Policy>, RealType>& c) | |
949 | { // Complemented Cumulative Distribution Function | |
950 | const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::cdf(%1%)"; | |
951 | non_central_chi_squared_distribution<RealType, Policy> const& dist = c.dist; | |
952 | RealType x = c.param; | |
953 | RealType k = dist.degrees_of_freedom(); | |
954 | RealType l = dist.non_centrality(); | |
955 | RealType r; | |
956 | if(!detail::check_df( | |
957 | function, | |
958 | k, &r, Policy()) | |
959 | || | |
960 | !detail::check_non_centrality( | |
961 | function, | |
962 | l, | |
963 | &r, | |
964 | Policy()) | |
965 | || | |
966 | !detail::check_positive_x( | |
967 | function, | |
968 | x, | |
969 | &r, | |
970 | Policy())) | |
971 | return r; | |
972 | ||
973 | return detail::non_central_chi_squared_cdf(x, k, l, true, Policy()); | |
974 | } // ccdf | |
975 | ||
976 | template <class RealType, class Policy> | |
977 | inline RealType quantile(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& p) | |
978 | { // Quantile (or Percent Point) function. | |
979 | return detail::nccs_quantile(dist, p, false); | |
980 | } // quantile | |
981 | ||
982 | template <class RealType, class Policy> | |
983 | inline RealType quantile(const complemented2_type<non_central_chi_squared_distribution<RealType, Policy>, RealType>& c) | |
984 | { // Quantile (or Percent Point) function. | |
985 | return detail::nccs_quantile(c.dist, c.param, true); | |
986 | } // quantile complement. | |
987 | ||
988 | } // namespace math | |
989 | } // namespace boost | |
990 | ||
991 | // This include must be at the end, *after* the accessors | |
992 | // for this distribution have been defined, in order to | |
993 | // keep compilers that support two-phase lookup happy. | |
994 | #include <boost/math/distributions/detail/derived_accessors.hpp> | |
995 | ||
996 | #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP | |
997 | ||
998 | ||
999 |