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
25namespace 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 coefficient.
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 diminishing 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 underflow 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 coefficient.
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 } // pdf
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
1000

source code of include/boost/math/distributions/non_central_chi_squared.hpp