1 // boost\math\distributions\non_central_t.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_T_HPP 11 #define BOOST_MATH_SPECIAL_NON_CENTRAL_T_HPP 12 13 #include <boost/math/distributions/fwd.hpp> 14 #include <boost/math/distributions/non_central_beta.hpp> // for nc beta 15 #include <boost/math/distributions/normal.hpp> // for normal CDF and quantile 16 #include <boost/math/distributions/students_t.hpp> 17 #include <boost/math/distributions/detail/generic_quantile.hpp> // quantile 18 19 namespace boost 20 { 21 namespace math 22 { 23 24 template <class RealType, class Policy> 25 class non_central_t_distribution; 26 27 namespace detail{ 28 29 template <class T, class Policy> 30 T non_central_t2_p(T n, T delta, T x, T y, const Policy& pol, T init_val) 31 { 32 BOOST_MATH_STD_USING 33 // 34 // Variables come first: 35 // 36 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); 37 T errtol = policies::get_epsilon<T, Policy>(); 38 T d2 = delta * delta / 2; 39 // 40 // k is the starting point for iteration, and is the 41 // maximum of the poisson weighting term: 42 // 43 int k = boost::math::itrunc(d2); 44 T pois; 45 if(k < 15) 46 { 47 // Since we'll likely need 30-40 terms anyway, start from zero 48 // since this simplifies the arithmetic, don't go too overboard though 49 // as this is the *unstable* direction: 50 k = 0; 51 // Starting Poisson weight: 52 pois = exp(-d2) * 2 / constants::root_pi<T>(); 53 pois *= delta / constants::root_two<T>(); 54 } 55 else 56 { 57 // Starting Poisson weight: 58 pois = gamma_p_derivative(T(k+1), d2, pol) 59 * tgamma_delta_ratio(T(k + 1), T(0.5f)) 60 * delta / constants::root_two<T>(); 61 } 62 if(pois == 0) 63 return init_val; 64 T xterm, beta; 65 // Recurrance & starting beta terms: 66 if(k == 0) 67 { 68 beta = -boost::math::powm1(y, n / 2, pol); 69 xterm = beta > 0.5f ? T(pow(y, n / 2)) : T(1 - beta); 70 } 71 else 72 { 73 beta = x < y 74 ? detail::ibeta_imp(T(k + 1), T(n / 2), x, pol, false, true, &xterm) 75 : detail::ibeta_imp(T(n / 2), T(k + 1), y, pol, true, true, &xterm); 76 xterm *= y / (n / 2 + k); 77 } 78 T poisf(pois), betaf(beta), xtermf(xterm); 79 T sum = init_val; 80 if((xterm == 0) && (beta == 0)) 81 return init_val; 82 83 // 84 // Backwards recursion first, this is the stable 85 // direction for recursion: 86 // 87 boost::uintmax_t count = 0; 88 for(int i = k; i >= 0; --i) 89 { 90 T term = beta * pois; 91 sum += term; 92 if(fabs(term/sum) < errtol) 93 break; 94 pois *= (i + 0.5f) / d2; 95 beta += xterm; 96 xterm *= (i) / (x * (n / 2 + i - 1)); 97 ++count; 98 } 99 for(int i = k + 1; ; ++i) 100 { 101 poisf *= d2 / (i + 0.5f); 102 xtermf *= (x * (n / 2 + i - 1)) / (i); 103 betaf -= xtermf; 104 T term = poisf * betaf; 105 sum += term; 106 if(fabs(term/sum) < errtol) 107 break; 108 ++count; 109 if(count > max_iter) 110 { 111 return policies::raise_evaluation_error( 112 "cdf(non_central_t_distribution<%1%>, %1%)", 113 "Series did not converge, closest value was %1%", sum, pol); 114 } 115 } 116 return sum; 117 } 118 119 template <class T, class Policy> 120 T non_central_t2_q(T n, T delta, T x, T y, const Policy& pol, T init_val) 121 { 122 BOOST_MATH_STD_USING 123 // 124 // Variables come first: 125 // 126 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); 127 T errtol = boost::math::policies::get_epsilon<T, Policy>(); 128 T d2 = delta * delta / 2; 129 // 130 // k is the starting point for iteration, and is the 131 // maximum of the poisson weighting term: 132 // 133 int k = boost::math::itrunc(d2); 134 if(k < 30) 135 { 136 // We typically need around 40 terms so may as well start at 0 137 // and gain faster computation of starting conditions: 138 k = 0; 139 } 140 // Starting Poisson weight: 141 T pois; 142 if(k == 0) 143 { 144 pois = exp(-d2) * 2 / constants::root_pi<T>(); 145 pois *= delta / constants::root_two<T>(); 146 } 147 else if((k < (int)(max_factorial<T>::value)) && (d2 < tools::log_max_value<T>()) && (log(d2) * k < tools::log_max_value<T>())) 148 { 149 // 150 // For small k we can optimise this calculation by using 151 // a simpler reduced formula: 152 // 153 pois = exp(-d2); 154 pois *= pow(d2, static_cast<T>(k)); 155 pois /= boost::math::tgamma(T(k + 1 + 0.5), pol); 156 pois *= delta / constants::root_two<T>(); 157 } 158 else 159 { 160 pois = gamma_p_derivative(T(k+1), d2, pol) 161 * tgamma_delta_ratio(T(k + 1), T(0.5f)) 162 * delta / constants::root_two<T>(); 163 } 164 if(pois == 0) 165 return init_val; 166 // Recurance term: 167 T xterm; 168 T beta; 169 // Starting beta term: 170 if(k != 0) 171 { 172 beta = x < y 173 ? detail::ibeta_imp(T(k + 1), T(n / 2), x, pol, true, true, &xterm) 174 : detail::ibeta_imp(T(n / 2), T(k + 1), y, pol, false, true, &xterm); 175 176 xterm *= y / (n / 2 + k); 177 } 178 else 179 { 180 beta = pow(y, n / 2); 181 xterm = beta; 182 } 183 T poisf(pois), betaf(beta), xtermf(xterm); 184 T sum = init_val; 185 if((xterm == 0) && (beta == 0)) 186 return init_val; 187 188 // 189 // Fused forward and backwards recursion: 190 // 191 boost::uintmax_t count = 0; 192 for(int i = k + 1, j = k; ; ++i, --j) 193 { 194 poisf *= d2 / (i + 0.5f); 195 xtermf *= (x * (n / 2 + i - 1)) / (i); 196 betaf += xtermf; 197 T term = poisf * betaf; 198 199 if(j >= 0) 200 { 201 term += beta * pois; 202 pois *= (j + 0.5f) / d2; 203 beta -= xterm; 204 xterm *= (j) / (x * (n / 2 + j - 1)); 205 } 206 207 sum += term; 208 if(fabs(term/sum) < errtol) 209 break; 210 if(count > max_iter) 211 { 212 return policies::raise_evaluation_error( 213 "cdf(non_central_t_distribution<%1%>, %1%)", 214 "Series did not converge, closest value was %1%", sum, pol); 215 } 216 ++count; 217 } 218 return sum; 219 } 220 221 template <class T, class Policy> 222 T non_central_t_cdf(T n, T delta, T t, bool invert, const Policy& pol) 223 { 224 // 225 // For t < 0 we have to use reflect: 226 // 227 if(t < 0) 228 { 229 t = -t; 230 delta = -delta; 231 invert = !invert; 232 } 233 // 234 // x and y are the corresponding random 235 // variables for the noncentral beta distribution, 236 // with y = 1 - x: 237 // 238 T x = t * t / (n + t * t); 239 T y = n / (n + t * t); 240 T d2 = delta * delta; 241 T a = 0.5f; 242 T b = n / 2; 243 T c = a + b + d2 / 2; 244 // 245 // Crossover point for calculating p or q is the same 246 // as for the noncentral beta: 247 // 248 T cross = 1 - (b / c) * (1 + d2 / (2 * c * c)); 249 T result; 250 if(x < cross) 251 { 252 // 253 // Calculate p: 254 // 255 if(x != 0) 256 { 257 result = non_central_beta_p(a, b, d2, x, y, pol); 258 result = non_central_t2_p(n, delta, x, y, pol, result); 259 result /= 2; 260 } 261 else 262 result = 0; 263 result += cdf(boost::math::normal_distribution<T, Policy>(), -delta); 264 } 265 else 266 { 267 // 268 // Calculate q: 269 // 270 invert = !invert; 271 if(x != 0) 272 { 273 result = non_central_beta_q(a, b, d2, x, y, pol); 274 result = non_central_t2_q(n, delta, x, y, pol, result); 275 result /= 2; 276 } 277 else 278 result = cdf(complement(boost::math::normal_distribution<T, Policy>(), -delta)); 279 } 280 if(invert) 281 result = 1 - result; 282 return result; 283 } 284 285 template <class T, class Policy> 286 T non_central_t_quantile(T v, T delta, T p, T q, const Policy&) 287 { 288 BOOST_MATH_STD_USING 289 static const char* function = "quantile(non_central_t_distribution<%1%>, %1%)"; 290 typedef typename policies::evaluation<T, Policy>::type value_type; 291 typedef typename policies::normalise< 292 Policy, 293 policies::promote_float<false>, 294 policies::promote_double<false>, 295 policies::discrete_quantile<>, 296 policies::assert_undefined<> >::type forwarding_policy; 297 298 T r; 299 if(!detail::check_df( 300 function, 301 v, &r, Policy()) 302 || 303 !detail::check_finite( 304 function, 305 delta, 306 &r, 307 Policy()) 308 || 309 !detail::check_probability( 310 function, 311 p, 312 &r, 313 Policy())) 314 return r; 315 316 value_type guess = 0; 317 if(v > 3) 318 { 319 value_type mean = delta * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, T(0.5f)); 320 value_type var = ((delta * delta + 1) * v) / (v - 2) - mean * mean; 321 if(p < q) 322 guess = quantile(normal_distribution<value_type, forwarding_policy>(mean, var), p); 323 else 324 guess = quantile(complement(normal_distribution<value_type, forwarding_policy>(mean, var), q)); 325 } 326 // 327 // We *must* get the sign of the initial guess correct, 328 // or our root-finder will fail, so double check it now: 329 // 330 value_type pzero = non_central_t_cdf( 331 static_cast<value_type>(v), 332 static_cast<value_type>(delta), 333 static_cast<value_type>(0), 334 !(p < q), 335 forwarding_policy()); 336 int s; 337 if(p < q) 338 s = boost::math::sign(p - pzero); 339 else 340 s = boost::math::sign(pzero - q); 341 if(s != boost::math::sign(guess)) 342 { 343 guess = s; 344 } 345 346 value_type result = detail::generic_quantile( 347 non_central_t_distribution<value_type, forwarding_policy>(v, delta), 348 (p < q ? p : q), 349 guess, 350 (p >= q), 351 function); 352 return policies::checked_narrowing_cast<T, forwarding_policy>( 353 result, 354 function); 355 } 356 357 template <class T, class Policy> 358 T non_central_t2_pdf(T n, T delta, T x, T y, const Policy& pol, T init_val) 359 { 360 BOOST_MATH_STD_USING 361 // 362 // Variables come first: 363 // 364 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); 365 T errtol = boost::math::policies::get_epsilon<T, Policy>(); 366 T d2 = delta * delta / 2; 367 // 368 // k is the starting point for iteration, and is the 369 // maximum of the poisson weighting term: 370 // 371 int k = boost::math::itrunc(d2); 372 T pois, xterm; 373 if(k < 30) 374 { 375 // 376 // Since we'll need at least 30-40 terms anyway, start from 0 377 // since this simplifies the starting arithmetic: 378 // 379 k = 0; 380 // Starting Poisson weight: 381 pois = exp(-d2) 382 * (2 / constants::root_pi<T>()) 383 * delta / constants::root_two<T>(); 384 // Starting beta term: 385 xterm = pow(y, n / 2 - 1) * n / 2; 386 } 387 else 388 { 389 // Starting Poisson weight: 390 pois = gamma_p_derivative(T(k+1), d2, pol) 391 * tgamma_delta_ratio(T(k + 1), T(0.5f)) 392 * delta / constants::root_two<T>(); 393 // Starting beta term: 394 xterm = x < y 395 ? ibeta_derivative(T(k + 1), n / 2, x, pol) 396 : ibeta_derivative(n / 2, T(k + 1), y, pol); 397 } 398 T poisf(pois), xtermf(xterm); 399 T sum = init_val; 400 if((pois == 0) || (xterm == 0)) 401 return init_val; 402 403 // 404 // Backwards recursion first, this is the stable 405 // direction for recursion: 406 // 407 boost::uintmax_t count = 0; 408 for(int i = k; i >= 0; --i) 409 { 410 T term = xterm * pois; 411 sum += term; 412 if((fabs(term/sum) < errtol) || (term == 0)) 413 break; 414 pois *= (i + 0.5f) / d2; 415 xterm *= (i) / (x * (n / 2 + i)); 416 ++count; 417 if(count > max_iter) 418 { 419 return policies::raise_evaluation_error( 420 "pdf(non_central_t_distribution<%1%>, %1%)", 421 "Series did not converge, closest value was %1%", sum, pol); 422 } 423 } 424 for(int i = k + 1; ; ++i) 425 { 426 poisf *= d2 / (i + 0.5f); 427 xtermf *= (x * (n / 2 + i)) / (i); 428 T term = poisf * xtermf; 429 sum += term; 430 if((fabs(term/sum) < errtol) || (term == 0)) 431 break; 432 ++count; 433 if(count > max_iter) 434 { 435 return policies::raise_evaluation_error( 436 "pdf(non_central_t_distribution<%1%>, %1%)", 437 "Series did not converge, closest value was %1%", sum, pol); 438 } 439 } 440 return sum; 441 } 442 443 template <class T, class Policy> 444 T non_central_t_pdf(T n, T delta, T t, const Policy& pol) 445 { 446 BOOST_MATH_STD_USING 447 // 448 // For t < 0 we have to use the reflection formula: 449 // 450 if(t < 0) 451 { 452 t = -t; 453 delta = -delta; 454 } 455 if(t == 0) 456 { 457 // 458 // Handle this as a special case, using the formula 459 // from Weisstein, Eric W. 460 // "Noncentral Student's t-Distribution." 461 // From MathWorld--A Wolfram Web Resource. 462 // http://mathworld.wolfram.com/NoncentralStudentst-Distribution.html 463 // 464 // The formula is simplified thanks to the relation 465 // 1F1(a,b,0) = 1. 466 // 467 return tgamma_delta_ratio(n / 2 + 0.5f, T(0.5f)) 468 * sqrt(n / constants::pi<T>()) 469 * exp(-delta * delta / 2) / 2; 470 } 471 // 472 // x and y are the corresponding random 473 // variables for the noncentral beta distribution, 474 // with y = 1 - x: 475 // 476 T x = t * t / (n + t * t); 477 T y = n / (n + t * t); 478 T a = 0.5f; 479 T b = n / 2; 480 T d2 = delta * delta; 481 // 482 // Calculate pdf: 483 // 484 T dt = n * t / (n * n + 2 * n * t * t + t * t * t * t); 485 T result = non_central_beta_pdf(a, b, d2, x, y, pol); 486 T tol = tools::epsilon<T>() * result * 500; 487 result = non_central_t2_pdf(n, delta, x, y, pol, result); 488 if(result <= tol) 489 result = 0; 490 result *= dt; 491 return result; 492 } 493 494 template <class T, class Policy> 495 T mean(T v, T delta, const Policy& pol) 496 { 497 BOOST_MATH_STD_USING 498 return delta * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, T(0.5f), pol); 499 } 500 501 template <class T, class Policy> 502 T variance(T v, T delta, const Policy& pol) 503 { 504 T result = ((delta * delta + 1) * v) / (v - 2); 505 T m = mean(v, delta, pol); 506 result -= m * m; 507 return result; 508 } 509 510 template <class T, class Policy> 511 T skewness(T v, T delta, const Policy& pol) 512 { 513 BOOST_MATH_STD_USING 514 T mean = boost::math::detail::mean(v, delta, pol); 515 T l2 = delta * delta; 516 T var = ((l2 + 1) * v) / (v - 2) - mean * mean; 517 T result = -2 * var; 518 result += v * (l2 + 2 * v - 3) / ((v - 3) * (v - 2)); 519 result *= mean; 520 result /= pow(var, T(1.5f)); 521 return result; 522 } 523 524 template <class T, class Policy> 525 T kurtosis_excess(T v, T delta, const Policy& pol) 526 { 527 BOOST_MATH_STD_USING 528 T mean = boost::math::detail::mean(v, delta, pol); 529 T l2 = delta * delta; 530 T var = ((l2 + 1) * v) / (v - 2) - mean * mean; 531 T result = -3 * var; 532 result += v * (l2 * (v + 1) + 3 * (3 * v - 5)) / ((v - 3) * (v - 2)); 533 result *= -mean * mean; 534 result += v * v * (l2 * l2 + 6 * l2 + 3) / ((v - 4) * (v - 2)); 535 result /= var * var; 536 return result; 537 } 538 539 #if 0 540 // 541 // This code is disabled, since there can be multiple answers to the 542 // question, and it's not clear how to find the "right" one. 543 // 544 template <class RealType, class Policy> 545 struct t_degrees_of_freedom_finder 546 { 547 t_degrees_of_freedom_finder( 548 RealType delta_, RealType x_, RealType p_, bool c) 549 : delta(delta_), x(x_), p(p_), comp(c) {} 550 551 RealType operator()(const RealType& v) 552 { 553 non_central_t_distribution<RealType, Policy> d(v, delta); 554 return comp ? 555 p - cdf(complement(d, x)) 556 : cdf(d, x) - p; 557 } 558 private: 559 RealType delta; 560 RealType x; 561 RealType p; 562 bool comp; 563 }; 564 565 template <class RealType, class Policy> 566 inline RealType find_t_degrees_of_freedom( 567 RealType delta, RealType x, RealType p, RealType q, const Policy& pol) 568 { 569 const char* function = "non_central_t<%1%>::find_degrees_of_freedom"; 570 if((p == 0) || (q == 0)) 571 { 572 // 573 // Can't a thing if one of p and q is zero: 574 // 575 return policies::raise_evaluation_error<RealType>(function, 576 "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%", 577 RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); 578 } 579 t_degrees_of_freedom_finder<RealType, Policy> f(delta, x, p < q ? p : q, p < q ? false : true); 580 tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>()); 581 boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); 582 // 583 // Pick an initial guess: 584 // 585 RealType guess = 200; 586 std::pair<RealType, RealType> ir = tools::bracket_and_solve_root( 587 f, guess, RealType(2), false, tol, max_iter, pol); 588 RealType result = ir.first + (ir.second - ir.first) / 2; 589 if(max_iter >= policies::get_max_root_iterations<Policy>()) 590 { 591 policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" 592 " or there is no answer to problem. Current best guess is %1%", result, Policy()); 593 } 594 return result; 595 } 596 597 template <class RealType, class Policy> 598 struct t_non_centrality_finder 599 { 600 t_non_centrality_finder( 601 RealType v_, RealType x_, RealType p_, bool c) 602 : v(v_), x(x_), p(p_), comp(c) {} 603 604 RealType operator()(const RealType& delta) 605 { 606 non_central_t_distribution<RealType, Policy> d(v, delta); 607 return comp ? 608 p - cdf(complement(d, x)) 609 : cdf(d, x) - p; 610 } 611 private: 612 RealType v; 613 RealType x; 614 RealType p; 615 bool comp; 616 }; 617 618 template <class RealType, class Policy> 619 inline RealType find_t_non_centrality( 620 RealType v, RealType x, RealType p, RealType q, const Policy& pol) 621 { 622 const char* function = "non_central_t<%1%>::find_t_non_centrality"; 623 if((p == 0) || (q == 0)) 624 { 625 // 626 // Can't do a thing if one of p and q is zero: 627 // 628 return policies::raise_evaluation_error<RealType>(function, 629 "Can't find non centrality parameter when the probability is 0 or 1, only possible answer is %1%", 630 RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); 631 } 632 t_non_centrality_finder<RealType, Policy> f(v, x, p < q ? p : q, p < q ? false : true); 633 tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>()); 634 boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); 635 // 636 // Pick an initial guess that we know is the right side of 637 // zero: 638 // 639 RealType guess; 640 if(f(0) < 0) 641 guess = 1; 642 else 643 guess = -1; 644 std::pair<RealType, RealType> ir = tools::bracket_and_solve_root( 645 f, guess, RealType(2), false, tol, max_iter, pol); 646 RealType result = ir.first + (ir.second - ir.first) / 2; 647 if(max_iter >= policies::get_max_root_iterations<Policy>()) 648 { 649 policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" 650 " or there is no answer to problem. Current best guess is %1%", result, Policy()); 651 } 652 return result; 653 } 654 #endif 655 } // namespace detail 656 657 template <class RealType = double, class Policy = policies::policy<> > 658 class non_central_t_distribution 659 { 660 public: 661 typedef RealType value_type; 662 typedef Policy policy_type; 663 non_central_t_distribution(RealType v_,RealType lambda)664 non_central_t_distribution(RealType v_, RealType lambda) : v(v_), ncp(lambda) 665 { 666 const char* function = "boost::math::non_central_t_distribution<%1%>::non_central_t_distribution(%1%,%1%)"; 667 RealType r; 668 detail::check_df( 669 function, 670 v, &r, Policy()); 671 detail::check_finite( 672 function, 673 lambda, 674 &r, 675 Policy()); 676 } // non_central_t_distribution constructor. 677 degrees_of_freedom() const678 RealType degrees_of_freedom() const 679 { // Private data getter function. 680 return v; 681 } non_centrality() const682 RealType non_centrality() const 683 { // Private data getter function. 684 return ncp; 685 } 686 #if 0 687 // 688 // This code is disabled, since there can be multiple answers to the 689 // question, and it's not clear how to find the "right" one. 690 // 691 static RealType find_degrees_of_freedom(RealType delta, RealType x, RealType p) 692 { 693 const char* function = "non_central_t<%1%>::find_degrees_of_freedom"; 694 typedef typename policies::evaluation<RealType, Policy>::type value_type; 695 typedef typename policies::normalise< 696 Policy, 697 policies::promote_float<false>, 698 policies::promote_double<false>, 699 policies::discrete_quantile<>, 700 policies::assert_undefined<> >::type forwarding_policy; 701 value_type result = detail::find_t_degrees_of_freedom( 702 static_cast<value_type>(delta), 703 static_cast<value_type>(x), 704 static_cast<value_type>(p), 705 static_cast<value_type>(1-p), 706 forwarding_policy()); 707 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 708 result, 709 function); 710 } 711 template <class A, class B, class C> 712 static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c) 713 { 714 const char* function = "non_central_t<%1%>::find_degrees_of_freedom"; 715 typedef typename policies::evaluation<RealType, Policy>::type value_type; 716 typedef typename policies::normalise< 717 Policy, 718 policies::promote_float<false>, 719 policies::promote_double<false>, 720 policies::discrete_quantile<>, 721 policies::assert_undefined<> >::type forwarding_policy; 722 value_type result = detail::find_t_degrees_of_freedom( 723 static_cast<value_type>(c.dist), 724 static_cast<value_type>(c.param1), 725 static_cast<value_type>(1-c.param2), 726 static_cast<value_type>(c.param2), 727 forwarding_policy()); 728 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 729 result, 730 function); 731 } 732 static RealType find_non_centrality(RealType v, RealType x, RealType p) 733 { 734 const char* function = "non_central_t<%1%>::find_t_non_centrality"; 735 typedef typename policies::evaluation<RealType, Policy>::type value_type; 736 typedef typename policies::normalise< 737 Policy, 738 policies::promote_float<false>, 739 policies::promote_double<false>, 740 policies::discrete_quantile<>, 741 policies::assert_undefined<> >::type forwarding_policy; 742 value_type result = detail::find_t_non_centrality( 743 static_cast<value_type>(v), 744 static_cast<value_type>(x), 745 static_cast<value_type>(p), 746 static_cast<value_type>(1-p), 747 forwarding_policy()); 748 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 749 result, 750 function); 751 } 752 template <class A, class B, class C> 753 static RealType find_non_centrality(const complemented3_type<A,B,C>& c) 754 { 755 const char* function = "non_central_t<%1%>::find_t_non_centrality"; 756 typedef typename policies::evaluation<RealType, Policy>::type value_type; 757 typedef typename policies::normalise< 758 Policy, 759 policies::promote_float<false>, 760 policies::promote_double<false>, 761 policies::discrete_quantile<>, 762 policies::assert_undefined<> >::type forwarding_policy; 763 value_type result = detail::find_t_non_centrality( 764 static_cast<value_type>(c.dist), 765 static_cast<value_type>(c.param1), 766 static_cast<value_type>(1-c.param2), 767 static_cast<value_type>(c.param2), 768 forwarding_policy()); 769 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 770 result, 771 function); 772 } 773 #endif 774 private: 775 // Data member, initialized by constructor. 776 RealType v; // degrees of freedom 777 RealType ncp; // non-centrality parameter 778 }; // template <class RealType, class Policy> class non_central_t_distribution 779 780 typedef non_central_t_distribution<double> non_central_t; // Reserved name of type double. 781 782 // Non-member functions to give properties of the distribution. 783 784 template <class RealType, class Policy> range(const non_central_t_distribution<RealType,Policy> &)785 inline const std::pair<RealType, RealType> range(const non_central_t_distribution<RealType, Policy>& /* dist */) 786 { // Range of permissible values for random variable k. 787 using boost::math::tools::max_value; 788 return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); 789 } 790 791 template <class RealType, class Policy> support(const non_central_t_distribution<RealType,Policy> &)792 inline const std::pair<RealType, RealType> support(const non_central_t_distribution<RealType, Policy>& /* dist */) 793 { // Range of supported values for random variable k. 794 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. 795 using boost::math::tools::max_value; 796 return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); 797 } 798 799 template <class RealType, class Policy> mode(const non_central_t_distribution<RealType,Policy> & dist)800 inline RealType mode(const non_central_t_distribution<RealType, Policy>& dist) 801 { // mode. 802 static const char* function = "mode(non_central_t_distribution<%1%> const&)"; 803 RealType v = dist.degrees_of_freedom(); 804 RealType l = dist.non_centrality(); 805 RealType r; 806 if(!detail::check_df( 807 function, 808 v, &r, Policy()) 809 || 810 !detail::check_finite( 811 function, 812 l, 813 &r, 814 Policy())) 815 return (RealType)r; 816 817 BOOST_MATH_STD_USING 818 819 RealType m = v < 3 ? 0 : detail::mean(v, l, Policy()); 820 RealType var = v < 4 ? 1 : detail::variance(v, l, Policy()); 821 822 return detail::generic_find_mode( 823 dist, 824 m, 825 function, 826 sqrt(var)); 827 } 828 829 template <class RealType, class Policy> mean(const non_central_t_distribution<RealType,Policy> & dist)830 inline RealType mean(const non_central_t_distribution<RealType, Policy>& dist) 831 { 832 BOOST_MATH_STD_USING 833 const char* function = "mean(const non_central_t_distribution<%1%>&)"; 834 typedef typename policies::evaluation<RealType, Policy>::type value_type; 835 typedef typename policies::normalise< 836 Policy, 837 policies::promote_float<false>, 838 policies::promote_double<false>, 839 policies::discrete_quantile<>, 840 policies::assert_undefined<> >::type forwarding_policy; 841 RealType v = dist.degrees_of_freedom(); 842 RealType l = dist.non_centrality(); 843 RealType r; 844 if(!detail::check_df( 845 function, 846 v, &r, Policy()) 847 || 848 !detail::check_finite( 849 function, 850 l, 851 &r, 852 Policy())) 853 return (RealType)r; 854 if(v <= 1) 855 return policies::raise_domain_error<RealType>( 856 function, 857 "The non central t distribution has no defined mean for degrees of freedom <= 1: got v=%1%.", v, Policy()); 858 // return l * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, RealType(0.5f)); 859 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 860 detail::mean(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function); 861 862 } // mean 863 864 template <class RealType, class Policy> variance(const non_central_t_distribution<RealType,Policy> & dist)865 inline RealType variance(const non_central_t_distribution<RealType, Policy>& dist) 866 { // variance. 867 const char* function = "variance(const non_central_t_distribution<%1%>&)"; 868 typedef typename policies::evaluation<RealType, Policy>::type value_type; 869 typedef typename policies::normalise< 870 Policy, 871 policies::promote_float<false>, 872 policies::promote_double<false>, 873 policies::discrete_quantile<>, 874 policies::assert_undefined<> >::type forwarding_policy; 875 BOOST_MATH_STD_USING 876 RealType v = dist.degrees_of_freedom(); 877 RealType l = dist.non_centrality(); 878 RealType r; 879 if(!detail::check_df( 880 function, 881 v, &r, Policy()) 882 || 883 !detail::check_finite( 884 function, 885 l, 886 &r, 887 Policy())) 888 return (RealType)r; 889 if(v <= 2) 890 return policies::raise_domain_error<RealType>( 891 function, 892 "The non central t distribution has no defined variance for degrees of freedom <= 2: got v=%1%.", v, Policy()); 893 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 894 detail::variance(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function); 895 } 896 897 // RealType standard_deviation(const non_central_t_distribution<RealType, Policy>& dist) 898 // standard_deviation provided by derived accessors. 899 900 template <class RealType, class Policy> skewness(const non_central_t_distribution<RealType,Policy> & dist)901 inline RealType skewness(const non_central_t_distribution<RealType, Policy>& dist) 902 { // skewness = sqrt(l). 903 const char* function = "skewness(const non_central_t_distribution<%1%>&)"; 904 typedef typename policies::evaluation<RealType, Policy>::type value_type; 905 typedef typename policies::normalise< 906 Policy, 907 policies::promote_float<false>, 908 policies::promote_double<false>, 909 policies::discrete_quantile<>, 910 policies::assert_undefined<> >::type forwarding_policy; 911 RealType v = dist.degrees_of_freedom(); 912 RealType l = dist.non_centrality(); 913 RealType r; 914 if(!detail::check_df( 915 function, 916 v, &r, Policy()) 917 || 918 !detail::check_finite( 919 function, 920 l, 921 &r, 922 Policy())) 923 return (RealType)r; 924 if(v <= 3) 925 return policies::raise_domain_error<RealType>( 926 function, 927 "The non central t distribution has no defined skewness for degrees of freedom <= 3: got v=%1%.", v, Policy());; 928 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 929 detail::skewness(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function); 930 } 931 932 template <class RealType, class Policy> kurtosis_excess(const non_central_t_distribution<RealType,Policy> & dist)933 inline RealType kurtosis_excess(const non_central_t_distribution<RealType, Policy>& dist) 934 { 935 const char* function = "kurtosis_excess(const non_central_t_distribution<%1%>&)"; 936 typedef typename policies::evaluation<RealType, Policy>::type value_type; 937 typedef typename policies::normalise< 938 Policy, 939 policies::promote_float<false>, 940 policies::promote_double<false>, 941 policies::discrete_quantile<>, 942 policies::assert_undefined<> >::type forwarding_policy; 943 RealType v = dist.degrees_of_freedom(); 944 RealType l = dist.non_centrality(); 945 RealType r; 946 if(!detail::check_df( 947 function, 948 v, &r, Policy()) 949 || 950 !detail::check_finite( 951 function, 952 l, 953 &r, 954 Policy())) 955 return (RealType)r; 956 if(v <= 4) 957 return policies::raise_domain_error<RealType>( 958 function, 959 "The non central t distribution has no defined kurtosis for degrees of freedom <= 4: got v=%1%.", v, Policy());; 960 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 961 detail::kurtosis_excess(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function); 962 } // kurtosis_excess 963 964 template <class RealType, class Policy> kurtosis(const non_central_t_distribution<RealType,Policy> & dist)965 inline RealType kurtosis(const non_central_t_distribution<RealType, Policy>& dist) 966 { 967 return kurtosis_excess(dist) + 3; 968 } 969 970 template <class RealType, class Policy> pdf(const non_central_t_distribution<RealType,Policy> & dist,const RealType & t)971 inline RealType pdf(const non_central_t_distribution<RealType, Policy>& dist, const RealType& t) 972 { // Probability Density/Mass Function. 973 const char* function = "cdf(non_central_t_distribution<%1%>, %1%)"; 974 typedef typename policies::evaluation<RealType, Policy>::type value_type; 975 typedef typename policies::normalise< 976 Policy, 977 policies::promote_float<false>, 978 policies::promote_double<false>, 979 policies::discrete_quantile<>, 980 policies::assert_undefined<> >::type forwarding_policy; 981 982 RealType v = dist.degrees_of_freedom(); 983 RealType l = dist.non_centrality(); 984 RealType r; 985 if(!detail::check_df( 986 function, 987 v, &r, Policy()) 988 || 989 !detail::check_finite( 990 function, 991 l, 992 &r, 993 Policy()) 994 || 995 !detail::check_x( 996 function, 997 t, 998 &r, 999 Policy())) 1000 return (RealType)r; 1001 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 1002 detail::non_central_t_pdf(static_cast<value_type>(v), 1003 static_cast<value_type>(l), 1004 static_cast<value_type>(t), 1005 Policy()), 1006 function); 1007 } // pdf 1008 1009 template <class RealType, class Policy> 1010 RealType cdf(const non_central_t_distribution<RealType, Policy>& dist, const RealType& x) 1011 { 1012 const char* function = "boost::math::non_central_t_distribution<%1%>::cdf(%1%)"; 1013 typedef typename policies::evaluation<RealType, Policy>::type value_type; 1014 typedef typename policies::normalise< 1015 Policy, 1016 policies::promote_float<false>, 1017 policies::promote_double<false>, 1018 policies::discrete_quantile<>, 1019 policies::assert_undefined<> >::type forwarding_policy; 1020 1021 RealType v = dist.degrees_of_freedom(); 1022 RealType l = dist.non_centrality(); 1023 RealType r; 1024 if(!detail::check_df( 1025 function, 1026 v, &r, Policy()) 1027 || 1028 !detail::check_finite( 1029 function, 1030 l, 1031 &r, 1032 Policy()) 1033 || 1034 !detail::check_x( 1035 function, 1036 x, 1037 &r, 1038 Policy())) 1039 return (RealType)r; 1040 1041 if(l == 0) 1042 return cdf(students_t_distribution<RealType, Policy>(v), x); 1043 1044 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 1045 detail::non_central_t_cdf( 1046 static_cast<value_type>(v), 1047 static_cast<value_type>(l), 1048 static_cast<value_type>(x), 1049 false, Policy()), 1050 function); 1051 } // cdf 1052 1053 template <class RealType, class Policy> 1054 RealType cdf(const complemented2_type<non_central_t_distribution<RealType, Policy>, RealType>& c) 1055 { // Complemented Cumulative Distribution Function 1056 const char* function = "boost::math::non_central_t_distribution<%1%>::cdf(%1%)"; 1057 typedef typename policies::evaluation<RealType, Policy>::type value_type; 1058 typedef typename policies::normalise< 1059 Policy, 1060 policies::promote_float<false>, 1061 policies::promote_double<false>, 1062 policies::discrete_quantile<>, 1063 policies::assert_undefined<> >::type forwarding_policy; 1064 1065 non_central_t_distribution<RealType, Policy> const& dist = c.dist; 1066 RealType x = c.param; 1067 RealType v = dist.degrees_of_freedom(); 1068 RealType l = dist.non_centrality(); 1069 RealType r; 1070 if(!detail::check_df( 1071 function, 1072 v, &r, Policy()) 1073 || 1074 !detail::check_finite( 1075 function, 1076 l, 1077 &r, 1078 Policy()) 1079 || 1080 !detail::check_x( 1081 function, 1082 x, 1083 &r, 1084 Policy())) 1085 return (RealType)r; 1086 1087 if(l == 0) 1088 return cdf(complement(students_t_distribution<RealType, Policy>(v), x)); 1089 1090 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 1091 detail::non_central_t_cdf( 1092 static_cast<value_type>(v), 1093 static_cast<value_type>(l), 1094 static_cast<value_type>(x), 1095 true, Policy()), 1096 function); 1097 } // ccdf 1098 1099 template <class RealType, class Policy> quantile(const non_central_t_distribution<RealType,Policy> & dist,const RealType & p)1100 inline RealType quantile(const non_central_t_distribution<RealType, Policy>& dist, const RealType& p) 1101 { // Quantile (or Percent Point) function. 1102 RealType v = dist.degrees_of_freedom(); 1103 RealType l = dist.non_centrality(); 1104 return detail::non_central_t_quantile(v, l, p, RealType(1-p), Policy()); 1105 } // quantile 1106 1107 template <class RealType, class Policy> quantile(const complemented2_type<non_central_t_distribution<RealType,Policy>,RealType> & c)1108 inline RealType quantile(const complemented2_type<non_central_t_distribution<RealType, Policy>, RealType>& c) 1109 { // Quantile (or Percent Point) function. 1110 non_central_t_distribution<RealType, Policy> const& dist = c.dist; 1111 RealType q = c.param; 1112 RealType v = dist.degrees_of_freedom(); 1113 RealType l = dist.non_centrality(); 1114 return detail::non_central_t_quantile(v, l, RealType(1-q), q, Policy()); 1115 } // quantile complement. 1116 1117 } // namespace math 1118 } // namespace boost 1119 1120 // This include must be at the end, *after* the accessors 1121 // for this distribution have been defined, in order to 1122 // keep compilers that support two-phase lookup happy. 1123 #include <boost/math/distributions/detail/derived_accessors.hpp> 1124 1125 #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_T_HPP 1126 1127