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