t_distribution_inv.hpp

来自「Boost provides free peer-reviewed portab」· HPP 代码 · 共 517 行 · 第 1/2 页

//  Copyright John Maddock 2007.
//  Copyright Paul A. Bristow 2007
//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_MATH_SF_DETAIL_INV_T_HPP
#define BOOST_MATH_SF_DETAIL_INV_T_HPP

#ifdef _MSC_VER
#pragma once
#endif

#include <boost/math/special_functions/cbrt.hpp>
#include <boost/math/special_functions/round.hpp>
#include <boost/math/special_functions/trunc.hpp>

namespace boost{ namespace math{ namespace detail{

//
// The main method used is due to Hill:
//
// G. W. Hill, Algorithm 396, Student's t-Quantiles,
// Communications of the ACM, 13(10): 619-620, Oct., 1970.
//
template <class T, class Policy>
T inverse_students_t_hill(T ndf, T u, const Policy& pol)
{
   BOOST_MATH_STD_USING
   BOOST_ASSERT(u <= 0.5);

   T a, b, c, d, q, x, y;

   if (ndf > 1e20f)
      return -boost::math::erfc_inv(2 * u, pol) * constants::root_two<T>();

   a = 1 / (ndf - 0.5f);
   b = 48 / (a * a);
   c = ((20700 * a / b - 98) * a - 16) * a + 96.36f;
   d = ((94.5f / (b + c) - 3) / b + 1) * sqrt(a * constants::pi<T>() / 2) * ndf;
   y = pow(d * 2 * u, 2 / ndf);

   if (y > (0.05f + a))
   {
      //
      // Asymptotic inverse expansion about normal:
      //
      x = -boost::math::erfc_inv(2 * u, pol) * constants::root_two<T>();
      y = x * x;

      if (ndf < 5)
         c += 0.3f * (ndf - 4.5f) * (x + 0.6f);
      c += (((0.05f * d * x - 5) * x - 7) * x - 2) * x + b;
      y = (((((0.4f * y + 6.3f) * y + 36) * y + 94.5f) / c - y - 3) / b + 1) * x;
      y = boost::math::expm1(a * y * y, pol);
   }
   else
   {
      y = ((1 / (((ndf + 6) / (ndf * y) - 0.089f * d - 0.822f)
              * (ndf + 2) * 3) + 0.5 / (ndf + 4)) * y - 1)
              * (ndf + 1) / (ndf + 2) + 1 / y;
   }
   q = sqrt(ndf * y);

   return -q;
}
//
// Tail and body series are due to Shaw:
//
// www.mth.kcl.ac.uk/~shaww/web_page/papers/Tdistribution06.pdf
//
// Shaw, W.T., 2006, "Sampling Student's T distribution - use of
// the inverse cumulative distribution function."
// Journal of Computational Finance, Vol 9 Issue 4, pp 37-73, Summer 2006
//
template <class T, class Policy>
T inverse_students_t_tail_series(T df, T v, const Policy& pol)
{
   BOOST_MATH_STD_USING
   // Tail series expansion, see section 6 of Shaw's paper.
   // w is calculated using Eq 60:
   T w = boost::math::tgamma_delta_ratio(df / 2, constants::half<T>(), pol)
      * sqrt(df * constants::pi<T>()) * v;
   // define some variables:
   T np2 = df + 2;
   T np4 = df + 4;
   T np6 = df + 6;
   //
   // Calculate the coefficients d(k), these depend only on the
   // number of degrees of freedom df, so at least in theory
   // we could tabulate these for fixed df, see p15 of Shaw:
   //
   T d[7] = { 1, };
   d[1] = -(df + 1) / (2 * np2);
   np2 *= (df + 2);
   d[2] = -df * (df + 1) * (df + 3) / (8 * np2 * np4);
   np2 *= df + 2;
   d[3] = -df * (df + 1) * (df + 5) * (((3 * df) + 7) * df -2) / (48 * np2 * np4 * np6);
   np2 *= (df + 2);
   np4 *= (df + 4);
   d[4] = -df * (df + 1) * (df + 7) *
      ( (((((15 * df) + 154) * df + 465) * df + 286) * df - 336) * df + 64 )
      / (384 * np2 * np4 * np6 * (df + 8));
   np2 *= (df + 2);
   d[5] = -df * (df + 1) * (df + 3) * (df + 9)
            * (((((((35 * df + 452) * df + 1573) * df + 600) * df - 2020) * df) + 928) * df -128)
            / (1280 * np2 * np4 * np6 * (df + 8) * (df + 10));
   np2 *= (df + 2);
   np4 *= (df + 4);
   np6 *= (df + 6);
   d[6] = -df * (df + 1) * (df + 11)
            * ((((((((((((945 * df) + 31506) * df + 425858) * df + 2980236) * df + 11266745) * df + 20675018) * df + 7747124) * df - 22574632) * df - 8565600) * df + 18108416) * df - 7099392) * df + 884736)
            / (46080 * np2 * np4 * np6 * (df + 8) * (df + 10) * (df +12));
   //
   // Now bring everthing together to provide the result,
   // this is Eq 62 of Shaw:
   //
   T rn = sqrt(df);
   T div = pow(rn * w, 1 / df);
   T power = div * div;
   T result = tools::evaluate_polynomial(d, power);
   result *= rn;
   result /= div;
   return -result;
}

template <class T, class Policy>
T inverse_students_t_body_series(T df, T u, const Policy& pol)
{
   BOOST_MATH_STD_USING
   //
   // Body series for small N:
   //
   // Start with Eq 56 of Shaw:
   //
   T v = boost::math::tgamma_delta_ratio(df / 2, constants::half<T>(), pol)
      * sqrt(df * constants::pi<T>()) * (u - constants::half<T>());
   //
   // Workspace for the polynomial coefficients:
   //
   T c[11] = { 0, 1, };
   //
   // Figure out what the coefficients are, note these depend
   // only on the degrees of freedom (Eq 57 of Shaw):
   //
   c[2] = T(1) / 6 + T(1) / (6 * df);
   T in = 1 / df;
   c[3] = (((T(1) / 120) * in) + (T(1) / 15)) * in + (T(7) / 120);
   c[4] = ((((T(1) / 5040) * in + (T(1) / 560)) * in + (T(3) / 112)) * in + T(127) / 5040);
   c[5] = ((((T(1) / 362880) * in + (T(17) / 45360)) * in + (T(67) / 60480)) * in + (T(479) / 45360)) * in + (T(4369) / 362880);
   c[6] = ((((((T(1) / 39916800) * in + (T(2503) / 39916800)) * in + (T(11867) / 19958400)) * in + (T(1285) / 798336)) * in + (T(153161) / 39916800)) * in + (T(34807) / 5702400));
   c[7] = (((((((T(1) / 6227020800LL) * in + (T(37) / 2402400)) * in +
      (T(339929) / 2075673600LL)) * in + (T(67217) / 97297200)) * in +
      (T(870341) / 691891200LL)) * in + (T(70691) / 64864800LL)) * in +
      (T(20036983LL) / 6227020800LL));
   c[8] = (((((((T(1) / 1307674368000LL) * in + (T(1042243LL) / 261534873600LL)) * in +
            (T(21470159) / 435891456000LL)) * in + (T(326228899LL) / 1307674368000LL)) * in +
            (T(843620579) / 1307674368000LL)) * in + (T(332346031LL) / 435891456000LL)) * in +
            (T(43847599) / 1307674368000LL)) * in + (T(2280356863LL) / 1307674368000LL);
   c[9] = (((((((((T(1) / 355687428096000LL)) * in + (T(24262727LL) / 22230464256000LL)) * in +
            (T(123706507LL) / 8083805184000LL)) * in + (T(404003599LL) / 4446092851200LL)) * in +
            (T(51811946317LL) / 177843714048000LL)) * in + (T(91423417LL) / 177843714048LL)) * in +
            (T(32285445833LL) / 88921857024000LL)) * in + (T(531839683LL) / 1710035712000LL)) * in +
            (T(49020204823LL) / 50812489728000LL);
   c[10] = (((((((((T(1) / 121645100408832000LL) * in +
            (T(4222378423LL) / 13516122267648000LL)) * in +
            (T(49573465457LL) / 10137091700736000LL)) * in +
            (T(176126809LL) / 5304600576000LL)) * in +
            (T(44978231873LL) / 355687428096000LL)) * in +
            (T(5816850595639LL) / 20274183401472000LL)) * in +
            (T(73989712601LL) / 206879422464000LL)) * in +
            (T(26591354017LL) / 259925428224000LL)) * in +
            (T(14979648446341LL) / 40548366802944000LL)) * in +
            (T(65967241200001LL) / 121645100408832000LL);
   //
   // The result is then a polynomial in v (see Eq 56 of Shaw):
   //
   return tools::evaluate_odd_polynomial(c, v);
}

template <class T, class Policy>
T inverse_students_t(T df, T u, T v, const Policy& pol, bool* pexact = 0)
{
   //
   // df = number of degrees of freedom.
   // u = probablity.
   // v = 1 - u.
   // l = lanczos type to use.
   //
   BOOST_MATH_STD_USING
   bool invert = false;
   T result = 0;
   if(pexact)
      *pexact = false;
   if(u > v)
   {
      // function is symmetric, invert it:
      std::swap(u, v);
      invert = true;
   }
   if((floor(df) == df) && (df < 20))
   {
      //
      // we have integer degrees of freedom, try for the special
      // cases first:
      //
      T tolerance = ldexp(1.0f, (2 * policies::digits<T, Policy>()) / 3);

      switch(itrunc(df, Policy()))
      {
      case 1:
         {
            //
            // df = 1 is the same as the Cauchy distribution, see
            // Shaw Eq 35:
            //
            if(u == 0.5)
               result = 0;
            else
               result = -cos(constants::pi<T>() * u) / sin(constants::pi<T>() * u);
            if(pexact)
               *pexact = true;
            break;
         }
      case 2:
         {
            //
            // df = 2 has an exact result, see Shaw Eq 36:
            //
            result =(2 * u - 1) / sqrt(2 * u * v);
            if(pexact)
               *pexact = true;
            break;
         }
      case 4:
         {
            //
            // df = 4 has an exact result, see Shaw Eq 38 & 39:
            //
            T alpha = 4 * u * v;
            T root_alpha = sqrt(alpha);
            T r = 4 * cos(acos(root_alpha) / 3) / root_alpha;
            T x = sqrt(r - 4);
            result = u - 0.5f < 0 ? -x : x;
            if(pexact)
               *pexact = true;
            break;
         }
      case 6:
         {
            //
            // We get numeric overflow in this area:
            //
            if(u < 1e-150)
               return (invert ? -1 : 1) * inverse_students_t_hill(df, u, pol);
            //
            // Newton-Raphson iteration of a polynomial case,
            // choice of seed value is taken from Shaw's online
            // supplement: