1 // boost\math\distributions\non_central_beta.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_BETA_HPP 11 #define BOOST_MATH_SPECIAL_NON_CENTRAL_BETA_HPP 12 13 #include <boost/math/distributions/fwd.hpp> 14 #include <boost/math/special_functions/beta.hpp> // for incomplete gamma. gamma_q 15 #include <boost/math/distributions/complement.hpp> // complements 16 #include <boost/math/distributions/beta.hpp> // central distribution 17 #include <boost/math/distributions/detail/generic_mode.hpp> 18 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks 19 #include <boost/math/special_functions/fpclassify.hpp> // isnan. 20 #include <boost/math/tools/roots.hpp> // for root finding. 21 #include <boost/math/tools/series.hpp> 22 23 namespace boost 24 { 25 namespace math 26 { 27 28 template <class RealType, class Policy> 29 class non_central_beta_distribution; 30 31 namespace detail{ 32 33 template <class T, class Policy> 34 T non_central_beta_p(T a, T b, T lam, T x, T y, const Policy& pol, T init_val = 0) 35 { 36 BOOST_MATH_STD_USING 37 using namespace boost::math; 38 // 39 // Variables come first: 40 // 41 std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); 42 T errtol = boost::math::policies::get_epsilon<T, Policy>(); 43 T l2 = lam / 2; 44 // 45 // k is the starting point for iteration, and is the 46 // maximum of the poisson weighting term, 47 // note that unlike other similar code, we do not set 48 // k to zero, when l2 is small, as forward iteration 49 // is unstable: 50 // 51 int k = itrunc(l2); 52 if(k == 0) 53 k = 1; 54 // Starting Poisson weight: 55 T pois = gamma_p_derivative(T(k+1), l2, pol); 56 if(pois == 0) 57 return init_val; 58 // recurance term: 59 T xterm; 60 // Starting beta term: 61 T beta = x < y 62 ? detail::ibeta_imp(T(a + k), b, x, pol, false, true, &xterm) 63 : detail::ibeta_imp(b, T(a + k), y, pol, true, true, &xterm); 64 65 xterm *= y / (a + b + k - 1); 66 T poisf(pois), betaf(beta), xtermf(xterm); 67 T sum = init_val; 68 69 if((beta == 0) && (xterm == 0)) 70 return init_val; 71 72 // 73 // Backwards recursion first, this is the stable 74 // direction for recursion: 75 // 76 T last_term = 0; 77 std::uintmax_t count = k; 78 for(int i = k; i >= 0; --i) 79 { 80 T term = beta * pois; 81 sum += term; 82 if(((fabs(term/sum) < errtol) && (last_term >= term)) || (term == 0)) 83 { 84 count = k - i; 85 break; 86 } 87 pois *= i / l2; 88 beta += xterm; 89 xterm *= (a + i - 1) / (x * (a + b + i - 2)); 90 last_term = term; 91 } 92 for(int i = k + 1; ; ++i) 93 { 94 poisf *= l2 / i; 95 xtermf *= (x * (a + b + i - 2)) / (a + i - 1); 96 betaf -= xtermf; 97 98 T term = poisf * betaf; 99 sum += term; 100 if((fabs(term/sum) < errtol) || (term == 0)) 101 { 102 break; 103 } 104 if(static_cast<std::uintmax_t>(count + i - k) > max_iter) 105 { 106 return policies::raise_evaluation_error( 107 "cdf(non_central_beta_distribution<%1%>, %1%)", 108 "Series did not converge, closest value was %1%", sum, pol); 109 } 110 } 111 return sum; 112 } 113 114 template <class T, class Policy> 115 T non_central_beta_q(T a, T b, T lam, T x, T y, const Policy& pol, T init_val = 0) 116 { 117 BOOST_MATH_STD_USING 118 using namespace boost::math; 119 // 120 // Variables come first: 121 // 122 std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); 123 T errtol = boost::math::policies::get_epsilon<T, Policy>(); 124 T l2 = lam / 2; 125 // 126 // k is the starting point for iteration, and is the 127 // maximum of the poisson weighting term: 128 // 129 int k = itrunc(l2); 130 T pois; 131 if(k <= 30) 132 { 133 // 134 // Might as well start at 0 since we'll likely have this number of terms anyway: 135 // 136 if(a + b > 1) 137 k = 0; 138 else if(k == 0) 139 k = 1; 140 } 141 if(k == 0) 142 { 143 // Starting Poisson weight: 144 pois = exp(-l2); 145 } 146 else 147 { 148 // Starting Poisson weight: 149 pois = gamma_p_derivative(T(k+1), l2, pol); 150 } 151 if(pois == 0) 152 return init_val; 153 // recurance term: 154 T xterm; 155 // Starting beta term: 156 T beta = x < y 157 ? detail::ibeta_imp(T(a + k), b, x, pol, true, true, &xterm) 158 : detail::ibeta_imp(b, T(a + k), y, pol, false, true, &xterm); 159 160 xterm *= y / (a + b + k - 1); 161 T poisf(pois), betaf(beta), xtermf(xterm); 162 T sum = init_val; 163 if((beta == 0) && (xterm == 0)) 164 return init_val; 165 // 166 // Forwards recursion first, this is the stable 167 // direction for recursion, and the location 168 // of the bulk of the sum: 169 // 170 T last_term = 0; 171 std::uintmax_t count = 0; 172 for(int i = k + 1; ; ++i) 173 { 174 poisf *= l2 / i; 175 xtermf *= (x * (a + b + i - 2)) / (a + i - 1); 176 betaf += xtermf; 177 178 T term = poisf * betaf; 179 sum += term; 180 if((fabs(term/sum) < errtol) && (last_term >= term)) 181 { 182 count = i - k; 183 break; 184 } 185 if(static_cast<std::uintmax_t>(i - k) > max_iter) 186 { 187 return policies::raise_evaluation_error( 188 "cdf(non_central_beta_distribution<%1%>, %1%)", 189 "Series did not converge, closest value was %1%", sum, pol); 190 } 191 last_term = term; 192 } 193 for(int i = k; i >= 0; --i) 194 { 195 T term = beta * pois; 196 sum += term; 197 if(fabs(term/sum) < errtol) 198 { 199 break; 200 } 201 if(static_cast<std::uintmax_t>(count + k - i) > max_iter) 202 { 203 return policies::raise_evaluation_error( 204 "cdf(non_central_beta_distribution<%1%>, %1%)", 205 "Series did not converge, closest value was %1%", sum, pol); 206 } 207 pois *= i / l2; 208 beta -= xterm; 209 xterm *= (a + i - 1) / (x * (a + b + i - 2)); 210 } 211 return sum; 212 } 213 214 template <class RealType, class Policy> non_central_beta_cdf(RealType x,RealType y,RealType a,RealType b,RealType l,bool invert,const Policy &)215 inline RealType non_central_beta_cdf(RealType x, RealType y, RealType a, RealType b, RealType l, bool invert, const Policy&) 216 { 217 typedef typename policies::evaluation<RealType, Policy>::type value_type; 218 typedef typename policies::normalise< 219 Policy, 220 policies::promote_float<false>, 221 policies::promote_double<false>, 222 policies::discrete_quantile<>, 223 policies::assert_undefined<> >::type forwarding_policy; 224 225 BOOST_MATH_STD_USING 226 227 if(x == 0) 228 return invert ? 1.0f : 0.0f; 229 if(y == 0) 230 return invert ? 0.0f : 1.0f; 231 value_type result; 232 value_type c = a + b + l / 2; 233 value_type cross = 1 - (b / c) * (1 + l / (2 * c * c)); 234 if(l == 0) 235 result = cdf(boost::math::beta_distribution<RealType, Policy>(a, b), x); 236 else if(x > cross) 237 { 238 // Complement is the smaller of the two: 239 result = detail::non_central_beta_q( 240 static_cast<value_type>(a), 241 static_cast<value_type>(b), 242 static_cast<value_type>(l), 243 static_cast<value_type>(x), 244 static_cast<value_type>(y), 245 forwarding_policy(), 246 static_cast<value_type>(invert ? 0 : -1)); 247 invert = !invert; 248 } 249 else 250 { 251 result = detail::non_central_beta_p( 252 static_cast<value_type>(a), 253 static_cast<value_type>(b), 254 static_cast<value_type>(l), 255 static_cast<value_type>(x), 256 static_cast<value_type>(y), 257 forwarding_policy(), 258 static_cast<value_type>(invert ? -1 : 0)); 259 } 260 if(invert) 261 result = -result; 262 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 263 result, 264 "boost::math::non_central_beta_cdf<%1%>(%1%, %1%, %1%)"); 265 } 266 267 template <class T, class Policy> 268 struct nc_beta_quantile_functor 269 { nc_beta_quantile_functorboost::math::detail::nc_beta_quantile_functor270 nc_beta_quantile_functor(const non_central_beta_distribution<T,Policy>& d, T t, bool c) 271 : dist(d), target(t), comp(c) {} 272 operator ()boost::math::detail::nc_beta_quantile_functor273 T operator()(const T& x) 274 { 275 return comp ? 276 T(target - cdf(complement(dist, x))) 277 : T(cdf(dist, x) - target); 278 } 279 280 private: 281 non_central_beta_distribution<T,Policy> dist; 282 T target; 283 bool comp; 284 }; 285 286 // 287 // This is more or less a copy of bracket_and_solve_root, but 288 // modified to search only the interval [0,1] using similar 289 // heuristics. 290 // 291 template <class F, class T, class Tol, class Policy> bracket_and_solve_root_01(F f,const T & guess,T factor,bool rising,Tol tol,std::uintmax_t & max_iter,const Policy & pol)292 std::pair<T, T> bracket_and_solve_root_01(F f, const T& guess, T factor, bool rising, Tol tol, std::uintmax_t& max_iter, const Policy& pol) 293 { 294 BOOST_MATH_STD_USING 295 static const char* function = "boost::math::tools::bracket_and_solve_root_01<%1%>"; 296 // 297 // Set up initial brackets: 298 // 299 T a = guess; 300 T b = a; 301 T fa = f(a); 302 T fb = fa; 303 // 304 // Set up invocation count: 305 // 306 std::uintmax_t count = max_iter - 1; 307 308 if((fa < 0) == (guess < 0 ? !rising : rising)) 309 { 310 // 311 // Zero is to the right of b, so walk upwards 312 // until we find it: 313 // 314 while((boost::math::sign)(fb) == (boost::math::sign)(fa)) 315 { 316 if(count == 0) 317 { 318 b = policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, pol); 319 return std::make_pair(a, b); 320 } 321 // 322 // Heuristic: every 20 iterations we double the growth factor in case the 323 // initial guess was *really* bad ! 324 // 325 if((max_iter - count) % 20 == 0) 326 factor *= 2; 327 // 328 // Now go ahead and move are guess by "factor", 329 // we do this by reducing 1-guess by factor: 330 // 331 a = b; 332 fa = fb; 333 b = 1 - ((1 - b) / factor); 334 fb = f(b); 335 --count; 336 BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count); 337 } 338 } 339 else 340 { 341 // 342 // Zero is to the left of a, so walk downwards 343 // until we find it: 344 // 345 while((boost::math::sign)(fb) == (boost::math::sign)(fa)) 346 { 347 if(fabs(a) < tools::min_value<T>()) 348 { 349 // Escape route just in case the answer is zero! 350 max_iter -= count; 351 max_iter += 1; 352 return a > 0 ? std::make_pair(T(0), T(a)) : std::make_pair(T(a), T(0)); 353 } 354 if(count == 0) 355 { 356 a = policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, pol); 357 return std::make_pair(a, b); 358 } 359 // 360 // Heuristic: every 20 iterations we double the growth factor in case the 361 // initial guess was *really* bad ! 362 // 363 if((max_iter - count) % 20 == 0) 364 factor *= 2; 365 // 366 // Now go ahead and move are guess by "factor": 367 // 368 b = a; 369 fb = fa; 370 a /= factor; 371 fa = f(a); 372 --count; 373 BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count); 374 } 375 } 376 max_iter -= count; 377 max_iter += 1; 378 std::pair<T, T> r = toms748_solve( 379 f, 380 (a < 0 ? b : a), 381 (a < 0 ? a : b), 382 (a < 0 ? fb : fa), 383 (a < 0 ? fa : fb), 384 tol, 385 count, 386 pol); 387 max_iter += count; 388 BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count); 389 return r; 390 } 391 392 template <class RealType, class Policy> 393 RealType nc_beta_quantile(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& p, bool comp) 394 { 395 static const char* function = "quantile(non_central_beta_distribution<%1%>, %1%)"; 396 typedef typename policies::evaluation<RealType, Policy>::type value_type; 397 typedef typename policies::normalise< 398 Policy, 399 policies::promote_float<false>, 400 policies::promote_double<false>, 401 policies::discrete_quantile<>, 402 policies::assert_undefined<> >::type forwarding_policy; 403 404 value_type a = dist.alpha(); 405 value_type b = dist.beta(); 406 value_type l = dist.non_centrality(); 407 value_type r; 408 if(!beta_detail::check_alpha( 409 function, 410 a, &r, Policy()) 411 || 412 !beta_detail::check_beta( 413 function, 414 b, &r, Policy()) 415 || 416 !detail::check_non_centrality( 417 function, 418 l, 419 &r, 420 Policy()) 421 || 422 !detail::check_probability( 423 function, 424 static_cast<value_type>(p), 425 &r, 426 Policy())) 427 return (RealType)r; 428 // 429 // Special cases first: 430 // 431 if(p == 0) 432 return comp 433 ? 1.0f 434 : 0.0f; 435 if(p == 1) 436 return !comp 437 ? 1.0f 438 : 0.0f; 439 440 value_type c = a + b + l / 2; 441 value_type mean = 1 - (b / c) * (1 + l / (2 * c * c)); 442 /* 443 // 444 // Calculate a normal approximation to the quantile, 445 // uses mean and variance approximations from: 446 // Algorithm AS 310: 447 // Computing the Non-Central Beta Distribution Function 448 // R. Chattamvelli; R. Shanmugam 449 // Applied Statistics, Vol. 46, No. 1. (1997), pp. 146-156. 450 // 451 // Unfortunately, when this is wrong it tends to be *very* 452 // wrong, so it's disabled for now, even though it often 453 // gets the initial guess quite close. Probably we could 454 // do much better by factoring in the skewness if only 455 // we could calculate it.... 456 // 457 value_type delta = l / 2; 458 value_type delta2 = delta * delta; 459 value_type delta3 = delta * delta2; 460 value_type delta4 = delta2 * delta2; 461 value_type G = c * (c + 1) + delta; 462 value_type alpha = a + b; 463 value_type alpha2 = alpha * alpha; 464 value_type eta = (2 * alpha + 1) * (2 * alpha + 1) + 1; 465 value_type H = 3 * alpha2 + 5 * alpha + 2; 466 value_type F = alpha2 * (alpha + 1) + H * delta 467 + (2 * alpha + 4) * delta2 + delta3; 468 value_type P = (3 * alpha + 1) * (9 * alpha + 17) 469 + 2 * alpha * (3 * alpha + 2) * (3 * alpha + 4) + 15; 470 value_type Q = 54 * alpha2 + 162 * alpha + 130; 471 value_type R = 6 * (6 * alpha + 11); 472 value_type D = delta 473 * (H * H + 2 * P * delta + Q * delta2 + R * delta3 + 9 * delta4); 474 value_type variance = (b / G) 475 * (1 + delta * (l * l + 3 * l + eta) / (G * G)) 476 - (b * b / F) * (1 + D / (F * F)); 477 value_type sd = sqrt(variance); 478 479 value_type guess = comp 480 ? quantile(complement(normal_distribution<RealType, Policy>(static_cast<RealType>(mean), static_cast<RealType>(sd)), p)) 481 : quantile(normal_distribution<RealType, Policy>(static_cast<RealType>(mean), static_cast<RealType>(sd)), p); 482 483 if(guess >= 1) 484 guess = mean; 485 if(guess <= tools::min_value<value_type>()) 486 guess = mean; 487 */ 488 value_type guess = mean; 489 detail::nc_beta_quantile_functor<value_type, Policy> 490 f(non_central_beta_distribution<value_type, Policy>(a, b, l), p, comp); 491 tools::eps_tolerance<value_type> tol(policies::digits<RealType, Policy>()); 492 std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); 493 494 std::pair<value_type, value_type> ir 495 = bracket_and_solve_root_01( 496 f, guess, value_type(2.5), true, tol, 497 max_iter, Policy()); 498 value_type result = ir.first + (ir.second - ir.first) / 2; 499 500 if(max_iter >= policies::get_max_root_iterations<Policy>()) 501 { 502 return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" 503 " either there is no answer to quantile of the non central beta distribution" 504 " or the answer is infinite. Current best guess is %1%", 505 policies::checked_narrowing_cast<RealType, forwarding_policy>( 506 result, 507 function), Policy()); 508 } 509 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 510 result, 511 function); 512 } 513 514 template <class T, class Policy> 515 T non_central_beta_pdf(T a, T b, T lam, T x, T y, const Policy& pol) 516 { 517 BOOST_MATH_STD_USING 518 // 519 // Special cases: 520 // 521 if((x == 0) || (y == 0)) 522 return 0; 523 // 524 // Variables come first: 525 // 526 std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); 527 T errtol = boost::math::policies::get_epsilon<T, Policy>(); 528 T l2 = lam / 2; 529 // 530 // k is the starting point for iteration, and is the 531 // maximum of the poisson weighting term: 532 // 533 int k = itrunc(l2); 534 // Starting Poisson weight: 535 T pois = gamma_p_derivative(T(k+1), l2, pol); 536 // Starting beta term: 537 T beta = x < y ? 538 ibeta_derivative(a + k, b, x, pol) 539 : ibeta_derivative(b, a + k, y, pol); 540 T sum = 0; 541 T poisf(pois); 542 T betaf(beta); 543 544 // 545 // Stable backwards recursion first: 546 // 547 std::uintmax_t count = k; 548 for(int i = k; i >= 0; --i) 549 { 550 T term = beta * pois; 551 sum += term; 552 if((fabs(term/sum) < errtol) || (term == 0)) 553 { 554 count = k - i; 555 break; 556 } 557 pois *= i / l2; 558 beta *= (a + i - 1) / (x * (a + i + b - 1)); 559 } 560 for(int i = k + 1; ; ++i) 561 { 562 poisf *= l2 / i; 563 betaf *= x * (a + b + i - 1) / (a + i - 1); 564 565 T term = poisf * betaf; 566 sum += term; 567 if((fabs(term/sum) < errtol) || (term == 0)) 568 { 569 break; 570 } 571 if(static_cast<std::uintmax_t>(count + i - k) > max_iter) 572 { 573 return policies::raise_evaluation_error( 574 "pdf(non_central_beta_distribution<%1%>, %1%)", 575 "Series did not converge, closest value was %1%", sum, pol); 576 } 577 } 578 return sum; 579 } 580 581 template <class RealType, class Policy> 582 RealType nc_beta_pdf(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& x) 583 { 584 BOOST_MATH_STD_USING 585 static const char* function = "pdf(non_central_beta_distribution<%1%>, %1%)"; 586 typedef typename policies::evaluation<RealType, Policy>::type value_type; 587 typedef typename policies::normalise< 588 Policy, 589 policies::promote_float<false>, 590 policies::promote_double<false>, 591 policies::discrete_quantile<>, 592 policies::assert_undefined<> >::type forwarding_policy; 593 594 value_type a = dist.alpha(); 595 value_type b = dist.beta(); 596 value_type l = dist.non_centrality(); 597 value_type r; 598 if(!beta_detail::check_alpha( 599 function, 600 a, &r, Policy()) 601 || 602 !beta_detail::check_beta( 603 function, 604 b, &r, Policy()) 605 || 606 !detail::check_non_centrality( 607 function, 608 l, 609 &r, 610 Policy()) 611 || 612 !beta_detail::check_x( 613 function, 614 static_cast<value_type>(x), 615 &r, 616 Policy())) 617 return (RealType)r; 618 619 if(l == 0) 620 return pdf(boost::math::beta_distribution<RealType, Policy>(dist.alpha(), dist.beta()), x); 621 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 622 non_central_beta_pdf(a, b, l, static_cast<value_type>(x), value_type(1 - static_cast<value_type>(x)), forwarding_policy()), 623 "function"); 624 } 625 626 template <class T> 627 struct hypergeometric_2F2_sum 628 { 629 typedef T result_type; hypergeometric_2F2_sumboost::math::detail::hypergeometric_2F2_sum630 hypergeometric_2F2_sum(T a1_, T a2_, T b1_, T b2_, T z_) : a1(a1_), a2(a2_), b1(b1_), b2(b2_), z(z_), term(1), k(0) {} operator ()boost::math::detail::hypergeometric_2F2_sum631 T operator()() 632 { 633 T result = term; 634 term *= a1 * a2 / (b1 * b2); 635 a1 += 1; 636 a2 += 1; 637 b1 += 1; 638 b2 += 1; 639 k += 1; 640 term /= k; 641 term *= z; 642 return result; 643 } 644 T a1, a2, b1, b2, z, term, k; 645 }; 646 647 template <class T, class Policy> 648 T hypergeometric_2F2(T a1, T a2, T b1, T b2, T z, const Policy& pol) 649 { 650 typedef typename policies::evaluation<T, Policy>::type value_type; 651 652 const char* function = "boost::math::detail::hypergeometric_2F2<%1%>(%1%,%1%,%1%,%1%,%1%)"; 653 654 hypergeometric_2F2_sum<value_type> s(a1, a2, b1, b2, z); 655 std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); 656 657 value_type result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<value_type, Policy>(), max_iter); 658 659 policies::check_series_iterations<T>(function, max_iter, pol); 660 return policies::checked_narrowing_cast<T, Policy>(result, function); 661 } 662 663 } // namespace detail 664 665 template <class RealType = double, class Policy = policies::policy<> > 666 class non_central_beta_distribution 667 { 668 public: 669 typedef RealType value_type; 670 typedef Policy policy_type; 671 non_central_beta_distribution(RealType a_,RealType b_,RealType lambda)672 non_central_beta_distribution(RealType a_, RealType b_, RealType lambda) : a(a_), b(b_), ncp(lambda) 673 { 674 const char* function = "boost::math::non_central_beta_distribution<%1%>::non_central_beta_distribution(%1%,%1%)"; 675 RealType r; 676 beta_detail::check_alpha( 677 function, 678 a, &r, Policy()); 679 beta_detail::check_beta( 680 function, 681 b, &r, Policy()); 682 detail::check_non_centrality( 683 function, 684 lambda, 685 &r, 686 Policy()); 687 } // non_central_beta_distribution constructor. 688 alpha() const689 RealType alpha() const 690 { // Private data getter function. 691 return a; 692 } beta() const693 RealType beta() const 694 { // Private data getter function. 695 return b; 696 } non_centrality() const697 RealType non_centrality() const 698 { // Private data getter function. 699 return ncp; 700 } 701 private: 702 // Data member, initialized by constructor. 703 RealType a; // alpha. 704 RealType b; // beta. 705 RealType ncp; // non-centrality parameter 706 }; // template <class RealType, class Policy> class non_central_beta_distribution 707 708 typedef non_central_beta_distribution<double> non_central_beta; // Reserved name of type double. 709 710 // Non-member functions to give properties of the distribution. 711 712 template <class RealType, class Policy> range(const non_central_beta_distribution<RealType,Policy> &)713 inline const std::pair<RealType, RealType> range(const non_central_beta_distribution<RealType, Policy>& /* dist */) 714 { // Range of permissible values for random variable k. 715 using boost::math::tools::max_value; 716 return std::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1)); 717 } 718 719 template <class RealType, class Policy> support(const non_central_beta_distribution<RealType,Policy> &)720 inline const std::pair<RealType, RealType> support(const non_central_beta_distribution<RealType, Policy>& /* dist */) 721 { // Range of supported values for random variable k. 722 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. 723 using boost::math::tools::max_value; 724 return std::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1)); 725 } 726 727 template <class RealType, class Policy> mode(const non_central_beta_distribution<RealType,Policy> & dist)728 inline RealType mode(const non_central_beta_distribution<RealType, Policy>& dist) 729 { // mode. 730 static const char* function = "mode(non_central_beta_distribution<%1%> const&)"; 731 732 RealType a = dist.alpha(); 733 RealType b = dist.beta(); 734 RealType l = dist.non_centrality(); 735 RealType r; 736 if(!beta_detail::check_alpha( 737 function, 738 a, &r, Policy()) 739 || 740 !beta_detail::check_beta( 741 function, 742 b, &r, Policy()) 743 || 744 !detail::check_non_centrality( 745 function, 746 l, 747 &r, 748 Policy())) 749 return (RealType)r; 750 RealType c = a + b + l / 2; 751 RealType mean = 1 - (b / c) * (1 + l / (2 * c * c)); 752 return detail::generic_find_mode_01( 753 dist, 754 mean, 755 function); 756 } 757 758 // 759 // We don't have the necessary information to implement 760 // these at present. These are just disabled for now, 761 // prototypes retained so we can fill in the blanks 762 // later: 763 // 764 template <class RealType, class Policy> mean(const non_central_beta_distribution<RealType,Policy> & dist)765 inline RealType mean(const non_central_beta_distribution<RealType, Policy>& dist) 766 { 767 BOOST_MATH_STD_USING 768 RealType a = dist.alpha(); 769 RealType b = dist.beta(); 770 RealType d = dist.non_centrality(); 771 RealType apb = a + b; 772 return exp(-d / 2) * a * detail::hypergeometric_2F2<RealType, Policy>(1 + a, apb, a, 1 + apb, d / 2, Policy()) / apb; 773 } // mean 774 775 template <class RealType, class Policy> variance(const non_central_beta_distribution<RealType,Policy> & dist)776 inline RealType variance(const non_central_beta_distribution<RealType, Policy>& dist) 777 { 778 // 779 // Relative error of this function may be arbitrarily large... absolute 780 // error will be small however... that's the best we can do for now. 781 // 782 BOOST_MATH_STD_USING 783 RealType a = dist.alpha(); 784 RealType b = dist.beta(); 785 RealType d = dist.non_centrality(); 786 RealType apb = a + b; 787 RealType result = detail::hypergeometric_2F2(RealType(1 + a), apb, a, RealType(1 + apb), RealType(d / 2), Policy()); 788 result *= result * -exp(-d) * a * a / (apb * apb); 789 result += exp(-d / 2) * a * (1 + a) * detail::hypergeometric_2F2(RealType(2 + a), apb, a, RealType(2 + apb), RealType(d / 2), Policy()) / (apb * (1 + apb)); 790 return result; 791 } 792 793 // RealType standard_deviation(const non_central_beta_distribution<RealType, Policy>& dist) 794 // standard_deviation provided by derived accessors. 795 template <class RealType, class Policy> skewness(const non_central_beta_distribution<RealType,Policy> &)796 inline RealType skewness(const non_central_beta_distribution<RealType, Policy>& /*dist*/) 797 { // skewness = sqrt(l). 798 const char* function = "boost::math::non_central_beta_distribution<%1%>::skewness()"; 799 typedef typename Policy::assert_undefined_type assert_type; 800 static_assert(assert_type::value == 0, "Assert type is undefined."); 801 802 return policies::raise_evaluation_error<RealType>( 803 function, 804 "This function is not yet implemented, the only sensible result is %1%.", 805 std::numeric_limits<RealType>::quiet_NaN(), Policy()); // infinity? 806 } 807 808 template <class RealType, class Policy> kurtosis_excess(const non_central_beta_distribution<RealType,Policy> &)809 inline RealType kurtosis_excess(const non_central_beta_distribution<RealType, Policy>& /*dist*/) 810 { 811 const char* function = "boost::math::non_central_beta_distribution<%1%>::kurtosis_excess()"; 812 typedef typename Policy::assert_undefined_type assert_type; 813 static_assert(assert_type::value == 0, "Assert type is undefined."); 814 815 return policies::raise_evaluation_error<RealType>( 816 function, 817 "This function is not yet implemented, the only sensible result is %1%.", 818 std::numeric_limits<RealType>::quiet_NaN(), Policy()); // infinity? 819 } // kurtosis_excess 820 821 template <class RealType, class Policy> kurtosis(const non_central_beta_distribution<RealType,Policy> & dist)822 inline RealType kurtosis(const non_central_beta_distribution<RealType, Policy>& dist) 823 { 824 return kurtosis_excess(dist) + 3; 825 } 826 827 template <class RealType, class Policy> pdf(const non_central_beta_distribution<RealType,Policy> & dist,const RealType & x)828 inline RealType pdf(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& x) 829 { // Probability Density/Mass Function. 830 return detail::nc_beta_pdf(dist, x); 831 } // pdf 832 833 template <class RealType, class Policy> 834 RealType cdf(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& x) 835 { 836 const char* function = "boost::math::non_central_beta_distribution<%1%>::cdf(%1%)"; 837 RealType a = dist.alpha(); 838 RealType b = dist.beta(); 839 RealType l = dist.non_centrality(); 840 RealType r; 841 if(!beta_detail::check_alpha( 842 function, 843 a, &r, Policy()) 844 || 845 !beta_detail::check_beta( 846 function, 847 b, &r, Policy()) 848 || 849 !detail::check_non_centrality( 850 function, 851 l, 852 &r, 853 Policy()) 854 || 855 !beta_detail::check_x( 856 function, 857 x, 858 &r, 859 Policy())) 860 return (RealType)r; 861 862 if(l == 0) 863 return cdf(beta_distribution<RealType, Policy>(a, b), x); 864 865 return detail::non_central_beta_cdf(x, RealType(1 - x), a, b, l, false, Policy()); 866 } // cdf 867 868 template <class RealType, class Policy> 869 RealType cdf(const complemented2_type<non_central_beta_distribution<RealType, Policy>, RealType>& c) 870 { // Complemented Cumulative Distribution Function 871 const char* function = "boost::math::non_central_beta_distribution<%1%>::cdf(%1%)"; 872 non_central_beta_distribution<RealType, Policy> const& dist = c.dist; 873 RealType a = dist.alpha(); 874 RealType b = dist.beta(); 875 RealType l = dist.non_centrality(); 876 RealType x = c.param; 877 RealType r; 878 if(!beta_detail::check_alpha( 879 function, 880 a, &r, Policy()) 881 || 882 !beta_detail::check_beta( 883 function, 884 b, &r, Policy()) 885 || 886 !detail::check_non_centrality( 887 function, 888 l, 889 &r, 890 Policy()) 891 || 892 !beta_detail::check_x( 893 function, 894 x, 895 &r, 896 Policy())) 897 return (RealType)r; 898 899 if(l == 0) 900 return cdf(complement(beta_distribution<RealType, Policy>(a, b), x)); 901 902 return detail::non_central_beta_cdf(x, RealType(1 - x), a, b, l, true, Policy()); 903 } // ccdf 904 905 template <class RealType, class Policy> quantile(const non_central_beta_distribution<RealType,Policy> & dist,const RealType & p)906 inline RealType quantile(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& p) 907 { // Quantile (or Percent Point) function. 908 return detail::nc_beta_quantile(dist, p, false); 909 } // quantile 910 911 template <class RealType, class Policy> quantile(const complemented2_type<non_central_beta_distribution<RealType,Policy>,RealType> & c)912 inline RealType quantile(const complemented2_type<non_central_beta_distribution<RealType, Policy>, RealType>& c) 913 { // Quantile (or Percent Point) function. 914 return detail::nc_beta_quantile(c.dist, c.param, true); 915 } // quantile complement. 916 917 } // namespace math 918 } // namespace boost 919 920 // This include must be at the end, *after* the accessors 921 // for this distribution have been defined, in order to 922 // keep compilers that support two-phase lookup happy. 923 #include <boost/math/distributions/detail/derived_accessors.hpp> 924 925 #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_BETA_HPP 926 927