boost/math/special_functions/detail/t_distribution_inv.hpp

   1 //  Copyright John Maddock 2007.
   2 //  Copyright Paul A. Bristow 2007
   3 //  Use, modification and distribution are subject to the
   4 //  Boost Software License, Version 1.0. (See accompanying file
   5 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #ifndef BOOST_MATH_SF_DETAIL_INV_T_HPP
   8 #define BOOST_MATH_SF_DETAIL_INV_T_HPP
   9
  10 #ifdef _MSC_VER
  11 #pragma once
  12 #endif
  13
  14 #include <boost/math/special_functions/cbrt.hpp>
  15 #include <boost/math/special_functions/round.hpp>
  16 #include <boost/math/special_functions/trunc.hpp>
  17
  18 namespace boost{ namespace math{ namespace detail{
  19
  20 //
  21 // The main method used is due to Hill:
  22 //
  23 // G. W. Hill, Algorithm 396, Student's t-Quantiles,
  24 // Communications of the ACM, 13(10): 619-620, Oct., 1970.
  25 //
  26 template <class T, class Policy>
  27 T inverse_students_t_hill(T ndf, T u, const Policy& pol)
  28 {
  29    BOOST_MATH_STD_USING
  30    BOOST_ASSERT(u <= 0.5);
  31
  32    T a, b, c, d, q, x, y;
  33
  34    if (ndf > 1e20f)
  35       return -boost::math::erfc_inv(2 * u, pol) * constants::root_two<T>();
  36
  37    a = 1 / (ndf - 0.5f);
  38    b = 48 / (a * a);
  39    c = ((20700 * a / b - 98) * a - 16) * a + 96.36f;
  40    d = ((94.5f / (b + c) - 3) / b + 1) * sqrt(a * constants::pi<T>() / 2) * ndf;
  41    y = pow(d * 2 * u, 2 / ndf);
  42
  43    if (y > (0.05f + a))
  44    {
  45       //
  46       // Asymptotic inverse expansion about normal:
  47       //
  48       x = -boost::math::erfc_inv(2 * u, pol) * constants::root_two<T>();
  49       y = x * x;
  50
  51       if (ndf < 5)
  52          c += 0.3f * (ndf - 4.5f) * (x + 0.6f);
  53       c += (((0.05f * d * x - 5) * x - 7) * x - 2) * x + b;
  54       y = (((((0.4f * y + 6.3f) * y + 36) * y + 94.5f) / c - y - 3) / b + 1) * x;
  55       y = boost::math::expm1(a * y * y, pol);
  56    }
  57    else
  58    {
  59       y = ((1 / (((ndf + 6) / (ndf * y) - 0.089f * d - 0.822f)
  60               * (ndf + 2) * 3) + 0.5 / (ndf + 4)) * y - 1)
  61               * (ndf + 1) / (ndf + 2) + 1 / y;
  62    }
  63    q = sqrt(ndf * y);
  64
  65    return -q;
  66 }
  67 //
  68 // Tail and body series are due to Shaw:
  69 //
  70 // www.mth.kcl.ac.uk/~shaww/web_page/papers/Tdistribution06.pdf
  71 //
  72 // Shaw, W.T., 2006, "Sampling Student's T distribution - use of
  73 // the inverse cumulative distribution function."
  74 // Journal of Computational Finance, Vol 9 Issue 4, pp 37-73, Summer 2006
  75 //
  76 template <class T, class Policy>
  77 T inverse_students_t_tail_series(T df, T v, const Policy& pol)
  78 {
  79    BOOST_MATH_STD_USING
  80    // Tail series expansion, see section 6 of Shaw's paper.
  81    // w is calculated using Eq 60:
  82    T w = boost::math::tgamma_delta_ratio(df / 2, constants::half<T>(), pol)
  83       * sqrt(df * constants::pi<T>()) * v;
  84    // define some variables:
  85    T np2 = df + 2;
  86    T np4 = df + 4;
  87    T np6 = df + 6;
  88    //
  89    // Calculate the coefficients d(k), these depend only on the
  90    // number of degrees of freedom df, so at least in theory
  91    // we could tabulate these for fixed df, see p15 of Shaw:
  92    //
  93    T d[7] = { 1, };
  94    d[1] = -(df + 1) / (2 * np2);
  95    np2 *= (df + 2);
  96    d[2] = -df * (df + 1) * (df + 3) / (8 * np2 * np4);
  97    np2 *= df + 2;
  98    d[3] = -df * (df + 1) * (df + 5) * (((3 * df) + 7) * df -2) / (48 * np2 * np4 * np6);
  99    np2 *= (df + 2);
 100    np4 *= (df + 4);
 101    d[4] = -df * (df + 1) * (df + 7) *
 102       ( (((((15 * df) + 154) * df + 465) * df + 286) * df - 336) * df + 64 )
 103       / (384 * np2 * np4 * np6 * (df + 8));
 104    np2 *= (df + 2);
 105    d[5] = -df * (df + 1) * (df + 3) * (df + 9)
 106             * (((((((35 * df + 452) * df + 1573) * df + 600) * df - 2020) * df) + 928) * df -128)
 107             / (1280 * np2 * np4 * np6 * (df + 8) * (df + 10));
 108    np2 *= (df + 2);
 109    np4 *= (df + 4);
 110    np6 *= (df + 6);
 111    d[6] = -df * (df + 1) * (df + 11)
 112             * ((((((((((((945 * df) + 31506) * df + 425858) * df + 2980236) * df + 11266745) * df + 20675018) * df + 7747124) * df - 22574632) * df - 8565600) * df + 18108416) * df - 7099392) * df + 884736)
 113             / (46080 * np2 * np4 * np6 * (df + 8) * (df + 10) * (df +12));
 114    //
 115    // Now bring everthing together to provide the result,
 116    // this is Eq 62 of Shaw:
 117    //
 118    T rn = sqrt(df);
 119    T div = pow(rn * w, 1 / df);
 120    T power = div * div;
 121    T result = tools::evaluate_polynomial<7, T, T>(d, power);
 122    result *= rn;
 123    result /= div;
 124    return -result;
 125 }
 126
 127 template <class T, class Policy>
 128 T inverse_students_t_body_series(T df, T u, const Policy& pol)
 129 {
 130    BOOST_MATH_STD_USING
 131    //
 132    // Body series for small N:
 133    //
 134    // Start with Eq 56 of Shaw:
 135    //
 136    T v = boost::math::tgamma_delta_ratio(df / 2, constants::half<T>(), pol)
 137       * sqrt(df * constants::pi<T>()) * (u - constants::half<T>());
 138    //
 139    // Workspace for the polynomial coefficients:
 140    //
 141    T c[11] = { 0, 1, };
 142    //
 143    // Figure out what the coefficients are, note these depend
 144    // only on the degrees of freedom (Eq 57 of Shaw):
 145    //
 146    T in = 1 / df;
 147    c[2] = 0.16666666666666666667 + 0.16666666666666666667 * in;
 148    c[3] = (0.0083333333333333333333 * in
 149       + 0.066666666666666666667) * in
 150       + 0.058333333333333333333;
 151    c[4] = ((0.00019841269841269841270 * in
 152       + 0.0017857142857142857143) * in
 153       + 0.026785714285714285714) * in
 154       + 0.025198412698412698413;
 155    c[5] = (((2.7557319223985890653e-6 * in
 156       + 0.00037477954144620811287) * in
 157       - 0.0011078042328042328042) * in
 158       + 0.010559964726631393298) * in
 159       + 0.012039792768959435626;
 160    c[6] = ((((2.5052108385441718775e-8 * in
 161       - 0.000062705427288760622094) * in
 162       + 0.00059458674042007375341) * in
 163       - 0.0016095979637646304313) * in
 164       + 0.0061039211560044893378) * in
 165       + 0.0038370059724226390893;
 166    c[7] = (((((1.6059043836821614599e-10 * in
 167       + 0.000015401265401265401265) * in
 168       - 0.00016376804137220803887) * in
 169       + 0.00069084207973096861986) * in
 170       - 0.0012579159844784844785) * in
 171       + 0.0010898206731540064873) * in
 172       + 0.0032177478835464946576;
 173    c[8] = ((((((7.6471637318198164759e-13 * in
 174       - 3.9851014346715404916e-6) * in
 175       + 0.000049255746366361445727) * in
 176       - 0.00024947258047043099953) * in
 177       + 0.00064513046951456342991) * in
 178       - 0.00076245135440323932387) * in
 179       + 0.000033530976880017885309) * in
 180       + 0.0017438262298340009980;
 181    c[9] = (((((((2.8114572543455207632e-15 * in
 182       + 1.0914179173496789432e-6) * in
 183       - 0.000015303004486655377567) * in
 184       + 0.000090867107935219902229) * in
 185       - 0.00029133414466938067350) * in
 186       + 0.00051406605788341121363) * in
 187       - 0.00036307660358786885787) * in
 188       - 0.00031101086326318780412) * in
 189       + 0.00096472747321388644237;
 190    c[10] = ((((((((8.2206352466243297170e-18 * in
 191       - 3.1239569599829868045e-7) * in
 192       + 4.8903045291975346210e-6) * in
 193       - 0.000033202652391372058698) * in
 194       + 0.00012645437628698076975) * in
 195       - 0.00028690924218514613987) * in
 196       + 0.00035764655430568632777) * in
 197       - 0.00010230378073700412687) * in
 198       - 0.00036942667800009661203) * in
 199       + 0.00054229262813129686486;
 200    //
 201    // The result is then a polynomial in v (see Eq 56 of Shaw):
 202    //
 203    return tools::evaluate_odd_polynomial<11, T, T>(c, v);
 204 }
 205
 206 template <class T, class Policy>
 207 T inverse_students_t(T df, T u, T v, const Policy& pol, bool* pexact = 0)
 208 {
 209    //
 210    // df = number of degrees of freedom.
 211    // u = probablity.
 212    // v = 1 - u.
 213    // l = lanczos type to use.
 214    //
 215    BOOST_MATH_STD_USING
 216    bool invert = false;
 217    T result = 0;
 218    if(pexact)
 219       *pexact = false;
 220    if(u > v)
 221    {
 222       // function is symmetric, invert it:
 223       std::swap(u, v);
 224       invert = true;
 225    }
 226    if((floor(df) == df) && (df < 20))
 227    {
 228       //
 229       // we have integer degrees of freedom, try for the special
 230       // cases first:
 231       //
 232       T tolerance = ldexp(1.0f, (2 * policies::digits<T, Policy>()) / 3);
 233
 234       switch(itrunc(df, Policy()))
 235       {
 236       case 1:
 237          {
 238             //
 239             // df = 1 is the same as the Cauchy distribution, see
 240             // Shaw Eq 35:
 241             //
 242             if(u == 0.5)
 243                result = 0;
 244             else
 245                result = -cos(constants::pi<T>() * u) / sin(constants::pi<T>() * u);
 246             if(pexact)
 247                *pexact = true;
 248             break;
 249          }
 250       case 2:
 251          {
 252             //
 253             // df = 2 has an exact result, see Shaw Eq 36:
 254             //
 255             result =(2 * u - 1) / sqrt(2 * u * v);
 256             if(pexact)
 257                *pexact = true;
 258             break;
 259          }
 260       case 4:
 261          {
 262             //
 263             // df = 4 has an exact result, see Shaw Eq 38 & 39:
 264             //
 265             T alpha = 4 * u * v;
 266             T root_alpha = sqrt(alpha);
 267             T r = 4 * cos(acos(root_alpha) / 3) / root_alpha;
 268             T x = sqrt(r - 4);
 269             result = u - 0.5f < 0 ? (T)-x : x;
 270             if(pexact)
 271                *pexact = true;
 272             break;
 273          }
 274       case 6:
 275          {
 276             //
 277             // We get numeric overflow in this area:
 278             //
 279             if(u < 1e-150)
 280                return (invert ? -1 : 1) * inverse_students_t_hill(df, u, pol);
 281             //
 282             // Newton-Raphson iteration of a polynomial case,
 283             // choice of seed value is taken from Shaw's online
 284             // supplement:
 285             //
 286             T a = 4 * (u - u * u);//1 - 4 * (u - 0.5f) * (u - 0.5f);
 287             T b = boost::math::cbrt(a);
 288             static const T c = 0.85498797333834849467655443627193;
 289             T p = 6 * (1 + c * (1 / b - 1));
 290             T p0;
 291             do{
 292                T p2 = p * p;
 293                T p4 = p2 * p2;
 294                T p5 = p * p4;
 295                p0 = p;
 296                // next term is given by Eq 41:
 297                p = 2 * (8 * a * p5 - 270 * p2 + 2187) / (5 * (4 * a * p4 - 216 * p - 243));
 298             }while(fabs((p - p0) / p) > tolerance);
 299             //
 300             // Use Eq 45 to extract the result:
 301             //
 302             p = sqrt(p - df);
 303             result = (u - 0.5f) < 0 ? (T)-p : p;
 304             break;
 305          }
 306 #if 0
 307          //
 308          // These are Shaw's "exact" but iterative solutions
 309          // for even df, the numerical accuracy of these is
 310          // rather less than Hill's method, so these are disabled
 311          // for now, which is a shame because they are reasonably
 312          // quick to evaluate...
 313          //
 314       case 8:
 315          {
 316             //
 317             // Newton-Raphson iteration of a polynomial case,
 318             // choice of seed value is taken from Shaw's online
 319             // supplement:
 320             //
 321             static const T c8 = 0.85994765706259820318168359251872L;
 322             T a = 4 * (u - u * u); //1 - 4 * (u - 0.5f) * (u - 0.5f);
 323             T b = pow(a, T(1) / 4);
 324             T p = 8 * (1 + c8 * (1 / b - 1));
 325             T p0 = p;
 326             do{
 327                T p5 = p * p;
 328                p5 *= p5 * p;
 329                p0 = p;
 330                // Next term is given by Eq 42:
 331                p = 2 * (3 * p + (640 * (160 + p * (24 + p * (p + 4)))) / (-5120 + p * (-2048 - 960 * p + a * p5))) / 7;
 332             }while(fabs((p - p0) / p) > tolerance);
 333             //
 334             // Use Eq 45 to extract the result:
 335             //
 336             p = sqrt(p - df);
 337             result = (u - 0.5f) < 0 ? -p : p;
 338             break;
 339          }
 340       case 10:
 341          {
 342             //
 343             // Newton-Raphson iteration of a polynomial case,
 344             // choice of seed value is taken from Shaw's online
 345             // supplement:
 346             //
 347             static const T c10 = 0.86781292867813396759105692122285L;
 348             T a = 4 * (u - u * u); //1 - 4 * (u - 0.5f) * (u - 0.5f);
 349             T b = pow(a, T(1) / 5);
 350             T p = 10 * (1 + c10 * (1 / b - 1));
 351             T p0;
 352             do{
 353                T p6 = p * p;
 354                p6 *= p6 * p6;
 355                p0 = p;
 356                // Next term given by Eq 43:
 357                p = (8 * p) / 9 + (218750 * (21875 + 4 * p * (625 + p * (75 + 2 * p * (5 + p))))) /
 358                   (9 * (-68359375 + 8 * p * (-2343750 + p * (-546875 - 175000 * p + 8 * a * p6))));
 359             }while(fabs((p - p0) / p) > tolerance);
 360             //
 361             // Use Eq 45 to extract the result:
 362             //
 363             p = sqrt(p - df);
 364             result = (u - 0.5f) < 0 ? -p : p;
 365             break;
 366          }
 367 #endif
 368       default:
 369          goto calculate_real;
 370       }
 371    }
 372    else
 373    {
 374 calculate_real:
 375       if(df > 0x10000000)
 376       {
 377          result = -boost::math::erfc_inv(2 * u, pol) * constants::root_two<T>();
 378          if((pexact) && (df >= 1e20))
 379             *pexact = true;
 380       }
 381       else if(df < 3)
 382       {
 383          //
 384          // Use a roughly linear scheme to choose between Shaw's
 385          // tail series and body series:
 386          //
 387          T crossover = 0.2742f - df * 0.0242143f;
 388          if(u > crossover)
 389          {
 390             result = boost::math::detail::inverse_students_t_body_series(df, u, pol);
 391          }
 392          else
 393          {
 394             result = boost::math::detail::inverse_students_t_tail_series(df, u, pol);
 395          }
 396       }
 397       else
 398       {
 399          //
 400          // Use Hill's method except in the exteme tails
 401          // where we use Shaw's tail series.
 402          // The crossover point is roughly exponential in -df:
 403          //
 404          T crossover = ldexp(1.0f, iround(T(df / -0.654f), typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type()));
 405          if(u > crossover)
 406          {
 407             result = boost::math::detail::inverse_students_t_hill(df, u, pol);
 408          }
 409          else
 410          {
 411             result = boost::math::detail::inverse_students_t_tail_series(df, u, pol);
 412          }
 413       }
 414    }
 415    return invert ? (T)-result : result;
 416 }
 417
 418 template <class T, class Policy>
 419 inline T find_ibeta_inv_from_t_dist(T a, T p, T /*q*/, T* py, const Policy& pol)
 420 {
 421    T u = p / 2;
 422    T v = 1 - u;
 423    T df = a * 2;
 424    T t = boost::math::detail::inverse_students_t(df, u, v, pol);
 425    *py = t * t / (df + t * t);
 426    return df / (df + t * t);
 427 }
 428
 429 template <class T, class Policy>
 430 inline T fast_students_t_quantile_imp(T df, T p, const Policy& pol, const mpl::false_*)
 431 {
 432    BOOST_MATH_STD_USING
 433    //
 434    // Need to use inverse incomplete beta to get
 435    // required precision so not so fast:
 436    //
 437    T probability = (p > 0.5) ? 1 - p : p;
 438    T t, x, y(0);
 439    x = ibeta_inv(df / 2, T(0.5), 2 * probability, &y, pol);
 440    if(df * y > tools::max_value<T>() * x)
 441       t = policies::raise_overflow_error<T>("boost::math::students_t_quantile<%1%>(%1%,%1%)", 0, pol);
 442    else
 443       t = sqrt(df * y / x);
 444    //
 445    // Figure out sign based on the size of p:
 446    //
 447    if(p < 0.5)
 448       t = -t;
 449    return t;
 450 }
 451
 452 template <class T, class Policy>
 453 T fast_students_t_quantile_imp(T df, T p, const Policy& pol, const mpl::true_*)
 454 {
 455    BOOST_MATH_STD_USING
 456    bool invert = false;
 457    if((df < 2) && (floor(df) != df))
 458       return boost::math::detail::fast_students_t_quantile_imp(df, p, pol, static_cast<mpl::false_*>(0));
 459    if(p > 0.5)
 460    {
 461       p = 1 - p;
 462       invert = true;
 463    }
 464    //
 465    // Get an estimate of the result:
 466    //
 467    bool exact;
 468    T t = inverse_students_t(df, p, T(1-p), pol, &exact);
 469    if((t == 0) || exact)
 470       return invert ? -t : t; // can't do better!
 471    //
 472    // Change variables to inverse incomplete beta:
 473    //
 474    T t2 = t * t;
 475    T xb = df / (df + t2);
 476    T y = t2 / (df + t2);
 477    T a = df / 2;
 478    //
 479    // t can be so large that x underflows,
 480    // just return our estimate in that case:
 481    //
 482    if(xb == 0)
 483       return t;
 484    //
 485    // Get incomplete beta and it's derivative:
 486    //
 487    T f1;
 488    T f0 = xb < y ? ibeta_imp(a, constants::half<T>(), xb, pol, false, true, &f1)
 489       : ibeta_imp(constants::half<T>(), a, y, pol, true, true, &f1);
 490
 491    // Get cdf from incomplete beta result:
 492    T p0 = f0 / 2  - p;
 493    // Get pdf from derivative:
 494    T p1 = f1 * sqrt(y * xb * xb * xb / df);
 495    //
 496    // Second derivative divided by p1:
 497    //
 498    // yacas gives:
 499    //
 500    // In> PrettyForm(Simplify(D(t) (1 + t^2/v) ^ (-(v+1)/2)))
 501    //
 502    //  |                        | v + 1     |     |
 503    //  |                       -| ----- + 1 |     |
 504    //  |                        |   2       |     |
 505    // -|             |  2     |                   |
 506    //  |             | t      |                   |
 507    //  |             | -- + 1 |                   |
 508    //  | ( v + 1 ) * | v      |               * t |
 509    // ---------------------------------------------
 510    //                       v
 511    //
 512    // Which after some manipulation is:
 513    //
 514    // -p1 * t * (df + 1) / (t^2 + df)
 515    //
 516    T p2 = t * (df + 1) / (t * t + df);
 517    // Halley step:
 518    t = fabs(t);
 519    t += p0 / (p1 + p0 * p2 / 2);
 520    return !invert ? -t : t;
 521 }
 522
 523 template <class T, class Policy>
 524 inline T fast_students_t_quantile(T df, T p, const Policy& pol)
 525 {
 526    typedef typename policies::evaluation<T, Policy>::type value_type;
 527    typedef typename policies::normalise<
 528       Policy,
 529       policies::promote_float<false>,
 530       policies::promote_double<false>,
 531       policies::discrete_quantile<>,
 532       policies::assert_undefined<> >::type forwarding_policy;
 533
 534    typedef mpl::bool_<
 535       (std::numeric_limits<T>::digits <= 53)
 536        &&
 537       (std::numeric_limits<T>::is_specialized)
 538        &&
 539       (std::numeric_limits<T>::radix == 2)
 540    > tag_type;
 541    return policies::checked_narrowing_cast<T, forwarding_policy>(fast_students_t_quantile_imp(static_cast<value_type>(df), static_cast<value_type>(p), pol, static_cast<tag_type*>(0)), "boost::math::students_t_quantile<%1%>(%1%,%1%,%1%)");
 542 }
 543
 544 }}} // namespaces
 545
 546 #endif // BOOST_MATH_SF_DETAIL_INV_T_HPP
 547
 548
 549