1 // boost\math\distributions\non_central_chi_squared.hpp 2 3 // Copyright John Maddock 2008. 4 5 // Use, modification and distribution are subject to the 6 // Boost Software License, Version 1.0. 7 // (See accompanying file LICENSE_1_0.txt 8 // or copy at http://www.boost.org/LICENSE_1_0.txt) 9 10 #ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP 11 #define BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP 12 13 #include <boost/math/distributions/fwd.hpp> 14 #include <boost/math/special_functions/gamma.hpp> // for incomplete gamma. gamma_q 15 #include <boost/math/special_functions/bessel.hpp> // for cyl_bessel_i 16 #include <boost/math/special_functions/round.hpp> // for iround 17 #include <boost/math/distributions/complement.hpp> // complements 18 #include <boost/math/distributions/chi_squared.hpp> // central distribution 19 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks 20 #include <boost/math/special_functions/fpclassify.hpp> // isnan. 21 #include <boost/math/tools/roots.hpp> // for root finding. 22 #include <boost/math/distributions/detail/generic_mode.hpp> 23 #include <boost/math/distributions/detail/generic_quantile.hpp> 24 25 namespace boost 26 { 27 namespace math 28 { 29 30 template <class RealType, class Policy> 31 class non_central_chi_squared_distribution; 32 33 namespace detail{ 34 35 template <class T, class Policy> 36 T non_central_chi_square_q(T x, T f, T theta, const Policy& pol, T init_sum = 0) 37 { 38 // 39 // Computes the complement of the Non-Central Chi-Square 40 // Distribution CDF by summing a weighted sum of complements 41 // of the central-distributions. The weighting factor is 42 // a Poisson Distribution. 43 // 44 // This is an application of the technique described in: 45 // 46 // Computing discrete mixtures of continuous 47 // distributions: noncentral chisquare, noncentral t 48 // and the distribution of the square of the sample 49 // multiple correlation coeficient. 50 // D. Benton, K. Krishnamoorthy. 51 // Computational Statistics & Data Analysis 43 (2003) 249 - 267 52 // 53 BOOST_MATH_STD_USING 54 55 // Special case: 56 if(x == 0) 57 return 1; 58 59 // 60 // Initialize the variables we'll be using: 61 // 62 T lambda = theta / 2; 63 T del = f / 2; 64 T y = x / 2; 65 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); 66 T errtol = boost::math::policies::get_epsilon<T, Policy>(); 67 T sum = init_sum; 68 // 69 // k is the starting location for iteration, we'll 70 // move both forwards and backwards from this point. 71 // k is chosen as the peek of the Poisson weights, which 72 // will occur *before* the largest term. 73 // 74 int k = iround(lambda, pol); 75 // Forwards and backwards Poisson weights: 76 T poisf = boost::math::gamma_p_derivative(1 + k, lambda, pol); 77 T poisb = poisf * k / lambda; 78 // Initial forwards central chi squared term: 79 T gamf = boost::math::gamma_q(del + k, y, pol); 80 // Forwards and backwards recursion terms on the central chi squared: 81 T xtermf = boost::math::gamma_p_derivative(del + 1 + k, y, pol); 82 T xtermb = xtermf * (del + k) / y; 83 // Initial backwards central chi squared term: 84 T gamb = gamf - xtermb; 85 86 // 87 // Forwards iteration first, this is the 88 // stable direction for the gamma function 89 // recurrences: 90 // 91 int i; 92 for(i = k; static_cast<boost::uintmax_t>(i-k) < max_iter; ++i) 93 { 94 T term = poisf * gamf; 95 sum += term; 96 poisf *= lambda / (i + 1); 97 gamf += xtermf; 98 xtermf *= y / (del + i + 1); 99 if(((sum == 0) || (fabs(term / sum) < errtol)) && (term >= poisf * gamf)) 100 break; 101 } 102 //Error check: 103 if(static_cast<boost::uintmax_t>(i-k) >= max_iter) 104 return policies::raise_evaluation_error( 105 "cdf(non_central_chi_squared_distribution<%1%>, %1%)", 106 "Series did not converge, closest value was %1%", sum, pol); 107 // 108 // Now backwards iteration: the gamma 109 // function recurrences are unstable in this 110 // direction, we rely on the terms deminishing in size 111 // faster than we introduce cancellation errors. 112 // For this reason it's very important that we start 113 // *before* the largest term so that backwards iteration 114 // is strictly converging. 115 // 116 for(i = k - 1; i >= 0; --i) 117 { 118 T term = poisb * gamb; 119 sum += term; 120 poisb *= i / lambda; 121 xtermb *= (del + i) / y; 122 gamb -= xtermb; 123 if((sum == 0) || (fabs(term / sum) < errtol)) 124 break; 125 } 126 127 return sum; 128 } 129 130 template <class T, class Policy> 131 T non_central_chi_square_p_ding(T x, T f, T theta, const Policy& pol, T init_sum = 0) 132 { 133 // 134 // This is an implementation of: 135 // 136 // Algorithm AS 275: 137 // Computing the Non-Central #2 Distribution Function 138 // Cherng G. Ding 139 // Applied Statistics, Vol. 41, No. 2. (1992), pp. 478-482. 140 // 141 // This uses a stable forward iteration to sum the 142 // CDF, unfortunately this can not be used for large 143 // values of the non-centrality parameter because: 144 // * The first term may underfow to zero. 145 // * We may need an extra-ordinary number of terms 146 // before we reach the first *significant* term. 147 // 148 BOOST_MATH_STD_USING 149 // Special case: 150 if(x == 0) 151 return 0; 152 T tk = boost::math::gamma_p_derivative(f/2 + 1, x/2, pol); 153 T lambda = theta / 2; 154 T vk = exp(-lambda); 155 T uk = vk; 156 T sum = init_sum + tk * vk; 157 if(sum == 0) 158 return sum; 159 160 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); 161 T errtol = boost::math::policies::get_epsilon<T, Policy>(); 162 163 int i; 164 T lterm(0), term(0); 165 for(i = 1; static_cast<boost::uintmax_t>(i) < max_iter; ++i) 166 { 167 tk = tk * x / (f + 2 * i); 168 uk = uk * lambda / i; 169 vk = vk + uk; 170 lterm = term; 171 term = vk * tk; 172 sum += term; 173 if((fabs(term / sum) < errtol) && (term <= lterm)) 174 break; 175 } 176 //Error check: 177 if(static_cast<boost::uintmax_t>(i) >= max_iter) 178 return policies::raise_evaluation_error( 179 "cdf(non_central_chi_squared_distribution<%1%>, %1%)", 180 "Series did not converge, closest value was %1%", sum, pol); 181 return sum; 182 } 183 184 185 template <class T, class Policy> 186 T non_central_chi_square_p(T y, T n, T lambda, const Policy& pol, T init_sum) 187 { 188 // 189 // This is taken more or less directly from: 190 // 191 // Computing discrete mixtures of continuous 192 // distributions: noncentral chisquare, noncentral t 193 // and the distribution of the square of the sample 194 // multiple correlation coeficient. 195 // D. Benton, K. Krishnamoorthy. 196 // Computational Statistics & Data Analysis 43 (2003) 249 - 267 197 // 198 // We're summing a Poisson weighting term multiplied by 199 // a central chi squared distribution. 200 // 201 BOOST_MATH_STD_USING 202 // Special case: 203 if(y == 0) 204 return 0; 205 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); 206 T errtol = boost::math::policies::get_epsilon<T, Policy>(); 207 T errorf(0), errorb(0); 208 209 T x = y / 2; 210 T del = lambda / 2; 211 // 212 // Starting location for the iteration, we'll iterate 213 // both forwards and backwards from this point. The 214 // location chosen is the maximum of the Poisson weight 215 // function, which ocurrs *after* the largest term in the 216 // sum. 217 // 218 int k = iround(del, pol); 219 T a = n / 2 + k; 220 // Central chi squared term for forward iteration: 221 T gamkf = boost::math::gamma_p(a, x, pol); 222 223 if(lambda == 0) 224 return gamkf; 225 // Central chi squared term for backward iteration: 226 T gamkb = gamkf; 227 // Forwards Poisson weight: 228 T poiskf = gamma_p_derivative(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(k + 1, l2, pol) * gamma_p_derivative(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> non_central_chi_squared_cdf(RealType x,RealType k,RealType l,bool invert,const 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 result = cdf(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 { nccs_quantile_functorboost::math::detail::nccs_quantile_functor384 nccs_quantile_functor(const non_central_chi_squared_distribution<T,Policy>& d, T t, bool c) 385 : dist(d), target(t), comp(c) {} 386 operator ()boost::math::detail::nccs_quantile_functor387 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 // 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 RealType nccs_pdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x) 490 { 491 BOOST_MATH_STD_USING 492 static const char* 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 (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 else 532 { 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 else 539 { 540 r = exp(r); 541 r = 0.5f * r 542 * 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_finder 552 { degrees_of_freedom_finderboost::math::detail::degrees_of_freedom_finder553 degrees_of_freedom_finder( 554 RealType lam_, RealType x_, RealType p_, bool c) 555 : lam(lam_), x(x_), p(p_), comp(c) {} 556 operator ()boost::math::detail::degrees_of_freedom_finder557 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> find_degrees_of_freedom(RealType lam,RealType x,RealType p,RealType q,const Policy & pol)572 inline RealType find_degrees_of_freedom( 573 RealType lam, RealType x, RealType p, RealType q, const Policy& pol) 574 { 575 const char* 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, 582 "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%", 583 RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); 584 } 585 degrees_of_freedom_finder<RealType, Policy> f(lam, x, p < q ? p : q, p < q ? false : true); 586 tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>()); 587 boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); 588 // 589 // Pick an initial guess that we know will give us a probability 590 // right around 0.5. 591 // 592 RealType guess = x - lam; 593 if(guess < 1) 594 guess = 1; 595 std::pair<RealType, RealType> ir = tools::bracket_and_solve_root( 596 f, guess, RealType(2), false, tol, max_iter, pol); 597 RealType result = ir.first + (ir.second - ir.first) / 2; 598 if(max_iter >= policies::get_max_root_iterations<Policy>()) 599 { 600 return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" 601 " or there is no answer to problem. Current best guess is %1%", result, Policy()); 602 } 603 return result; 604 } 605 606 template <class RealType, class Policy> 607 struct non_centrality_finder 608 { non_centrality_finderboost::math::detail::non_centrality_finder609 non_centrality_finder( 610 RealType v_, RealType x_, RealType p_, bool c) 611 : v(v_), x(x_), p(p_), comp(c) {} 612 operator ()boost::math::detail::non_centrality_finder613 RealType operator()(const RealType& lam) 614 { 615 non_central_chi_squared_distribution<RealType, Policy> d(v, lam); 616 return comp ? 617 RealType(p - cdf(complement(d, x))) 618 : RealType(cdf(d, x) - p); 619 } 620 private: 621 RealType v; 622 RealType x; 623 RealType p; 624 bool comp; 625 }; 626 627 template <class RealType, class Policy> find_non_centrality(RealType v,RealType x,RealType p,RealType q,const Policy & pol)628 inline RealType find_non_centrality( 629 RealType v, RealType x, RealType p, RealType q, const Policy& pol) 630 { 631 const char* function = "non_central_chi_squared<%1%>::find_non_centrality"; 632 if((p == 0) || (q == 0)) 633 { 634 // 635 // Can't do a thing if one of p and q is zero: 636 // 637 return policies::raise_evaluation_error<RealType>(function, 638 "Can't find non centrality parameter when the probability is 0 or 1, only possible answer is %1%", 639 RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); 640 } 641 non_centrality_finder<RealType, Policy> f(v, x, p < q ? p : q, p < q ? false : true); 642 tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>()); 643 boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); 644 // 645 // Pick an initial guess that we know will give us a probability 646 // right around 0.5. 647 // 648 RealType guess = x - v; 649 if(guess < 1) 650 guess = 1; 651 std::pair<RealType, RealType> ir = tools::bracket_and_solve_root( 652 f, guess, RealType(2), false, tol, max_iter, pol); 653 RealType result = ir.first + (ir.second - ir.first) / 2; 654 if(max_iter >= policies::get_max_root_iterations<Policy>()) 655 { 656 return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" 657 " or there is no answer to problem. Current best guess is %1%", result, Policy()); 658 } 659 return result; 660 } 661 662 } 663 664 template <class RealType = double, class Policy = policies::policy<> > 665 class non_central_chi_squared_distribution 666 { 667 public: 668 typedef RealType value_type; 669 typedef Policy policy_type; 670 non_central_chi_squared_distribution(RealType df_,RealType lambda)671 non_central_chi_squared_distribution(RealType df_, RealType lambda) : df(df_), ncp(lambda) 672 { 673 const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::non_central_chi_squared_distribution(%1%,%1%)"; 674 RealType r; 675 detail::check_df( 676 function, 677 df, &r, Policy()); 678 detail::check_non_centrality( 679 function, 680 ncp, 681 &r, 682 Policy()); 683 } // non_central_chi_squared_distribution constructor. 684 degrees_of_freedom() const685 RealType degrees_of_freedom() const 686 { // Private data getter function. 687 return df; 688 } non_centrality() const689 RealType non_centrality() const 690 { // Private data getter function. 691 return ncp; 692 } find_degrees_of_freedom(RealType lam,RealType x,RealType p)693 static RealType find_degrees_of_freedom(RealType lam, RealType x, RealType p) 694 { 695 const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom"; 696 typedef typename policies::evaluation<RealType, Policy>::type value_type; 697 typedef typename policies::normalise< 698 Policy, 699 policies::promote_float<false>, 700 policies::promote_double<false>, 701 policies::discrete_quantile<>, 702 policies::assert_undefined<> >::type forwarding_policy; 703 value_type result = detail::find_degrees_of_freedom( 704 static_cast<value_type>(lam), 705 static_cast<value_type>(x), 706 static_cast<value_type>(p), 707 static_cast<value_type>(1-p), 708 forwarding_policy()); 709 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 710 result, 711 function); 712 } 713 template <class A, class B, class C> find_degrees_of_freedom(const complemented3_type<A,B,C> & c)714 static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c) 715 { 716 const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom"; 717 typedef typename policies::evaluation<RealType, Policy>::type value_type; 718 typedef typename policies::normalise< 719 Policy, 720 policies::promote_float<false>, 721 policies::promote_double<false>, 722 policies::discrete_quantile<>, 723 policies::assert_undefined<> >::type forwarding_policy; 724 value_type result = detail::find_degrees_of_freedom( 725 static_cast<value_type>(c.dist), 726 static_cast<value_type>(c.param1), 727 static_cast<value_type>(1-c.param2), 728 static_cast<value_type>(c.param2), 729 forwarding_policy()); 730 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 731 result, 732 function); 733 } find_non_centrality(RealType v,RealType x,RealType p)734 static RealType find_non_centrality(RealType v, RealType x, RealType p) 735 { 736 const char* function = "non_central_chi_squared<%1%>::find_non_centrality"; 737 typedef typename policies::evaluation<RealType, Policy>::type value_type; 738 typedef typename policies::normalise< 739 Policy, 740 policies::promote_float<false>, 741 policies::promote_double<false>, 742 policies::discrete_quantile<>, 743 policies::assert_undefined<> >::type forwarding_policy; 744 value_type result = detail::find_non_centrality( 745 static_cast<value_type>(v), 746 static_cast<value_type>(x), 747 static_cast<value_type>(p), 748 static_cast<value_type>(1-p), 749 forwarding_policy()); 750 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 751 result, 752 function); 753 } 754 template <class A, class B, class C> find_non_centrality(const complemented3_type<A,B,C> & c)755 static RealType find_non_centrality(const complemented3_type<A,B,C>& c) 756 { 757 const char* function = "non_central_chi_squared<%1%>::find_non_centrality"; 758 typedef typename policies::evaluation<RealType, Policy>::type value_type; 759 typedef typename policies::normalise< 760 Policy, 761 policies::promote_float<false>, 762 policies::promote_double<false>, 763 policies::discrete_quantile<>, 764 policies::assert_undefined<> >::type forwarding_policy; 765 value_type result = detail::find_non_centrality( 766 static_cast<value_type>(c.dist), 767 static_cast<value_type>(c.param1), 768 static_cast<value_type>(1-c.param2), 769 static_cast<value_type>(c.param2), 770 forwarding_policy()); 771 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 772 result, 773 function); 774 } 775 private: 776 // Data member, initialized by constructor. 777 RealType df; // degrees of freedom. 778 RealType ncp; // non-centrality parameter 779 }; // template <class RealType, class Policy> class non_central_chi_squared_distribution 780 781 typedef non_central_chi_squared_distribution<double> non_central_chi_squared; // Reserved name of type double. 782 783 // Non-member functions to give properties of the distribution. 784 785 template <class RealType, class Policy> range(const non_central_chi_squared_distribution<RealType,Policy> &)786 inline const std::pair<RealType, RealType> range(const non_central_chi_squared_distribution<RealType, Policy>& /* dist */) 787 { // Range of permissible values for random variable k. 788 using boost::math::tools::max_value; 789 return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // Max integer? 790 } 791 792 template <class RealType, class Policy> support(const non_central_chi_squared_distribution<RealType,Policy> &)793 inline const std::pair<RealType, RealType> support(const non_central_chi_squared_distribution<RealType, Policy>& /* dist */) 794 { // Range of supported values for random variable k. 795 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. 796 using boost::math::tools::max_value; 797 return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); 798 } 799 800 template <class RealType, class Policy> mean(const non_central_chi_squared_distribution<RealType,Policy> & dist)801 inline RealType mean(const non_central_chi_squared_distribution<RealType, Policy>& dist) 802 { // Mean of poisson distribution = lambda. 803 const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::mean()"; 804 RealType k = dist.degrees_of_freedom(); 805 RealType l = dist.non_centrality(); 806 RealType r; 807 if(!detail::check_df( 808 function, 809 k, &r, Policy()) 810 || 811 !detail::check_non_centrality( 812 function, 813 l, 814 &r, 815 Policy())) 816 return r; 817 return k + l; 818 } // mean 819 820 template <class RealType, class Policy> mode(const non_central_chi_squared_distribution<RealType,Policy> & dist)821 inline RealType mode(const non_central_chi_squared_distribution<RealType, Policy>& dist) 822 { // mode. 823 static const char* function = "mode(non_central_chi_squared_distribution<%1%> const&)"; 824 825 RealType k = dist.degrees_of_freedom(); 826 RealType l = dist.non_centrality(); 827 RealType r; 828 if(!detail::check_df( 829 function, 830 k, &r, Policy()) 831 || 832 !detail::check_non_centrality( 833 function, 834 l, 835 &r, 836 Policy())) 837 return (RealType)r; 838 return detail::generic_find_mode(dist, 1 + k, function); 839 } 840 841 template <class RealType, class Policy> variance(const non_central_chi_squared_distribution<RealType,Policy> & dist)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 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> skewness(const non_central_chi_squared_distribution<RealType,Policy> & dist)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 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> kurtosis_excess(const non_central_chi_squared_distribution<RealType,Policy> & dist)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 r; 902 return 12 * (k + 4 * l) / ((k + 2 * l) * (k + 2 * l)); 903 } // kurtosis_excess 904 905 template <class RealType, class Policy> kurtosis(const non_central_chi_squared_distribution<RealType,Policy> & dist)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> pdf(const non_central_chi_squared_distribution<RealType,Policy> & dist,const RealType & x)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 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 r; 969 970 return detail::non_central_chi_squared_cdf(x, k, l, true, Policy()); 971 } // ccdf 972 973 template <class RealType, class Policy> quantile(const non_central_chi_squared_distribution<RealType,Policy> & dist,const RealType & p)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> quantile(const complemented2_type<non_central_chi_squared_distribution<RealType,Policy>,RealType> & c)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 995 996 997