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

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