special_functions/detail/t_distribution_inv.hpp

0001 //  Copyright John Maddock 2007.
0002 //  Copyright Paul A. Bristow 2007
0003 //  Use, modification and distribution are subject to the
0004 //  Boost Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_MATH_SF_DETAIL_INV_T_HPP
0008 #define BOOST_MATH_SF_DETAIL_INV_T_HPP
0009
0010 #ifdef _MSC_VER
0011 #pragma once
0012 #endif
0013
0014 #include <boost/math/special_functions/cbrt.hpp>
0015 #include <boost/math/special_functions/round.hpp>
0016 #include <boost/math/special_functions/trunc.hpp>
0017
0018 namespace boost{ namespace math{ namespace detail{
0019
0020 //
0021 // The main method used is due to Hill:
0022 //
0023 // G. W. Hill, Algorithm 396, Student's t-Quantiles,
0024 // Communications of the ACM, 13(10): 619-620, Oct., 1970.
0025 //
0026 template <class T, class Policy>
0027 T inverse_students_t_hill(T ndf, T u, const Policy& pol)
0028 {
0029    BOOST_MATH_STD_USING
0030    BOOST_MATH_ASSERT(u <= 0.5);
0031
0032    T a, b, c, d, q, x, y;
0033
0034    if (ndf > 1e20f)
0035       return -boost::math::erfc_inv(2 * u, pol) * constants::root_two<T>();
0036
0037    a = 1 / (ndf - 0.5f);
0038    b = 48 / (a * a);
0039    c = ((20700 * a / b - 98) * a - 16) * a + 96.36f;
0040    d = ((94.5f / (b + c) - 3) / b + 1) * sqrt(a * constants::pi<T>() / 2) * ndf;
0041    y = pow(d * 2 * u, 2 / ndf);
0042
0043    if (y > (0.05f + a))
0044    {
0045       //
0046       // Asymptotic inverse expansion about normal:
0047       //
0048       x = -boost::math::erfc_inv(2 * u, pol) * constants::root_two<T>();
0049       y = x * x;
0050
0051       if (ndf < 5)
0052          c += 0.3f * (ndf - 4.5f) * (x + 0.6f);
0053       c += (((0.05f * d * x - 5) * x - 7) * x - 2) * x + b;
0054       y = (((((0.4f * y + 6.3f) * y + 36) * y + 94.5f) / c - y - 3) / b + 1) * x;
0055       y = boost::math::expm1(a * y * y, pol);
0056    }
0057    else
0058    {
0059       y = static_cast<T>(((1 / (((ndf + 6) / (ndf * y) - 0.089f * d - 0.822f)
0060               * (ndf + 2) * 3) + 0.5 / (ndf + 4)) * y - 1)
0061               * (ndf + 1) / (ndf + 2) + 1 / y);
0062    }
0063    q = sqrt(ndf * y);
0064
0065    return -q;
0066 }
0067 //
0068 // Tail and body series are due to Shaw:
0069 //
0070 // www.mth.kcl.ac.uk/~shaww/web_page/papers/Tdistribution06.pdf
0071 //
0072 // Shaw, W.T., 2006, "Sampling Student's T distribution - use of
0073 // the inverse cumulative distribution function."
0074 // Journal of Computational Finance, Vol 9 Issue 4, pp 37-73, Summer 2006
0075 //
0076 template <class T, class Policy>
0077 T inverse_students_t_tail_series(T df, T v, const Policy& pol)
0078 {
0079    BOOST_MATH_STD_USING
0080    // Tail series expansion, see section 6 of Shaw's paper.
0081    // w is calculated using Eq 60:
0082    T w = boost::math::tgamma_delta_ratio(df / 2, constants::half<T>(), pol)
0083       * sqrt(df * constants::pi<T>()) * v;
0084    // define some variables:
0085    T np2 = df + 2;
0086    T np4 = df + 4;
0087    T np6 = df + 6;
0088    //
0089    // Calculate the coefficients d(k), these depend only on the
0090    // number of degrees of freedom df, so at least in theory
0091    // we could tabulate these for fixed df, see p15 of Shaw:
0092    //
0093    T d[7] = { 1, };
0094    d[1] = -(df + 1) / (2 * np2);
0095    np2 *= (df + 2);
0096    d[2] = -df * (df + 1) * (df + 3) / (8 * np2 * np4);
0097    np2 *= df + 2;
0098    d[3] = -df * (df + 1) * (df + 5) * (((3 * df) + 7) * df -2) / (48 * np2 * np4 * np6);
0099    np2 *= (df + 2);
0100    np4 *= (df + 4);
0101    d[4] = -df * (df + 1) * (df + 7) *
0102       ( (((((15 * df) + 154) * df + 465) * df + 286) * df - 336) * df + 64 )
0103       / (384 * np2 * np4 * np6 * (df + 8));
0104    np2 *= (df + 2);
0105    d[5] = -df * (df + 1) * (df + 3) * (df + 9)
0106             * (((((((35 * df + 452) * df + 1573) * df + 600) * df - 2020) * df) + 928) * df -128)
0107             / (1280 * np2 * np4 * np6 * (df + 8) * (df + 10));
0108    np2 *= (df + 2);
0109    np4 *= (df + 4);
0110    np6 *= (df + 6);
0111    d[6] = -df * (df + 1) * (df + 11)
0112             * ((((((((((((945 * df) + 31506) * df + 425858) * df + 2980236) * df + 11266745) * df + 20675018) * df + 7747124) * df - 22574632) * df - 8565600) * df + 18108416) * df - 7099392) * df + 884736)
0113             / (46080 * np2 * np4 * np6 * (df + 8) * (df + 10) * (df +12));
0114    //
0115    // Now bring everything together to provide the result,
0116    // this is Eq 62 of Shaw:
0117    //
0118    T rn = sqrt(df);
0119    T div = pow(rn * w, 1 / df);
0120    T power = div * div;
0121    T result = tools::evaluate_polynomial<7, T, T>(d, power);
0122    result *= rn;
0123    result /= div;
0124    return -result;
0125 }
0126
0127 template <class T, class Policy>
0128 T inverse_students_t_body_series(T df, T u, const Policy& pol)
0129 {
0130    BOOST_MATH_STD_USING
0131    //
0132    // Body series for small N:
0133    //
0134    // Start with Eq 56 of Shaw:
0135    //
0136    T v = boost::math::tgamma_delta_ratio(df / 2, constants::half<T>(), pol)
0137       * sqrt(df * constants::pi<T>()) * (u - constants::half<T>());
0138    //
0139    // Workspace for the polynomial coefficients:
0140    //
0141    T c[11] = { 0, 1, };
0142    //
0143    // Figure out what the coefficients are, note these depend
0144    // only on the degrees of freedom (Eq 57 of Shaw):
0145    //
0146    T in = 1 / df;
0147    c[2] = static_cast<T>(0.16666666666666666667 + 0.16666666666666666667 * in);
0148    c[3] = static_cast<T>((0.0083333333333333333333 * in
0149       + 0.066666666666666666667) * in
0150       + 0.058333333333333333333);
0151    c[4] = static_cast<T>(((0.00019841269841269841270 * in
0152       + 0.0017857142857142857143) * in
0153       + 0.026785714285714285714) * in
0154       + 0.025198412698412698413);
0155    c[5] = static_cast<T>((((2.7557319223985890653e-6 * in
0156       + 0.00037477954144620811287) * in
0157       - 0.0011078042328042328042) * in
0158       + 0.010559964726631393298) * in
0159       + 0.012039792768959435626);
0160    c[6] = static_cast<T>(((((2.5052108385441718775e-8 * in
0161       - 0.000062705427288760622094) * in
0162       + 0.00059458674042007375341) * in
0163       - 0.0016095979637646304313) * in
0164       + 0.0061039211560044893378) * in
0165       + 0.0038370059724226390893);
0166    c[7] = static_cast<T>((((((1.6059043836821614599e-10 * in
0167       + 0.000015401265401265401265) * in
0168       - 0.00016376804137220803887) * in
0169       + 0.00069084207973096861986) * in
0170       - 0.0012579159844784844785) * in
0171       + 0.0010898206731540064873) * in
0172       + 0.0032177478835464946576);
0173    c[8] = static_cast<T>(((((((7.6471637318198164759e-13 * in
0174       - 3.9851014346715404916e-6) * in
0175       + 0.000049255746366361445727) * in
0176       - 0.00024947258047043099953) * in
0177       + 0.00064513046951456342991) * in
0178       - 0.00076245135440323932387) * in
0179       + 0.000033530976880017885309) * in
0180       + 0.0017438262298340009980);
0181    c[9] = static_cast<T>((((((((2.8114572543455207632e-15 * in
0182       + 1.0914179173496789432e-6) * in
0183       - 0.000015303004486655377567) * in
0184       + 0.000090867107935219902229) * in
0185       - 0.00029133414466938067350) * in
0186       + 0.00051406605788341121363) * in
0187       - 0.00036307660358786885787) * in
0188       - 0.00031101086326318780412) * in
0189       + 0.00096472747321388644237);
0190    c[10] = static_cast<T>(((((((((8.2206352466243297170e-18 * in
0191       - 3.1239569599829868045e-7) * in
0192       + 4.8903045291975346210e-6) * in
0193       - 0.000033202652391372058698) * in
0194       + 0.00012645437628698076975) * in
0195       - 0.00028690924218514613987) * in
0196       + 0.00035764655430568632777) * in
0197       - 0.00010230378073700412687) * in
0198       - 0.00036942667800009661203) * in
0199       + 0.00054229262813129686486);
0200    //
0201    // The result is then a polynomial in v (see Eq 56 of Shaw):
0202    //
0203    return tools::evaluate_odd_polynomial<11, T, T>(c, v);
0204 }
0205
0206 template <class T, class Policy>
0207 T inverse_students_t(T df, T u, T v, const Policy& pol, bool* pexact = nullptr)
0208 {
0209    //
0210    // df = number of degrees of freedom.
0211    // u = probability.
0212    // v = 1 - u.
0213    // l = lanczos type to use.
0214    //
0215    BOOST_MATH_STD_USING
0216    bool invert = false;
0217    T result = 0;
0218    if(pexact)
0219       *pexact = false;
0220    if(u > v)
0221    {
0222       // function is symmetric, invert it:
0223       std::swap(u, v);
0224       invert = true;
0225    }
0226    if((floor(df) == df) && (df < 20))
0227    {
0228       //
0229       // we have integer degrees of freedom, try for the special
0230       // cases first:
0231       //
0232       T tolerance = ldexp(1.0f, (2 * policies::digits<T, Policy>()) / 3);
0233
0234       switch(itrunc(df, Policy()))
0235       {
0236       case 1:
0237          {
0238             //
0239             // df = 1 is the same as the Cauchy distribution, see
0240             // Shaw Eq 35:
0241             //
0242             if(u == 0.5)
0243                result = 0;
0244             else
0245                result = -cos(constants::pi<T>() * u) / sin(constants::pi<T>() * u);
0246             if(pexact)
0247                *pexact = true;
0248             break;
0249          }
0250       case 2:
0251          {
0252             //
0253             // df = 2 has an exact result, see Shaw Eq 36:
0254             //
0255             result =(2 * u - 1) / sqrt(2 * u * v);
0256             if(pexact)
0257                *pexact = true;
0258             break;
0259          }
0260       case 4:
0261          {
0262             //
0263             // df = 4 has an exact result, see Shaw Eq 38 & 39:
0264             //
0265             T alpha = 4 * u * v;
0266             T root_alpha = sqrt(alpha);
0267             T r = 4 * cos(acos(root_alpha) / 3) / root_alpha;
0268             T x = sqrt(r - 4);
0269             result = u - 0.5f < 0 ? (T)-x : x;
0270             if(pexact)
0271                *pexact = true;
0272             break;
0273          }
0274       case 6:
0275          {
0276             //
0277             // We get numeric overflow in this area:
0278             //
0279             if(u < 1e-150)
0280                return (invert ? -1 : 1) * inverse_students_t_hill(df, u, pol);
0281             //
0282             // Newton-Raphson iteration of a polynomial case,
0283             // choice of seed value is taken from Shaw's online
0284             // supplement:
0285             //
0286             T a = 4 * (u - u * u);//1 - 4 * (u - 0.5f) * (u - 0.5f);
0287             T b = boost::math::cbrt(a, pol);
0288             static const T c = static_cast<T>(0.85498797333834849467655443627193);
0289             T p = 6 * (1 + c * (1 / b - 1));
0290             T p0;
0291             do{
0292                T p2 = p * p;
0293                T p4 = p2 * p2;
0294                T p5 = p * p4;
0295                p0 = p;
0296                // next term is given by Eq 41:
0297                p = 2 * (8 * a * p5 - 270 * p2 + 2187) / (5 * (4 * a * p4 - 216 * p - 243));
0298             }while(fabs((p - p0) / p) > tolerance);
0299             //
0300             // Use Eq 45 to extract the result:
0301             //
0302             p = sqrt(p - df);
0303             result = (u - 0.5f) < 0 ? (T)-p : p;
0304             break;
0305          }
0306 #if 0
0307          //
0308          // These are Shaw's "exact" but iterative solutions
0309          // for even df, the numerical accuracy of these is
0310          // rather less than Hill's method, so these are disabled
0311          // for now, which is a shame because they are reasonably
0312          // quick to evaluate...
0313          //
0314       case 8:
0315          {
0316             //
0317             // Newton-Raphson iteration of a polynomial case,
0318             // choice of seed value is taken from Shaw's online
0319             // supplement:
0320             //
0321             static const T c8 = 0.85994765706259820318168359251872L;
0322             T a = 4 * (u - u * u); //1 - 4 * (u - 0.5f) * (u - 0.5f);
0323             T b = pow(a, T(1) / 4);
0324             T p = 8 * (1 + c8 * (1 / b - 1));
0325             T p0 = p;
0326             do{
0327                T p5 = p * p;
0328                p5 *= p5 * p;
0329                p0 = p;
0330                // Next term is given by Eq 42:
0331                p = 2 * (3 * p + (640 * (160 + p * (24 + p * (p + 4)))) / (-5120 + p * (-2048 - 960 * p + a * p5))) / 7;
0332             }while(fabs((p - p0) / p) > tolerance);
0333             //
0334             // Use Eq 45 to extract the result:
0335             //
0336             p = sqrt(p - df);
0337             result = (u - 0.5f) < 0 ? -p : p;
0338             break;
0339          }
0340       case 10:
0341          {
0342             //
0343             // Newton-Raphson iteration of a polynomial case,
0344             // choice of seed value is taken from Shaw's online
0345             // supplement:
0346             //
0347             static const T c10 = 0.86781292867813396759105692122285L;
0348             T a = 4 * (u - u * u); //1 - 4 * (u - 0.5f) * (u - 0.5f);
0349             T b = pow(a, T(1) / 5);
0350             T p = 10 * (1 + c10 * (1 / b - 1));
0351             T p0;
0352             do{
0353                T p6 = p * p;
0354                p6 *= p6 * p6;
0355                p0 = p;
0356                // Next term given by Eq 43:
0357                p = (8 * p) / 9 + (218750 * (21875 + 4 * p * (625 + p * (75 + 2 * p * (5 + p))))) /
0358                   (9 * (-68359375 + 8 * p * (-2343750 + p * (-546875 - 175000 * p + 8 * a * p6))));
0359             }while(fabs((p - p0) / p) > tolerance);
0360             //
0361             // Use Eq 45 to extract the result:
0362             //
0363             p = sqrt(p - df);
0364             result = (u - 0.5f) < 0 ? -p : p;
0365             break;
0366          }
0367 #endif
0368       default:
0369          goto calculate_real;
0370       }
0371    }
0372    else
0373    {
0374 calculate_real:
0375       if(df > 0x10000000)
0376       {
0377          result = -boost::math::erfc_inv(2 * u, pol) * constants::root_two<T>();
0378          if((pexact) && (df >= 1e20))
0379             *pexact = true;
0380       }
0381       else if(df < 3)
0382       {
0383          //
0384          // Use a roughly linear scheme to choose between Shaw's
0385          // tail series and body series:
0386          //
0387          T crossover = 0.2742f - df * 0.0242143f;
0388          if(u > crossover)
0389          {
0390             result = boost::math::detail::inverse_students_t_body_series(df, u, pol);
0391          }
0392          else
0393          {
0394             result = boost::math::detail::inverse_students_t_tail_series(df, u, pol);
0395          }
0396       }
0397       else
0398       {
0399          //
0400          // Use Hill's method except in the extreme tails
0401          // where we use Shaw's tail series.
0402          // The crossover point is roughly exponential in -df:
0403          //
0404          T crossover = ldexp(1.0f, iround(T(df / -0.654f), typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type()));
0405          if(u > crossover)
0406          {
0407             result = boost::math::detail::inverse_students_t_hill(df, u, pol);
0408          }
0409          else
0410          {
0411             result = boost::math::detail::inverse_students_t_tail_series(df, u, pol);
0412          }
0413       }
0414    }
0415    return invert ? (T)-result : result;
0416 }
0417
0418 template <class T, class Policy>
0419 inline T find_ibeta_inv_from_t_dist(T a, T p, T /*q*/, T* py, const Policy& pol)
0420 {
0421    T u = p / 2;
0422    T v = 1 - u;
0423    T df = a * 2;
0424    T t = boost::math::detail::inverse_students_t(df, u, v, pol);
0425    *py = t * t / (df + t * t);
0426    return df / (df + t * t);
0427 }
0428
0429 template <class T, class Policy>
0430 inline T fast_students_t_quantile_imp(T df, T p, const Policy& pol, const std::false_type*)
0431 {
0432    BOOST_MATH_STD_USING
0433    //
0434    // Need to use inverse incomplete beta to get
0435    // required precision so not so fast:
0436    //
0437    T probability = (p > 0.5) ? 1 - p : p;
0438    T t, x, y(0);
0439    x = ibeta_inv(df / 2, T(0.5), 2 * probability, &y, pol);
0440    if(df * y > tools::max_value<T>() * x)
0441       t = policies::raise_overflow_error<T>("boost::math::students_t_quantile<%1%>(%1%,%1%)", nullptr, pol);
0442    else
0443       t = sqrt(df * y / x);
0444    //
0445    // Figure out sign based on the size of p:
0446    //
0447    if(p < 0.5)
0448       t = -t;
0449    return t;
0450 }
0451
0452 template <class T, class Policy>
0453 T fast_students_t_quantile_imp(T df, T p, const Policy& pol, const std::true_type*)
0454 {
0455    BOOST_MATH_STD_USING
0456    bool invert = false;
0457    if((df < 2) && (floor(df) != df))
0458       return boost::math::detail::fast_students_t_quantile_imp(df, p, pol, static_cast<std::false_type*>(nullptr));
0459    if(p > 0.5)
0460    {
0461       p = 1 - p;
0462       invert = true;
0463    }
0464    //
0465    // Get an estimate of the result:
0466    //
0467    bool exact;
0468    T t = inverse_students_t(df, p, T(1-p), pol, &exact);
0469    if((t == 0) || exact)
0470       return invert ? -t : t; // can't do better!
0471    //
0472    // Change variables to inverse incomplete beta:
0473    //
0474    T t2 = t * t;
0475    T xb = df / (df + t2);
0476    T y = t2 / (df + t2);
0477    T a = df / 2;
0478    //
0479    // t can be so large that x underflows,
0480    // just return our estimate in that case:
0481    //
0482    if(xb == 0)
0483       return t;
0484    //
0485    // Get incomplete beta and it's derivative:
0486    //
0487    T f1;
0488    T f0 = xb < y ? ibeta_imp(a, constants::half<T>(), xb, pol, false, true, &f1)
0489       : ibeta_imp(constants::half<T>(), a, y, pol, true, true, &f1);
0490
0491    // Get cdf from incomplete beta result:
0492    T p0 = f0 / 2  - p;
0493    // Get pdf from derivative:
0494    T p1 = f1 * sqrt(y * xb * xb * xb / df);
0495    //
0496    // Second derivative divided by p1:
0497    //
0498    // yacas gives:
0499    //
0500    // In> PrettyForm(Simplify(D(t) (1 + t^2/v) ^ (-(v+1)/2)))
0501    //
0502    //  |                        | v + 1     |     |
0503    //  |                       -| ----- + 1 |     |
0504    //  |                        |   2       |     |
0505    // -|             |  2     |                   |
0506    //  |             | t      |                   |
0507    //  |             | -- + 1 |                   |
0508    //  | ( v + 1 ) * | v      |               * t |
0509    // ---------------------------------------------
0510    //                       v
0511    //
0512    // Which after some manipulation is:
0513    //
0514    // -p1 * t * (df + 1) / (t^2 + df)
0515    //
0516    T p2 = t * (df + 1) / (t * t + df);
0517    // Halley step:
0518    t = fabs(t);
0519    t += p0 / (p1 + p0 * p2 / 2);
0520    return !invert ? -t : t;
0521 }
0522
0523 template <class T, class Policy>
0524 inline T fast_students_t_quantile(T df, T p, const Policy& pol)
0525 {
0526    typedef typename policies::evaluation<T, Policy>::type value_type;
0527    typedef typename policies::normalise<
0528       Policy,
0529       policies::promote_float<false>,
0530       policies::promote_double<false>,
0531       policies::discrete_quantile<>,
0532       policies::assert_undefined<> >::type forwarding_policy;
0533
0534    typedef std::integral_constant<bool,
0535       (std::numeric_limits<T>::digits <= 53)
0536        &&
0537       (std::numeric_limits<T>::is_specialized)
0538        &&
0539       (std::numeric_limits<T>::radix == 2)
0540    > tag_type;
0541    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*>(nullptr)), "boost::math::students_t_quantile<%1%>(%1%,%1%,%1%)");
0542 }
0543
0544 }}} // namespaces
0545
0546 #endif // BOOST_MATH_SF_DETAIL_INV_T_HPP
0547
0548
0549