Boost logo

Boost-Commit :

From: johnmaddock_at_[hidden]
Date: 2007-06-07 13:42:23


Author: johnmaddock
Date: 2007-06-07 13:42:22 EDT (Thu, 07 Jun 2007)
New Revision: 4489
URL: http://svn.boost.org/trac/boost/changeset/4489

Log:
Performance optimisations: added Cornish-Fisher style expansions to give much better initial approximations to the inverses. Optimised Student's T quantile so that only one Halley step is taken at double precision.

Added:
   sandbox/math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp
Text files modified:
   sandbox/math_toolkit/boost/math/distributions/binomial.hpp | 109 +++++++++++++++++++++++++++++++++++----
   sandbox/math_toolkit/boost/math/distributions/students_t.hpp | 25 ++++++--
   sandbox/math_toolkit/boost/math/special_functions/detail/gamma_inva.hpp | 65 ++++++++++++++++++++++-
   sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inv_ab.hpp | 80 +++++++++++++++++++++++++++--
   sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inverse.hpp | 21 ++++++
   5 files changed, 270 insertions(+), 30 deletions(-)

Modified: sandbox/math_toolkit/boost/math/distributions/binomial.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/binomial.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/binomial.hpp 2007-06-07 13:42:22 EDT (Thu, 07 Jun 2007)
@@ -160,6 +160,42 @@
               return false;
            return true;
         }
+
+ template <class T>
+ T inverse_binomial_cornish_fisher(T n, T sf, T p, T q)
+ {
+ using namespace std;
+ // mean:
+ T m = n * sf;
+ // standard deviation:
+ T sigma = sqrt(n * sf * (1 - sf));
+ // skewness
+ T sk = (1 - 2 * sf) / sigma;
+ // kurtosis:
+ // T k = (1 - 6 * sf * (1 - sf) ) / (n * sf * (1 - sf));
+ // Get the inverse of a std normal distribution:
+ T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p) * constants::root_two<T>();
+ // Set the sign:
+ if(p < 0.5)
+ x = -x;
+ T x2 = x * x;
+ // w is correction term due to skewness
+ T w = x + sk * (x2 - 1) / 6;
+ /*
+ // Add on correction due to kurtosis.
+ // Disabled for now, seems to make things worse?
+ //
+ if(n >= 10)
+ w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;
+ */
+ w = m + sigma * w;
+ if(w < tools::min_value<T>())
+ return sqrt(tools::min_value<T>());
+ if(w > n)
+ return n;
+ return w;
+ }
+
      }
 
     // typedef binomial_distribution<double> binomial;
@@ -566,9 +602,10 @@
         // Error checks:
         using namespace std; // ADL of std names
         RealType result;
+ RealType trials = dist.trials();
         if(false == binomial_detail::check_dist_and_prob(
            BOOST_CURRENT_FUNCTION,
- dist.trials(),
+ trials,
            dist.success_fraction(),
            p,
            &result))
@@ -589,20 +626,40 @@
            // so n is the most sensible answer here:
            return dist.trials();
         }
- if (p <= pow(1 - dist.success_fraction(), dist.trials()))
+ if (p <= pow(1 - dist.success_fraction(), trials))
         { // p <= pdf(dist, 0) == cdf(dist, 0)
                                         return 0; // So the only reasonable result is zero.
                                 } // And root finder would fail otherwise.
 
         // Solve for quantile numerically:
         //
+ RealType guess = binomial_detail::inverse_binomial_cornish_fisher(trials, dist.success_fraction(), p, 1-p);
+ RealType factor = 8;
+ if(trials > 100)
+ factor = 1.01f; // guess is pretty accurate
+ else if((trials > 10) && (trials - 1 > guess) && (guess > 3))
+ factor = 1.15f; // less accurate but OK.
+ else if(trials < 10)
+ {
+ // pretty inaccurate guess in this area:
+ if(guess > trials / 64)
+ {
+ guess = trials / 4;
+ factor = 2;
+ }
+ else
+ guess = trials / 1024;
+ }
+ else
+ factor = 2; // trials largish, but in far tails.
+
         detail::binomial_functor<RealType> f(dist, p);
         tools::eps_tolerance<RealType> tol(tools::digits<RealType>());
         boost::uintmax_t max_iter = 1000;
         std::pair<RealType, RealType> r = tools::bracket_and_solve_root(
            f,
- dist.trials() / 2,
- static_cast<RealType>(8),
+ guess,
+ factor,
            true,
            tol,
            max_iter);
@@ -621,10 +678,11 @@
         // Error checks:
         RealType q = c.param;
         const binomial_distribution<RealType>& dist = c.dist;
+ RealType trials = dist.trials();
         RealType result;
         if(false == binomial_detail::check_dist_and_prob(
            BOOST_CURRENT_FUNCTION,
- dist.trials(),
+ trials,
            dist.success_fraction(),
            q,
            &result))
@@ -645,23 +703,50 @@
            // so n is the most sensible answer here:
            return dist.trials();
         }
-
- if (-q <= powm1(1 - dist.success_fraction(), dist.trials()))
- { // // q <= cdf(complement(dist, 0)) == pdf(dist, 0)
- return 0; // So the only reasonable result is zero.
- } // And root finder would fail otherwise.
+ if(trials == 1)
+ {
+ if(-q <= -dist.success_fraction())
+ return 0;
+ }
+ else
+ {
+ if (-q <= expm1(log1p(-dist.success_fraction()) * trials))
+ { // // q <= cdf(complement(dist, 0)) == pdf(dist, 0)
+ return 0; // So the only reasonable result is zero.
+ } // And root finder would fail otherwise.
+ }
 
         // Need to consider the case where
         //
         // Solve for quantile numerically:
         //
+ RealType guess = binomial_detail::inverse_binomial_cornish_fisher(trials, dist.success_fraction(), 1-q, q);
+ RealType factor = 8;
+ if(trials > 100)
+ factor = 1.01f; // guess is pretty accurate
+ else if((trials > 10) && (trials - 1 > guess) && (guess > 3))
+ factor = 1.15f; // less accurate but OK.
+ else if(trials < 10)
+ {
+ // pretty inaccurate guess in this area:
+ if(guess > trials / 64)
+ {
+ guess = trials / 4;
+ factor = 2;
+ }
+ else
+ guess = trials / 1024;
+ }
+ else
+ factor = 2; // trials largish, but in far tails.
+
         detail::binomial_functor<RealType> f(dist, q, true);
         tools::eps_tolerance<RealType> tol(tools::digits<RealType>());
         boost::uintmax_t max_iter = 1000;
         std::pair<RealType, RealType> r = tools::bracket_and_solve_root(
            f,
- dist.trials() / 2,
- static_cast<RealType>(8),
+ guess,
+ factor,
            true,
            tol,
            max_iter);

Modified: sandbox/math_toolkit/boost/math/distributions/students_t.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/students_t.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/students_t.hpp 2007-06-07 13:42:22 EDT (Thu, 07 Jun 2007)
@@ -85,11 +85,9 @@
    if(false == detail::check_df(
          BOOST_CURRENT_FUNCTION, degrees_of_freedom, &error_result))
       return error_result;
- // Might conceivably permit df = +infinity and use normal distribution.
- // TODO fails for t == 0 and df >=1e16 for ALL fp types.
- // - probably need to use normal distribution - when available.
- RealType basem1 = t * t / degrees_of_freedom;
+ // Might conceivably permit df = +infinity and use normal distribution.
    RealType result;
+ RealType basem1 = t * t / degrees_of_freedom;
    if(basem1 < 0.125)
    {
       result = exp(-boost::math::log1p(basem1) * (1+degrees_of_freedom) / 2);
@@ -99,7 +97,6 @@
       result = pow(1 / (1 + basem1), (degrees_of_freedom + 1) / 2);
    }
    result /= sqrt(degrees_of_freedom) * boost::math::beta(degrees_of_freedom / 2, RealType(0.5f));
-
    return result;
 } // pdf
 
@@ -139,12 +136,12 @@
    RealType probability;
    if(degrees_of_freedom > 2 * t2)
    {
- RealType z = t2 / (degrees_of_freedom + t * t);
+ RealType z = t2 / (degrees_of_freedom + t2);
       probability = ibetac(static_cast<RealType>(0.5), degrees_of_freedom / 2, z) / 2;
    }
    else
    {
- RealType z = degrees_of_freedom / (degrees_of_freedom + t * t);
+ RealType z = degrees_of_freedom / (degrees_of_freedom + t2);
       probability = ibeta(degrees_of_freedom / 2, static_cast<RealType>(0.5), z) / 2;
    }
    // Check 0 <= probability probability <= 1.
@@ -189,6 +186,11 @@
    if (probability == static_cast<RealType>(0.5))
      return 0;
    //
+ // This next block is disabled in favour of a faster method than
+ // incomplete beta inverse, code retained for future reference:
+ //
+#if 0
+ //
    // Calculate quantile of Student's t using the incomplete beta function inverse:
    //
    probability = (probability > 0.5) ? 1 - probability : probability;
@@ -205,6 +207,15 @@
       t = -t;
 
    return t;
+#endif
+ //
+ // Depending on how many digits RealType has, this may forward
+ // to the incomplete beta inverse as above. Otherwise uses a
+ // faster method that is accurate to ~15 digits everywhere
+ // and a couple of epsilon at double precision and in the central
+ // region where most use cases will occur...
+ //
+ return boost::math::detail::fast_students_t_quantile(degrees_of_freedom, probability);
 } // quantile
 
 template <class RealType>

Modified: sandbox/math_toolkit/boost/math/special_functions/detail/gamma_inva.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/gamma_inva.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/gamma_inva.hpp 2007-06-07 13:42:22 EDT (Thu, 07 Jun 2007)
@@ -32,6 +32,37 @@
 };
 
 template <class T>
+T inverse_poisson_cornish_fisher(T lambda, T p, T q)
+{
+ using namespace std;
+ // mean:
+ T m = lambda;
+ // standard deviation:
+ T sigma = sqrt(lambda);
+ // skewness
+ T sk = 1 / sigma;
+ // kurtosis:
+ // T k = 1/lambda;
+ // Get the inverse of a std normal distribution:
+ T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p) * constants::root_two<T>();
+ // Set the sign:
+ if(p < 0.5)
+ x = -x;
+ T x2 = x * x;
+ // w is correction term due to skewness
+ T w = x + sk * (x2 - 1) / 6;
+ /*
+ // Add on correction due to kurtosis.
+ // Disabled for now, seems to make things worse?
+ //
+ if(lambda >= 10)
+ w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;
+ */
+ w = m + sigma * w;
+ return w > tools::min_value<T>() ? w : tools::min_value<T>();
+}
+
+template <class T>
 T gamma_inva_imp(const T& z, const T& p, const T& q)
 {
    using namespace std; // for ADL of std lib math functions
@@ -63,9 +94,31 @@
    // to bracket the root from there:
    //
    T guess;
- if(z > 1.1)
+ T factor = 8;
+ if(z >= 1)
    {
- guess = z;
+ //
+ // We can use the relationship between the incomplete
+ // gamma function and the poisson distribution to
+ // calculate an approximate inverse, for large z
+ // this is actually pretty accurate, but it fails badly
+ // when z is very small. Also set our step-factor according
+ // to how accurate we think the result is likely to be:
+ //
+ guess = 1 + inverse_poisson_cornish_fisher(z, q, p);
+ if(z > 5)
+ {
+ if(z > 1000)
+ factor = 1.01f;
+ else if(z > 50)
+ factor = 1.1f;
+ else if(guess > 10)
+ factor = 1.25f;
+ else
+ factor = 2;
+ if(guess < 1.1)
+ factor = 8;
+ }
    }
    else if(z > 0.5)
    {
@@ -79,7 +132,13 @@
    // Max iterations permitted:
    //
    boost::uintmax_t max_iter = 200;
- std::pair<T, T> r = bracket_and_solve_root(f, guess, static_cast<T>(2), false, tol, max_iter);
+ //
+ // Use our generic derivative-free root finding procedure.
+ // We could use Newton steps here, taking the PDF of the
+ // Poisson distribution as our derivative, but that's
+ // even worse performance-wise than the generic method :-(
+ //
+ std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, false, tol, max_iter);
    if(max_iter >= 200)
       tools::logic_error<T>(BOOST_CURRENT_FUNCTION, "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first);
    return (r.first + r.second) / 2;

Modified: sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inv_ab.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inv_ab.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inv_ab.hpp 2007-06-07 13:42:22 EDT (Thu, 07 Jun 2007)
@@ -34,6 +34,39 @@
 };
 
 template <class T>
+T inverse_negative_binomial_cornish_fisher(T n, T sf, T sfc, T p, T q)
+{
+ using namespace std;
+ // mean:
+ T m = n * (sfc) / sf;
+ T t = sqrt(n * (sfc));
+ // standard deviation:
+ T sigma = t / sf;
+ // skewness
+ T sk = (1 + sfc) / t;
+ // kurtosis:
+ T k = (6 - sf * (5+sfc)) / (n * (sfc));
+ // Get the inverse of a std normal distribution:
+ T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p) * constants::root_two<T>();
+ // Set the sign:
+ if(p < 0.5)
+ x = -x;
+ T x2 = x * x;
+ // w is correction term due to skewness
+ T w = x + sk * (x2 - 1) / 6;
+ //
+ // Add on correction due to kurtosis.
+ //
+ if(n >= 10)
+ w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;
+
+ w = m + sigma * w;
+ if(w < tools::min_value<T>())
+ return tools::min_value<T>();
+ return w;
+}
+
+template <class T>
 T ibeta_inv_ab_imp(const T& b, const T& z, const T& p, const T& q, bool swap_ab)
 {
    using namespace std; // for ADL of std lib math functions
@@ -65,21 +98,56 @@
    // are quite sensitive so we should need too many steps
    // to bracket the root from there:
    //
- T guess;
- if((p < q) != swap_ab)
+ T guess = 0;
+ T factor = 5;
+ //
+ // Convert variables to parameters of a negative binomial distribution:
+ //
+ T n = b;
+ T sf = swap_ab ? z : 1-z;
+ T sfc = swap_ab ? 1-z : z;
+ T u = swap_ab ? p : q;
+ T v = swap_ab ? q : p;
+ if(u <= pow(sf, n))
    {
- guess = b * 2;
+ //
+ // Result is less than 1, negative binomial approximation
+ // is useless....
+ //
+ if((p < q) != swap_ab)
+ {
+ guess = (std::min)(b * 2, T(1));
+ }
+ else
+ {
+ guess = (std::min)(b / 2, T(1));
+ }
    }
- else
+ if(n * n * n * u * sf > 0.005)
+ guess = 1 + inverse_negative_binomial_cornish_fisher(n, sf, sfc, u, v);
+
+ if(guess < 10)
    {
- guess = b / 2;
+ //
+ // Negative binomial approximation not accurate in this area:
+ //
+ if((p < q) != swap_ab)
+ {
+ guess = (std::min)(b * 2, T(10));
+ }
+ else
+ {
+ guess = (std::min)(b / 2, T(10));
+ }
    }
+ else
+ factor = (v < sqrt(tools::epsilon<T>())) ? 2 : (guess < 20 ? 1.2f : 1.1f);
    BOOST_MATH_INSTRUMENT_CODE("guess = " << guess);
    //
    // Max iterations permitted:
    //
    boost::uintmax_t max_iter = 200;
- std::pair<T, T> r = bracket_and_solve_root(f, guess, static_cast<T>(5), swap_ab ? true : false, tol, max_iter);
+ std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, swap_ab ? true : false, tol, max_iter);
    if(max_iter >= 200)
       tools::logic_error<T>(BOOST_CURRENT_FUNCTION, "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first);
    return (r.first + r.second) / 2;

Modified: sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inverse.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inverse.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inverse.hpp 2007-06-07 13:42:22 EDT (Thu, 07 Jun 2007)
@@ -9,6 +9,7 @@
 #include <boost/math/special_functions/beta.hpp>
 #include <boost/math/special_functions/erf.hpp>
 #include <boost/math/tools/roots.hpp>
+#include <boost/math/special_functions/detail/t_distribution_inv.hpp>
 
 namespace boost{ namespace math{ namespace detail{
 
@@ -444,7 +445,7 @@
 };
 
 template <class T, class L>
-T ibeta_inv_imp(T a, T b, T p, T q, const L& /* l */, T* py)
+T ibeta_inv_imp(T a, T b, T p, T q, const L& l, T* py)
 {
    using namespace std; // For ADL of math functions.
 
@@ -465,6 +466,16 @@
    T lower = 0;
    T upper = 1;
    //
+ // Student's T with b = 0.5 gets handled as a special case, swap
+ // around if the arguments are in the "wrong" order:
+ //
+ if(a == 0.5f)
+ {
+ std::swap(a, b);
+ std::swap(p, q);
+ invert = !invert;
+ }
+ //
    // Handle trivial cases first:
    //
    if(q == 0)
@@ -482,6 +493,12 @@
       if(py) *py = 1 - p;
       return p;
    }
+ else if((b == 0.5f) && (a >= 0.5f))
+ {
+ //
+ // We have a Student's T distribution:
+ x = estimate_ibeta_inv_from_t_dist(a, p, q, &y, l);
+ }
    else if(a + b > 5)
    {
       //
@@ -734,7 +751,7 @@
    //
    if(lower == 0)
    {
- if(invert)
+ if(invert && (py == 0))
       {
          //
          // We're not interested in answers smaller than machine epsilon:

Added: sandbox/math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp 2007-06-07 13:42:22 EDT (Thu, 07 Jun 2007)
@@ -0,0 +1,502 @@
+// (C) Copyright John Maddock 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
+
+#include <boost/math/special_functions/cbrt.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>
+T inverse_students_t_hill(T ndf, T u)
+{
+ using namespace std;
+ BOOST_ASSERT(u <= 0.5);
+
+ T a, b, c, d, q, x, y;
+
+ if (ndf > 1e20f)
+ return -boost::math::erfc_inv(2 * u) * 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) * 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);
+ }
+ 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 L>
+T inverse_students_t_tail_series(T df, T v, T u, const L& l)
+{
+ using namespace std;
+ // Tail series expansion, see section 6 of Shaw's paper.
+ // w is calculated using Eq 60:
+ T w = detail::tgamma_delta_ratio_imp(df / 2, constants::half<T>(), l)
+ * 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 L>
+T inverse_students_t_body_series(T df, T u, const L& l)
+{
+ using namespace std;
+ //
+ // Body series for small N:
+ //
+ // Start with Eq 56 of Shaw:
+ //
+ T v = detail::tgamma_delta_ratio_imp(df / 2, constants::half<T>(), l)
+ * 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 L>
+T inverse_students_t(T df, T u, T v, const L& l, bool* pexact = 0)
+{
+ //
+ // df = number of degrees of freedom.
+ // u = probablity.
+ // v = 1 - u.
+ // l = lanczos type to use.
+ //
+ using namespace std;
+ 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 * tools::digits<T>()) / 3);
+
+ switch(boost::math::tools::real_cast<int>(df))
+ {
+ case 1:
+ {
+ //
+ // df = 1 is the same as the Cauchy distribution, see
+ // Shaw Eq 35:
+ //
+ // FIXME: fails when u is small!!!
+ result = tan(constants::pi<T>() * (u - constants::half<T>()));
+ 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);
+ //
+ // Newton-Raphson iteration of a polynomial case,
+ // choice of seed value is taken from Shaw's online
+ // supplement:
+ //
+ T a = 4 * (u - u * u);//1 - 4 * (u - 0.5f) * (u - 0.5f);
+ T b = boost::math::cbrt(a);
+ static const T c = 0.85498797333834849467655443627193L;
+ T p = 6 * (1 + c * (1 / b - 1));
+ T p0;
+ do{
+ T p2 = p * p;
+ T p4 = p2 * p2;
+ T p5 = p * p4;
+ p0 = p;
+ // next term is given by Eq 41:
+ p = 2 * (8 * a * p5 - 270 * p2 + 2187) / (5 * (4 * a * p4 - 216 * p - 243));
+ }while(fabs((p - p0) / p) > tolerance);
+ //
+ // Use Eq 45 to extract the result:
+ //
+ p = sqrt(p - df);
+ result = (u - 0.5f) < 0 ? -p : p;
+ break;
+ }
+#if 0
+ //
+ // These are Shaw's "exact" but iterative solutions
+ // for even df, the numerical accuracy of these is
+ // rather less than Hill's method, so these are disabled
+ // for now, which is a shame because they are reasonably
+ // quick to evaluate...
+ //
+ case 8:
+ {
+ //
+ // Newton-Raphson iteration of a polynomial case,
+ // choice of seed value is taken from Shaw's online
+ // supplement:
+ //
+ static const T c8 = 0.85994765706259820318168359251872L;
+ T a = 4 * (u - u * u); //1 - 4 * (u - 0.5f) * (u - 0.5f);
+ T b = pow(a, T(1) / 4);
+ T p = 8 * (1 + c8 * (1 / b - 1));
+ T p0 = p;
+ do{
+ T p5 = p * p;
+ p5 *= p5 * p;
+ p0 = p;
+ // Next term is given by Eq 42:
+ p = 2 * (3 * p + (640 * (160 + p * (24 + p * (p + 4)))) / (-5120 + p * (-2048 - 960 * p + a * p5))) / 7;
+ }while(fabs((p - p0) / p) > tolerance);
+ //
+ // Use Eq 45 to extract the result:
+ //
+ p = sqrt(p - df);
+ result = (u - 0.5f) < 0 ? -p : p;
+ break;
+ }
+ case 10:
+ {
+ //
+ // Newton-Raphson iteration of a polynomial case,
+ // choice of seed value is taken from Shaw's online
+ // supplement:
+ //
+ static const T c10 = 0.86781292867813396759105692122285L;
+ T a = 4 * (u - u * u); //1 - 4 * (u - 0.5f) * (u - 0.5f);
+ T b = pow(a, T(1) / 5);
+ T p = 10 * (1 + c10 * (1 / b - 1));
+ T p0;
+ do{
+ T p6 = p * p;
+ p6 *= p6 * p6;
+ p0 = p;
+ // Next term given by Eq 43:
+ p = (8 * p) / 9 + (218750 * (21875 + 4 * p * (625 + p * (75 + 2 * p * (5 + p))))) /
+ (9 * (-68359375 + 8 * p * (-2343750 + p * (-546875 - 175000 * p + 8 * a * p6))));
+ }while(fabs((p - p0) / p) > tolerance);
+ //
+ // Use Eq 45 to extract the result:
+ //
+ p = sqrt(p - df);
+ result = (u - 0.5f) < 0 ? -p : p;
+ break;
+ }
+#endif
+ default:
+ goto calculate_real;
+ }
+ }
+ else
+ {
+calculate_real:
+ if(df < 3)
+ {
+ //
+ // Use a roughly linear scheme to choose between Shaw's
+ // tail series and body series:
+ //
+ T crossover = 0.2742f - df * 0.0242143f;
+ if(u > crossover)
+ {
+ result = boost::math::detail::inverse_students_t_body_series(df, u, l);
+ }
+ else
+ {
+ result = boost::math::detail::inverse_students_t_tail_series(df, u, v, l);
+ }
+ }
+ else
+ {
+ //
+ // Use Hill's method except in the exteme tails
+ // where we use Shaw's tail series.
+ // The crossover point is roughly exponential in -df:
+ //
+ T crossover = ldexp(1.0f, tools::real_cast<int>(df / -0.654f));
+ if(u > crossover)
+ {
+ result = boost::math::detail::inverse_students_t_hill(df, u);
+ }
+ else
+ {
+ result = boost::math::detail::inverse_students_t_tail_series(df, u, v, l);
+ }
+ }
+ }
+ return invert ? -result : result;
+}
+
+template <class T, class L>
+inline T estimate_ibeta_inv_from_t_dist(T a, T p, T q, T* py, const L& l)
+{
+ T u = (p > q) ? 0.5f - q / 2 : p / 2;
+ T v = 1 - u; // u < 0.5 so no cancellation error
+ T df = a * 2;
+ T t = boost::math::detail::inverse_students_t(df, u, v, l);
+ T x = df / (df + t * t);
+ *py = t * t / (df + t * t);
+ return x;
+}
+
+template <class T, class L>
+inline T fast_students_t_quantile_imp(T df, T p, const L& l, const mpl::false_*)
+{
+ using namespace std;
+ //
+ // Need to use inverse incomplete beta to get
+ // required precision so not so fast:
+ //
+ T probability = (p > 0.5) ? 1 - p : p;
+ T t, x, y;
+ x = ibeta_inv(df / 2, T(0.5), 2 * probability, &y);
+ if(df * y > tools::max_value<T>() * x)
+ t = tools::overflow_error<T>(BOOST_CURRENT_FUNCTION);
+ else
+ t = sqrt(df * y / x);
+ //
+ // Figure out sign based on the size of p:
+ //
+ if(p < 0.5)
+ t = -t;
+ return t;
+}
+
+template <class T, class L>
+T fast_students_t_quantile_imp(T df, T p, const L& l, const mpl::true_*)
+{
+ using namespace std;
+ bool invert = false;
+ if((df < 2) && (floor(df) != df))
+ return boost::math::detail::fast_students_t_quantile_imp(df, p, l, static_cast<mpl::false_*>(0));
+ if(p > 0.5)
+ {
+ p = 1 - p;
+ invert = true;
+ }
+ //
+ // Get an estimate of the result:
+ //
+ bool exact;
+ T t = inverse_students_t(df, p, 1-p, l, &exact);
+ if((t == 0) || exact)
+ return invert ? -t : t; // can't do better!
+ //
+ // Change variables to inverse incomplete beta:
+ //
+ T t2 = t * t;
+ T xb = df / (df + t2);
+ T y = t2 / (df + t2);
+ T a = df / 2;
+ //
+ // t can be so large that x underflows,
+ // just return our estimate in that case:
+ //
+ if(xb == 0)
+ return t;
+ //
+ // Get incomplete beta and it's derivative:
+ //
+ T f1;
+ T f0 = xb < y ? ibeta_imp(a, constants::half<T>(), xb, l, false, true, &f1)
+ : ibeta_imp(constants::half<T>(), a, y, l, true, true, &f1);
+
+ // Get cdf from incomplete beta result:
+ T p0 = f0 / 2 - p;
+ // Get pdf from derivative:
+ T p1 = f1 * sqrt(y * xb * xb * xb / df);
+ //
+ // Second derivative divided by p1:
+ //
+ // yacas gives:
+ //
+ // In> PrettyForm(Simplify(D(t) (1 + t^2/v) ^ (-(v+1)/2)))
+ //
+ // | | v + 1 | |
+ // | -| ----- + 1 | |
+ // | | 2 | |
+ // -| | 2 | |
+ // | | t | |
+ // | | -- + 1 | |
+ // | ( v + 1 ) * | v | * t |
+ // ---------------------------------------------
+ // v
+ //
+ // Which after some manipulation is:
+ //
+ // -p1 * t * (df + 1) / (t^2 + df)
+ //
+ T p2 = t * (df + 1) / (t * t + df);
+ // Halley step:
+ t = fabs(t);
+ t += p0 / (p1 + p0 * p2 / 2);
+ return !invert ? -t : t;
+}
+
+template <class T>
+inline T fast_students_t_quantile(T df, T p)
+{
+ typedef typename lanczos::lanczos_traits<T>::value_type value_type;
+ typedef typename lanczos::lanczos_traits<T>::evaluation_type evaluation_type;
+ typedef mpl::bool_<
+ (std::numeric_limits<T>::digits <= 53)
+ &&
+ (std::numeric_limits<T>::is_specialized)> tag_type;
+ return tools::checked_narrowing_cast<T>(fast_students_t_quantile_imp(static_cast<value_type>(df), static_cast<value_type>(p), evaluation_type(), static_cast<tag_type*>(0)), BOOST_CURRENT_FUNCTION);
+}
+
+}}} // namespaces
+
+#endif // BOOST_MATH_SF_DETAIL_INV_T_HPP


Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk