| // 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<7, T, T>(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): |
| // |
| T in = 1 / df; |
| c[2] = 0.16666666666666666667 + 0.16666666666666666667 * in; |
| c[3] = (0.0083333333333333333333 * in |
| + 0.066666666666666666667) * in |
| + 0.058333333333333333333; |
| c[4] = ((0.00019841269841269841270 * in |
| + 0.0017857142857142857143) * in |
| + 0.026785714285714285714) * in |
| + 0.025198412698412698413; |
| c[5] = (((2.7557319223985890653e-6 * in |
| + 0.00037477954144620811287) * in |
| - 0.0011078042328042328042) * in |
| + 0.010559964726631393298) * in |
| + 0.012039792768959435626; |
| c[6] = ((((2.5052108385441718775e-8 * in |
| - 0.000062705427288760622094) * in |
| + 0.00059458674042007375341) * in |
| - 0.0016095979637646304313) * in |
| + 0.0061039211560044893378) * in |
| + 0.0038370059724226390893; |
| c[7] = (((((1.6059043836821614599e-10 * in |
| + 0.000015401265401265401265) * in |
| - 0.00016376804137220803887) * in |
| + 0.00069084207973096861986) * in |
| - 0.0012579159844784844785) * in |
| + 0.0010898206731540064873) * in |
| + 0.0032177478835464946576; |
| c[8] = ((((((7.6471637318198164759e-13 * in |
| - 3.9851014346715404916e-6) * in |
| + 0.000049255746366361445727) * in |
| - 0.00024947258047043099953) * in |
| + 0.00064513046951456342991) * in |
| - 0.00076245135440323932387) * in |
| + 0.000033530976880017885309) * in |
| + 0.0017438262298340009980; |
| c[9] = (((((((2.8114572543455207632e-15 * in |
| + 1.0914179173496789432e-6) * in |
| - 0.000015303004486655377567) * in |
| + 0.000090867107935219902229) * in |
| - 0.00029133414466938067350) * in |
| + 0.00051406605788341121363) * in |
| - 0.00036307660358786885787) * in |
| - 0.00031101086326318780412) * in |
| + 0.00096472747321388644237; |
| c[10] = ((((((((8.2206352466243297170e-18 * in |
| - 3.1239569599829868045e-7) * in |
| + 4.8903045291975346210e-6) * in |
| - 0.000033202652391372058698) * in |
| + 0.00012645437628698076975) * in |
| - 0.00028690924218514613987) * in |
| + 0.00035764655430568632777) * in |
| - 0.00010230378073700412687) * in |
| - 0.00036942667800009661203) * in |
| + 0.00054229262813129686486; |
| // |
| // The result is then a polynomial in v (see Eq 56 of Shaw): |
| // |
| return tools::evaluate_odd_polynomial<11, T, T>(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 ? (T)-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: |
| // |
| 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.85498797333834849467655443627193; |
| 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 ? (T)-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, pol); |
| } |
| else |
| { |
| result = boost::math::detail::inverse_students_t_tail_series(df, u, pol); |
| } |
| } |
| 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, iround(T(df / -0.654f), pol)); |
| if(u > crossover) |
| { |
| result = boost::math::detail::inverse_students_t_hill(df, u, pol); |
| } |
| else |
| { |
| result = boost::math::detail::inverse_students_t_tail_series(df, u, pol); |
| } |
| } |
| } |
| return invert ? (T)-result : result; |
| } |
| |
| template <class T, class Policy> |
| inline T find_ibeta_inv_from_t_dist(T a, T p, T q, T* py, const Policy& pol) |
| { |
| T u = (p > q) ? T(0.5f - q) / T(2) : T(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, pol); |
| T x = df / (df + t * t); |
| *py = t * t / (df + t * t); |
| return x; |
| } |
| |
| template <class T, class Policy> |
| inline T fast_students_t_quantile_imp(T df, T p, const Policy& pol, const mpl::false_*) |
| { |
| BOOST_MATH_STD_USING |
| // |
| // 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(0); |
| x = ibeta_inv(df / 2, T(0.5), 2 * probability, &y, pol); |
| if(df * y > tools::max_value<T>() * x) |
| t = policies::raise_overflow_error<T>("boost::math::students_t_quantile<%1%>(%1%,%1%)", 0, pol); |
| 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 Policy> |
| T fast_students_t_quantile_imp(T df, T p, const Policy& pol, const mpl::true_*) |
| { |
| BOOST_MATH_STD_USING |
| bool invert = false; |
| if((df < 2) && (floor(df) != df)) |
| return boost::math::detail::fast_students_t_quantile_imp(df, p, pol, 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, pol, &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, pol, false, true, &f1) |
| : ibeta_imp(constants::half<T>(), a, y, pol, 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, class Policy> |
| inline T fast_students_t_quantile(T df, T p, const Policy& pol) |
| { |
| typedef typename policies::evaluation<T, Policy>::type value_type; |
| typedef typename policies::normalise< |
| Policy, |
| policies::promote_float<false>, |
| policies::promote_double<false>, |
| policies::discrete_quantile<>, |
| policies::assert_undefined<> >::type forwarding_policy; |
| |
| typedef mpl::bool_< |
| (std::numeric_limits<T>::digits <= 53) |
| && |
| (std::numeric_limits<T>::is_specialized)> tag_type; |
| 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%)"); |
| } |
| |
| }}} // namespaces |
| |
| #endif // BOOST_MATH_SF_DETAIL_INV_T_HPP |
| |
| |
| |