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