1// Copyright John Maddock 2007.
2// Copyright Paul A. Bristow 2007
3// Use, modification and distribution are subject to the
4// Boost Software License, Version 1.0. (See accompanying file
5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
6
7#ifndef BOOST_MATH_SF_DETAIL_INV_T_HPP
8#define BOOST_MATH_SF_DETAIL_INV_T_HPP
9
10#ifdef _MSC_VER
11#pragma once
12#endif
13
14#include <boost/math/special_functions/cbrt.hpp>
15#include <boost/math/special_functions/round.hpp>
16#include <boost/math/special_functions/trunc.hpp>
17
18namespace boost{ namespace math{ namespace detail{
19
20//
21// The main method used is due to Hill:
22//
23// G. W. Hill, Algorithm 396, Student's t-Quantiles,
24// Communications of the ACM, 13(10): 619-620, Oct., 1970.
25//
26template <class T, class Policy>
27T inverse_students_t_hill(T ndf, T u, const Policy& pol)
28{
29 BOOST_MATH_STD_USING
30 BOOST_ASSERT(u <= 0.5);
31
32 T a, b, c, d, q, x, y;
33
34 if (ndf > 1e20f)
35 return -boost::math::erfc_inv(2 * u, pol) * constants::root_two<T>();
36
37 a = 1 / (ndf - 0.5f);
38 b = 48 / (a * a);
39 c = ((20700 * a / b - 98) * a - 16) * a + 96.36f;
40 d = ((94.5f / (b + c) - 3) / b + 1) * sqrt(a * constants::pi<T>() / 2) * ndf;
41 y = pow(d * 2 * u, 2 / ndf);
42
43 if (y > (0.05f + a))
44 {
45 //
46 // Asymptotic inverse expansion about normal:
47 //
48 x = -boost::math::erfc_inv(2 * u, pol) * constants::root_two<T>();
49 y = x * x;
50
51 if (ndf < 5)
52 c += 0.3f * (ndf - 4.5f) * (x + 0.6f);
53 c += (((0.05f * d * x - 5) * x - 7) * x - 2) * x + b;
54 y = (((((0.4f * y + 6.3f) * y + 36) * y + 94.5f) / c - y - 3) / b + 1) * x;
55 y = boost::math::expm1(a * y * y, pol);
56 }
57 else
58 {
59 y = static_cast<T>(((1 / (((ndf + 6) / (ndf * y) - 0.089f * d - 0.822f)
60 * (ndf + 2) * 3) + 0.5 / (ndf + 4)) * y - 1)
61 * (ndf + 1) / (ndf + 2) + 1 / y);
62 }
63 q = sqrt(ndf * y);
64
65 return -q;
66}
67//
68// Tail and body series are due to Shaw:
69//
70// www.mth.kcl.ac.uk/~shaww/web_page/papers/Tdistribution06.pdf
71//
72// Shaw, W.T., 2006, "Sampling Student's T distribution - use of
73// the inverse cumulative distribution function."
74// Journal of Computational Finance, Vol 9 Issue 4, pp 37-73, Summer 2006
75//
76template <class T, class Policy>
77T inverse_students_t_tail_series(T df, T v, const Policy& pol)
78{
79 BOOST_MATH_STD_USING
80 // Tail series expansion, see section 6 of Shaw's paper.
81 // w is calculated using Eq 60:
82 T w = boost::math::tgamma_delta_ratio(df / 2, constants::half<T>(), pol)
83 * sqrt(df * constants::pi<T>()) * v;
84 // define some variables:
85 T np2 = df + 2;
86 T np4 = df + 4;
87 T np6 = df + 6;
88 //
89 // Calculate the coefficients d(k), these depend only on the
90 // number of degrees of freedom df, so at least in theory
91 // we could tabulate these for fixed df, see p15 of Shaw:
92 //
93 T d[7] = { 1, };
94 d[1] = -(df + 1) / (2 * np2);
95 np2 *= (df + 2);
96 d[2] = -df * (df + 1) * (df + 3) / (8 * np2 * np4);
97 np2 *= df + 2;
98 d[3] = -df * (df + 1) * (df + 5) * (((3 * df) + 7) * df -2) / (48 * np2 * np4 * np6);
99 np2 *= (df + 2);
100 np4 *= (df + 4);
101 d[4] = -df * (df + 1) * (df + 7) *
102 ( (((((15 * df) + 154) * df + 465) * df + 286) * df - 336) * df + 64 )
103 / (384 * np2 * np4 * np6 * (df + 8));
104 np2 *= (df + 2);
105 d[5] = -df * (df + 1) * (df + 3) * (df + 9)
106 * (((((((35 * df + 452) * df + 1573) * df + 600) * df - 2020) * df) + 928) * df -128)
107 / (1280 * np2 * np4 * np6 * (df + 8) * (df + 10));
108 np2 *= (df + 2);
109 np4 *= (df + 4);
110 np6 *= (df + 6);
111 d[6] = -df * (df + 1) * (df + 11)
112 * ((((((((((((945 * df) + 31506) * df + 425858) * df + 2980236) * df + 11266745) * df + 20675018) * df + 7747124) * df - 22574632) * df - 8565600) * df + 18108416) * df - 7099392) * df + 884736)
113 / (46080 * np2 * np4 * np6 * (df + 8) * (df + 10) * (df +12));
114 //
115 // Now bring everything together to provide the result,
116 // this is Eq 62 of Shaw:
117 //
118 T rn = sqrt(df);
119 T div = pow(rn * w, 1 / df);
120 T power = div * div;
121 T result = tools::evaluate_polynomial<7, T, T>(d, power);
122 result *= rn;
123 result /= div;
124 return -result;
125}
126
127template <class T, class Policy>
128T inverse_students_t_body_series(T df, T u, const Policy& pol)
129{
130 BOOST_MATH_STD_USING
131 //
132 // Body series for small N:
133 //
134 // Start with Eq 56 of Shaw:
135 //
136 T v = boost::math::tgamma_delta_ratio(df / 2, constants::half<T>(), pol)
137 * sqrt(df * constants::pi<T>()) * (u - constants::half<T>());
138 //
139 // Workspace for the polynomial coefficients:
140 //
141 T c[11] = { 0, 1, };
142 //
143 // Figure out what the coefficients are, note these depend
144 // only on the degrees of freedom (Eq 57 of Shaw):
145 //
146 T in = 1 / df;
147 c[2] = static_cast<T>(0.16666666666666666667 + 0.16666666666666666667 * in);
148 c[3] = static_cast<T>((0.0083333333333333333333 * in
149 + 0.066666666666666666667) * in
150 + 0.058333333333333333333);
151 c[4] = static_cast<T>(((0.00019841269841269841270 * in
152 + 0.0017857142857142857143) * in
153 + 0.026785714285714285714) * in
154 + 0.025198412698412698413);
155 c[5] = static_cast<T>((((2.7557319223985890653e-6 * in
156 + 0.00037477954144620811287) * in
157 - 0.0011078042328042328042) * in
158 + 0.010559964726631393298) * in
159 + 0.012039792768959435626);
160 c[6] = static_cast<T>(((((2.5052108385441718775e-8 * in
161 - 0.000062705427288760622094) * in
162 + 0.00059458674042007375341) * in
163 - 0.0016095979637646304313) * in
164 + 0.0061039211560044893378) * in
165 + 0.0038370059724226390893);
166 c[7] = static_cast<T>((((((1.6059043836821614599e-10 * in
167 + 0.000015401265401265401265) * in
168 - 0.00016376804137220803887) * in
169 + 0.00069084207973096861986) * in
170 - 0.0012579159844784844785) * in
171 + 0.0010898206731540064873) * in
172 + 0.0032177478835464946576);
173 c[8] = static_cast<T>(((((((7.6471637318198164759e-13 * in
174 - 3.9851014346715404916e-6) * in
175 + 0.000049255746366361445727) * in
176 - 0.00024947258047043099953) * in
177 + 0.00064513046951456342991) * in
178 - 0.00076245135440323932387) * in
179 + 0.000033530976880017885309) * in
180 + 0.0017438262298340009980);
181 c[9] = static_cast<T>((((((((2.8114572543455207632e-15 * in
182 + 1.0914179173496789432e-6) * in
183 - 0.000015303004486655377567) * in
184 + 0.000090867107935219902229) * in
185 - 0.00029133414466938067350) * in
186 + 0.00051406605788341121363) * in
187 - 0.00036307660358786885787) * in
188 - 0.00031101086326318780412) * in
189 + 0.00096472747321388644237);
190 c[10] = static_cast<T>(((((((((8.2206352466243297170e-18 * in
191 - 3.1239569599829868045e-7) * in
192 + 4.8903045291975346210e-6) * in
193 - 0.000033202652391372058698) * in
194 + 0.00012645437628698076975) * in
195 - 0.00028690924218514613987) * in
196 + 0.00035764655430568632777) * in
197 - 0.00010230378073700412687) * in
198 - 0.00036942667800009661203) * in
199 + 0.00054229262813129686486);
200 //
201 // The result is then a polynomial in v (see Eq 56 of Shaw):
202 //
203 return tools::evaluate_odd_polynomial<11, T, T>(c, v);
204}
205
206template <class T, class Policy>
207T inverse_students_t(T df, T u, T v, const Policy& pol, bool* pexact = 0)
208{
209 //
210 // df = number of degrees of freedom.
211 // u = probability.
212 // v = 1 - u.
213 // l = lanczos type to use.
214 //
215 BOOST_MATH_STD_USING
216 bool invert = false;
217 T result = 0;
218 if(pexact)
219 *pexact = false;
220 if(u > v)
221 {
222 // function is symmetric, invert it:
223 std::swap(u, v);
224 invert = true;
225 }
226 if((floor(df) == df) && (df < 20))
227 {
228 //
229 // we have integer degrees of freedom, try for the special
230 // cases first:
231 //
232 T tolerance = ldexp(1.0f, (2 * policies::digits<T, Policy>()) / 3);
233
234 switch(itrunc(df, Policy()))
235 {
236 case 1:
237 {
238 //
239 // df = 1 is the same as the Cauchy distribution, see
240 // Shaw Eq 35:
241 //
242 if(u == 0.5)
243 result = 0;
244 else
245 result = -cos(constants::pi<T>() * u) / sin(constants::pi<T>() * u);
246 if(pexact)
247 *pexact = true;
248 break;
249 }
250 case 2:
251 {
252 //
253 // df = 2 has an exact result, see Shaw Eq 36:
254 //
255 result =(2 * u - 1) / sqrt(2 * u * v);
256 if(pexact)
257 *pexact = true;
258 break;
259 }
260 case 4:
261 {
262 //
263 // df = 4 has an exact result, see Shaw Eq 38 & 39:
264 //
265 T alpha = 4 * u * v;
266 T root_alpha = sqrt(alpha);
267 T r = 4 * cos(acos(root_alpha) / 3) / root_alpha;
268 T x = sqrt(r - 4);
269 result = u - 0.5f < 0 ? (T)-x : x;
270 if(pexact)
271 *pexact = true;
272 break;
273 }
274 case 6:
275 {
276 //
277 // We get numeric overflow in this area:
278 //
279 if(u < 1e-150)
280 return (invert ? -1 : 1) * inverse_students_t_hill(df, u, pol);
281 //
282 // Newton-Raphson iteration of a polynomial case,
283 // choice of seed value is taken from Shaw's online
284 // supplement:
285 //
286 T a = 4 * (u - u * u);//1 - 4 * (u - 0.5f) * (u - 0.5f);
287 T b = boost::math::cbrt(a);
288 static const T c = static_cast<T>(0.85498797333834849467655443627193);
289 T p = 6 * (1 + c * (1 / b - 1));
290 T p0;
291 do{
292 T p2 = p * p;
293 T p4 = p2 * p2;
294 T p5 = p * p4;
295 p0 = p;
296 // next term is given by Eq 41:
297 p = 2 * (8 * a * p5 - 270 * p2 + 2187) / (5 * (4 * a * p4 - 216 * p - 243));
298 }while(fabs((p - p0) / p) > tolerance);
299 //
300 // Use Eq 45 to extract the result:
301 //
302 p = sqrt(p - df);
303 result = (u - 0.5f) < 0 ? (T)-p : p;
304 break;
305 }
306#if 0
307 //
308 // These are Shaw's "exact" but iterative solutions
309 // for even df, the numerical accuracy of these is
310 // rather less than Hill's method, so these are disabled
311 // for now, which is a shame because they are reasonably
312 // quick to evaluate...
313 //
314 case 8:
315 {
316 //
317 // Newton-Raphson iteration of a polynomial case,
318 // choice of seed value is taken from Shaw's online
319 // supplement:
320 //
321 static const T c8 = 0.85994765706259820318168359251872L;
322 T a = 4 * (u - u * u); //1 - 4 * (u - 0.5f) * (u - 0.5f);
323 T b = pow(a, T(1) / 4);
324 T p = 8 * (1 + c8 * (1 / b - 1));
325 T p0 = p;
326 do{
327 T p5 = p * p;
328 p5 *= p5 * p;
329 p0 = p;
330 // Next term is given by Eq 42:
331 p = 2 * (3 * p + (640 * (160 + p * (24 + p * (p + 4)))) / (-5120 + p * (-2048 - 960 * p + a * p5))) / 7;
332 }while(fabs((p - p0) / p) > tolerance);
333 //
334 // Use Eq 45 to extract the result:
335 //
336 p = sqrt(p - df);
337 result = (u - 0.5f) < 0 ? -p : p;
338 break;
339 }
340 case 10:
341 {
342 //
343 // Newton-Raphson iteration of a polynomial case,
344 // choice of seed value is taken from Shaw's online
345 // supplement:
346 //
347 static const T c10 = 0.86781292867813396759105692122285L;
348 T a = 4 * (u - u * u); //1 - 4 * (u - 0.5f) * (u - 0.5f);
349 T b = pow(a, T(1) / 5);
350 T p = 10 * (1 + c10 * (1 / b - 1));
351 T p0;
352 do{
353 T p6 = p * p;
354 p6 *= p6 * p6;
355 p0 = p;
356 // Next term given by Eq 43:
357 p = (8 * p) / 9 + (218750 * (21875 + 4 * p * (625 + p * (75 + 2 * p * (5 + p))))) /
358 (9 * (-68359375 + 8 * p * (-2343750 + p * (-546875 - 175000 * p + 8 * a * p6))));
359 }while(fabs((p - p0) / p) > tolerance);
360 //
361 // Use Eq 45 to extract the result:
362 //
363 p = sqrt(p - df);
364 result = (u - 0.5f) < 0 ? -p : p;
365 break;
366 }
367#endif
368 default:
369 goto calculate_real;
370 }
371 }
372 else
373 {
374calculate_real:
375 if(df > 0x10000000)
376 {
377 result = -boost::math::erfc_inv(2 * u, pol) * constants::root_two<T>();
378 if((pexact) && (df >= 1e20))
379 *pexact = true;
380 }
381 else if(df < 3)
382 {
383 //
384 // Use a roughly linear scheme to choose between Shaw's
385 // tail series and body series:
386 //
387 T crossover = 0.2742f - df * 0.0242143f;
388 if(u > crossover)
389 {
390 result = boost::math::detail::inverse_students_t_body_series(df, u, pol);
391 }
392 else
393 {
394 result = boost::math::detail::inverse_students_t_tail_series(df, u, pol);
395 }
396 }
397 else
398 {
399 //
400 // Use Hill's method except in the extreme tails
401 // where we use Shaw's tail series.
402 // The crossover point is roughly exponential in -df:
403 //
404 T crossover = ldexp(1.0f, iround(T(df / -0.654f), typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type()));
405 if(u > crossover)
406 {
407 result = boost::math::detail::inverse_students_t_hill(df, u, pol);
408 }
409 else
410 {
411 result = boost::math::detail::inverse_students_t_tail_series(df, u, pol);
412 }
413 }
414 }
415 return invert ? (T)-result : result;
416}
417
418template <class T, class Policy>
419inline T find_ibeta_inv_from_t_dist(T a, T p, T /*q*/, T* py, const Policy& pol)
420{
421 T u = p / 2;
422 T v = 1 - u;
423 T df = a * 2;
424 T t = boost::math::detail::inverse_students_t(df, u, v, pol);
425 *py = t * t / (df + t * t);
426 return df / (df + t * t);
427}
428
429template <class T, class Policy>
430inline T fast_students_t_quantile_imp(T df, T p, const Policy& pol, const boost::false_type*)
431{
432 BOOST_MATH_STD_USING
433 //
434 // Need to use inverse incomplete beta to get
435 // required precision so not so fast:
436 //
437 T probability = (p > 0.5) ? 1 - p : p;
438 T t, x, y(0);
439 x = ibeta_inv(df / 2, T(0.5), 2 * probability, &y, pol);
440 if(df * y > tools::max_value<T>() * x)
441 t = policies::raise_overflow_error<T>("boost::math::students_t_quantile<%1%>(%1%,%1%)", 0, pol);
442 else
443 t = sqrt(df * y / x);
444 //
445 // Figure out sign based on the size of p:
446 //
447 if(p < 0.5)
448 t = -t;
449 return t;
450}
451
452template <class T, class Policy>
453T fast_students_t_quantile_imp(T df, T p, const Policy& pol, const boost::true_type*)
454{
455 BOOST_MATH_STD_USING
456 bool invert = false;
457 if((df < 2) && (floor(df) != df))
458 return boost::math::detail::fast_students_t_quantile_imp(df, p, pol, static_cast<boost::false_type*>(0));
459 if(p > 0.5)
460 {
461 p = 1 - p;
462 invert = true;
463 }
464 //
465 // Get an estimate of the result:
466 //
467 bool exact;
468 T t = inverse_students_t(df, p, T(1-p), pol, &exact);
469 if((t == 0) || exact)
470 return invert ? -t : t; // can't do better!
471 //
472 // Change variables to inverse incomplete beta:
473 //
474 T t2 = t * t;
475 T xb = df / (df + t2);
476 T y = t2 / (df + t2);
477 T a = df / 2;
478 //
479 // t can be so large that x underflows,
480 // just return our estimate in that case:
481 //
482 if(xb == 0)
483 return t;
484 //
485 // Get incomplete beta and it's derivative:
486 //
487 T f1;
488 T f0 = xb < y ? ibeta_imp(a, constants::half<T>(), xb, pol, false, true, &f1)
489 : ibeta_imp(constants::half<T>(), a, y, pol, true, true, &f1);
490
491 // Get cdf from incomplete beta result:
492 T p0 = f0 / 2 - p;
493 // Get pdf from derivative:
494 T p1 = f1 * sqrt(y * xb * xb * xb / df);
495 //
496 // Second derivative divided by p1:
497 //
498 // yacas gives:
499 //
500 // In> PrettyForm(Simplify(D(t) (1 + t^2/v) ^ (-(v+1)/2)))
501 //
502 // | | v + 1 | |
503 // | -| ----- + 1 | |
504 // | | 2 | |
505 // -| | 2 | |
506 // | | t | |
507 // | | -- + 1 | |
508 // | ( v + 1 ) * | v | * t |
509 // ---------------------------------------------
510 // v
511 //
512 // Which after some manipulation is:
513 //
514 // -p1 * t * (df + 1) / (t^2 + df)
515 //
516 T p2 = t * (df + 1) / (t * t + df);
517 // Halley step:
518 t = fabs(t);
519 t += p0 / (p1 + p0 * p2 / 2);
520 return !invert ? -t : t;
521}
522
523template <class T, class Policy>
524inline T fast_students_t_quantile(T df, T p, const Policy& pol)
525{
526 typedef typename policies::evaluation<T, Policy>::type value_type;
527 typedef typename policies::normalise<
528 Policy,
529 policies::promote_float<false>,
530 policies::promote_double<false>,
531 policies::discrete_quantile<>,
532 policies::assert_undefined<> >::type forwarding_policy;
533
534 typedef boost::integral_constant<bool,
535 (std::numeric_limits<T>::digits <= 53)
536 &&
537 (std::numeric_limits<T>::is_specialized)
538 &&
539 (std::numeric_limits<T>::radix == 2)
540 > tag_type;
541 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%)");
542}
543
544}}} // namespaces
545
546#endif // BOOST_MATH_SF_DETAIL_INV_T_HPP
547
548
549
550

source code of include/boost/math/special_functions/detail/t_distribution_inv.hpp