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