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