1 /*
2 * Mathlib : A C Library of Special Functions
3 * Copyright (C) 1998 Ross Ihaka
4 * Copyright (C) 2000-2013 The R Core Team
5 * Copyright (C) 2003-2013 The R Foundation
6 *
7 * This program is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program; if not, a copy is available at
19 * http://www.r-project.org/Licenses/
20 *
21 * DESCRIPTION
22 *
23 * The "Student" t distribution quantile function.
24 *
25 * NOTES
26 *
27 * This is a C translation of the Fortran routine given in:
28 * Hill, G.W (1970) "Algorithm 396: Student's t-quantiles"
29 * CACM 13(10), 619-620.
30 *
31 * Supplemented by inversion for 0 < ndf < 1.
32 *
33 * ADDITIONS:
34 * - lower_tail, log_p
35 * - using expm1() : takes care of Lozy (1979) "Remark on Algo.", TOMS
36 * - Apply 2-term Taylor expansion as in
37 * Hill, G.W (1981) "Remark on Algo.396", ACM TOMS 7, 250-1
38 * - Improve the formula decision for 1 < df < 2
39 */
40
41 #include "nmath.h"
42 #include "dpq.h"
43
qt(double p,double ndf,int lower_tail,int log_p)44 double qt(double p, double ndf, int lower_tail, int log_p)
45 {
46 const static double eps = 1.e-12;
47
48 double P, q;
49
50 #ifdef IEEE_754
51 if (ISNAN(p) || ISNAN(ndf))
52 return p + ndf;
53 #endif
54
55 R_Q_P01_boundaries(p, ML_NEGINF, ML_POSINF);
56
57 if (ndf <= 0) ML_ERR_return_NAN;
58
59 if (ndf < 1) { /* based on qnt */
60 const static double accu = 1e-13;
61 const static double Eps = 1e-11; /* must be > accu */
62
63 double ux, lx, nx, pp;
64
65 int iter = 0;
66
67 p = R_DT_qIv(p);
68
69 /* Invert pt(.) :
70 * 1. finding an upper and lower bound */
71 if(p > 1 - DBL_EPSILON) return ML_POSINF;
72 pp = fmin2(1 - DBL_EPSILON, p * (1 + Eps));
73 for(ux = 1.; ux < DBL_MAX && pt(ux, ndf, TRUE, FALSE) < pp; ux *= 2);
74 pp = p * (1 - Eps);
75 for(lx =-1.; lx > -DBL_MAX && pt(lx, ndf, TRUE, FALSE) > pp; lx *= 2);
76
77 /* 2. interval (lx,ux) halving
78 regula falsi failed on qt(0.1, 0.1)
79 */
80 do {
81 nx = 0.5 * (lx + ux);
82 if (pt(nx, ndf, TRUE, FALSE) > p) ux = nx; else lx = nx;
83 } while ((ux - lx) / fabs(nx) > accu && ++iter < 1000);
84
85 if(iter >= 1000) ML_ERROR(ME_PRECISION, "qt");
86
87 return 0.5 * (lx + ux);
88 }
89
90 /* Old comment:
91 * FIXME: "This test should depend on ndf AND p !!
92 * ----- and in fact should be replaced by
93 * something like Abramowitz & Stegun 26.7.5 (p.949)"
94 *
95 * That would say that if the qnorm value is x then
96 * the result is about x + (x^3+x)/4df + (5x^5+16x^3+3x)/96df^2
97 * The differences are tiny even if x ~ 1e5, and qnorm is not
98 * that accurate in the extreme tails.
99 */
100 if (ndf > 1e20) return qnorm(p, 0., 1., lower_tail, log_p);
101
102 P = R_D_qIv(p); /* if exp(p) underflows, we fix below */
103
104 Rboolean neg = (!lower_tail || P < 0.5) && (lower_tail || P > 0.5),
105 is_neg_lower = (lower_tail == neg); /* both TRUE or FALSE == !xor */
106 if(neg)
107 P = 2 * (log_p ? (lower_tail ? P : -expm1(p)) : R_D_Lval(p));
108 else
109 P = 2 * (log_p ? (lower_tail ? -expm1(p) : P) : R_D_Cval(p));
110 /* 0 <= P <= 1 ; P = 2*min(P', 1 - P') in all cases */
111
112 if (fabs(ndf - 2) < eps) { /* df ~= 2 */
113 if(P > DBL_MIN) {
114 if(3* P < DBL_EPSILON) /* P ~= 0 */
115 q = 1 / sqrt(P);
116 else if (P > 0.9) /* P ~= 1 */
117 q = (1 - P) * sqrt(2 /(P * (2 - P)));
118 else /* eps/3 <= P <= 0.9 */
119 q = sqrt(2 / (P * (2 - P)) - 2);
120 }
121 else { /* P << 1, q = 1/sqrt(P) = ... */
122 if(log_p)
123 q = is_neg_lower ? exp(- p/2) / M_SQRT2 : 1/sqrt(-expm1(p));
124 else
125 q = ML_POSINF;
126 }
127 }
128 else if (ndf < 1 + eps) { /* df ~= 1 (df < 1 excluded above): Cauchy */
129 if(P > 0)
130 q = 1/tan(P * M_PI_2);/* == - tan((P+1) * M_PI_2) -- suffers for P ~= 0 */
131
132 else { /* P = 0, but maybe = 2*exp(p) ! */
133 if(log_p) /* 1/tan(e) ~ 1/e */
134 q = is_neg_lower ? M_1_PI * exp(-p) : -1./(M_PI * expm1(p));
135 else
136 q = ML_POSINF;
137 }
138 }
139 else { /*-- usual case; including, e.g., df = 1.1 */
140 double x = 0., y, log_P2 = 0./* -Wall */,
141 a = 1 / (ndf - 0.5),
142 b = 48 / (a * a),
143 c = ((20700 * a / b - 98) * a - 16) * a + 96.36,
144 d = ((94.5 / (b + c) - 3) / b + 1) * sqrt(a * M_PI_2) * ndf;
145
146 Rboolean P_ok1 = P > DBL_MIN || !log_p, P_ok = P_ok1;
147 if(P_ok1) {
148 y = pow(d * P, 2 / ndf);
149 P_ok = (y >= DBL_EPSILON);
150 }
151 if(!P_ok) {// log.p && P very.small || (d*P)^(2/df) =: y < eps_c
152 log_P2 = is_neg_lower ? R_D_log(p) : R_D_LExp(p); /* == log(P / 2) */
153 x = (log(d) + M_LN2 + log_P2) / ndf;
154 y = exp(2 * x);
155 }
156
157 if ((ndf < 2.1 && P > 0.5) || y > 0.05 + a) { /* P > P0(df) */
158 /* Asymptotic inverse expansion about normal */
159 if(P_ok)
160 x = qnorm(0.5 * P, 0., 1., /*lower_tail*/TRUE, /*log_p*/FALSE);
161 else /* log_p && P underflowed */
162 x = qnorm(log_P2, 0., 1., lower_tail, /*log_p*/ TRUE);
163
164 y = x * x;
165 if (ndf < 5)
166 c += 0.3 * (ndf - 4.5) * (x + 0.6);
167 c = (((0.05 * d * x - 5) * x - 7) * x - 2) * x + b + c;
168 y = (((((0.4 * y + 6.3) * y + 36) * y + 94.5) / c
169 - y - 3) / b + 1) * x;
170 y = expm1(a * y * y);
171 q = sqrt(ndf * y);
172 } else if(!P_ok && x < - M_LN2 * DBL_MANT_DIG) {/* 0.5* log(DBL_EPSILON) */
173 /* y above might have underflown */
174 q = sqrt(ndf) * exp(-x);
175 }
176 else { /* re-use 'y' from above */
177 y = ((1 / (((ndf + 6) / (ndf * y) - 0.089 * d - 0.822)
178 * (ndf + 2) * 3) + 0.5 / (ndf + 4))
179 * y - 1) * (ndf + 1) / (ndf + 2) + 1 / y;
180 q = sqrt(ndf * y);
181 }
182
183
184 /* Now apply 2-term Taylor expansion improvement (1-term = Newton):
185 * as by Hill (1981) [ref.above] */
186
187 /* FIXME: This can be far from optimal when log_p = TRUE
188 * but is still needed, e.g. for qt(-2, df=1.01, log=TRUE).
189 * Probably also improvable when lower_tail = FALSE */
190
191 if(P_ok1) {
192 int it=0;
193 while(it++ < 10 && (y = dt(q, ndf, FALSE)) > 0 &&
194 R_FINITE(x = (pt(q, ndf, FALSE, FALSE) - P/2) / y) &&
195 fabs(x) > 1e-14*fabs(q))
196 /* Newton (=Taylor 1 term):
197 * q += x;
198 * Taylor 2-term : */
199 q += x * (1. + x * q * (ndf + 1) / (2 * (q * q + ndf)));
200 }
201 }
202 if(neg) q = -q;
203 return q;
204 }
205