1 /*
2 * Mathlib : A C Library of Special Functions
3 * Copyright (C) 1998 Ross Ihaka
4 * Copyright (C) 2000--2011 The R Core Team
5 * Copyright (C) 2004--2009 The R Foundation
6 * based on AS 91 (C) 1979 Royal Statistical Society
7 *
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful, but
14 * WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 * General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program; if not, a copy is available at
20 * http://www.r-project.org/Licenses/
21 *
22 * DESCRIPTION
23 *
24 * Compute the quantile function of the gamma distribution.
25 *
26 * NOTES
27 *
28 * This function is based on the Applied Statistics
29 * Algorithm AS 91 ("ppchi2") and via pgamma(.) AS 239.
30 *
31 * R core improvements:
32 * o lower_tail, log_p
33 * o non-trivial result for p outside [0.000002, 0.999998]
34 * o p ~ 1 no longer gives +Inf; final Newton step(s)
35 *
36 * REFERENCES
37 *
38 * Best, D. J. and D. E. Roberts (1975).
39 * Percentage Points of the Chi-Squared Distribution.
40 * Applied Statistics 24, page 385. */
41
42 #include "nmath.h"
43 #include "dpq.h"
44
45 #ifdef DEBUG_qgamma
46 # define DEBUG_q
47 #endif
48
49 attribute_hidden
qchisq_appr(double p,double nu,double g,int lower_tail,int log_p,double tol)50 double qchisq_appr(double p, double nu, double g /* = log Gamma(nu/2) */,
51 int lower_tail, int log_p, double tol /* EPS1 */)
52 {
53 #define C7 4.67
54 #define C8 6.66
55 #define C9 6.73
56 #define C10 13.32
57
58 double alpha, a, c, ch, p1;
59 double p2, q, t, x;
60
61 /* test arguments and initialise */
62
63 #ifdef IEEE_754
64 if (ISNAN(p) || ISNAN(nu))
65 return p + nu;
66 #endif
67 R_Q_P01_check(p);
68 if (nu <= 0) ML_ERR_return_NAN;
69
70 alpha = 0.5 * nu;/* = [pq]gamma() shape */
71 c = alpha-1;
72
73 if(nu < (-1.24)*(p1 = R_DT_log(p))) { /* for small chi-squared */
74 /* log(alpha) + g = log(alpha) + log(gamma(alpha)) =
75 * = log(alpha*gamma(alpha)) = lgamma(alpha+1) suffers from
76 * catastrophic cancellation when alpha << 1
77 */
78 double lgam1pa = (alpha < 0.5) ? lgamma1p(alpha) : (log(alpha) + g);
79 ch = exp((lgam1pa + p1)/alpha + M_LN2);
80 #ifdef DEBUG_qgamma
81 REprintf(" small chi-sq., ch0 = %g\n", ch);
82 #endif
83
84 } else if(nu > 0.32) { /* using Wilson and Hilferty estimate */
85
86 x = qnorm(p, 0, 1, lower_tail, log_p);
87 p1 = 2./(9*nu);
88 ch = nu*pow(x*sqrt(p1) + 1-p1, 3);
89
90 #ifdef DEBUG_qgamma
91 REprintf(" nu > .32: Wilson-Hilferty; x = %7g\n", x);
92 #endif
93 /* approximation for p tending to 1: */
94 if( ch > 2.2*nu + 6 )
95 ch = -2*(R_DT_Clog(p) - c*log(0.5*ch) + g);
96
97 } else { /* "small nu" : 1.24*(-log(p)) <= nu <= 0.32 */
98
99 ch = 0.4;
100 a = R_DT_Clog(p) + g + c*M_LN2;
101 #ifdef DEBUG_qgamma
102 REprintf(" nu <= .32: a = %7g\n", a);
103 #endif
104 do {
105 q = ch;
106 p1 = 1. / (1+ch*(C7+ch));
107 p2 = ch*(C9+ch*(C8+ch));
108 t = -0.5 +(C7+2*ch)*p1 - (C9+ch*(C10+3*ch))/p2;
109 ch -= (1- exp(a+0.5*ch)*p2*p1)/t;
110 } while(fabs(q - ch) > tol * fabs(ch));
111 }
112
113 return ch;
114 }
115
qgamma(double p,double alpha,double scale,int lower_tail,int log_p)116 double qgamma(double p, double alpha, double scale, int lower_tail, int log_p)
117 /* shape = alpha */
118 {
119 #define EPS1 1e-2
120 #define EPS2 5e-7/* final precision of AS 91 */
121 #define EPS_N 1e-15/* precision of Newton step / iterations */
122 #define LN_EPS -36.043653389117156 /* = log(.Machine$double.eps) iff IEEE_754 */
123
124 #define MAXIT 1000/* was 20 */
125
126 #define pMIN 1e-100 /* was 0.000002 = 2e-6 */
127 #define pMAX (1-1e-14)/* was (1-1e-12) and 0.999998 = 1 - 2e-6 */
128
129 const static double
130 i420 = 1./ 420.,
131 i2520 = 1./ 2520.,
132 i5040 = 1./ 5040;
133
134 double p_, a, b, c, g, ch, ch0, p1;
135 double p2, q, s1, s2, s3, s4, s5, s6, t, x;
136 int i, max_it_Newton = 1;
137
138 /* test arguments and initialise */
139
140 #ifdef IEEE_754
141 if (ISNAN(p) || ISNAN(alpha) || ISNAN(scale))
142 return p + alpha + scale;
143 #endif
144 R_Q_P01_boundaries(p, 0., ML_POSINF);
145
146 if (alpha < 0 || scale <= 0) ML_ERR_return_NAN;
147
148 if (alpha == 0) /* all mass at 0 : */ return 0.;
149
150 if (alpha < 1e-10) {
151 /* Warning seems unnecessary now: */
152 #ifdef _DO_WARN_qgamma_
153 MATHLIB_WARNING("value of shape (%g) is extremely small: results may be unreliable", alpha);
154 #endif
155 max_it_Newton = 7;/* may still be increased below */
156 }
157
158 p_ = R_DT_qIv(p);/* lower_tail prob (in any case) */
159
160 #ifdef DEBUG_qgamma
161 REprintf("qgamma(p=%7g, alpha=%7g, scale=%7g, l.t.=%2d, log_p=%2d): ",
162 p,alpha,scale, lower_tail, log_p);
163 #endif
164 g = lgammafn(alpha);/* log Gamma(v/2) */
165
166 /*----- Phase I : Starting Approximation */
167 ch = qchisq_appr(p, /* nu= 'df' = */ 2*alpha, /* lgamma(nu/2)= */ g,
168 lower_tail, log_p, /* tol= */ EPS1);
169 if(!R_FINITE(ch)) {
170 /* forget about all iterations! */
171 max_it_Newton = 0; goto END;
172 }
173 if(ch < EPS2) {/* Corrected according to AS 91; MM, May 25, 1999 */
174 max_it_Newton = 20;
175 goto END;/* and do Newton steps */
176 }
177
178 /* FIXME: This (cutoff to {0, +Inf}) is far from optimal
179 * ----- when log_p or !lower_tail, but NOT doing it can be even worse */
180 if(p_ > pMAX || p_ < pMIN) {
181 /* did return ML_POSINF or 0.; much better: */
182 max_it_Newton = 20;
183 goto END;/* and do Newton steps */
184 }
185
186 #ifdef DEBUG_qgamma
187 REprintf("\t==> ch = %10g:", ch);
188 #endif
189
190 /*----- Phase II: Iteration
191 * Call pgamma() [AS 239] and calculate seven term taylor series
192 */
193 c = alpha-1;
194 s6 = (120+c*(346+127*c)) * i5040; /* used below, is "const" */
195
196 ch0 = ch;/* save initial approx. */
197 for(i=1; i <= MAXIT; i++ ) {
198 q = ch;
199 p1 = 0.5*ch;
200 p2 = p_ - pgamma_raw(p1, alpha, /*lower_tail*/TRUE, /*log_p*/FALSE);
201 #ifdef DEBUG_qgamma
202 if(i == 1) REprintf(" Ph.II iter; ch=%g, p2=%g\n", ch, p2);
203 if(i >= 2) REprintf(" it=%d, ch=%g, p2=%g\n", i, ch, p2);
204 #endif
205 #ifdef IEEE_754
206 if(!R_FINITE(p2) || ch <= 0)
207 #else
208 if(errno != 0 || ch <= 0)
209 #endif
210 { ch = ch0; max_it_Newton = 27; goto END; }/*was return ML_NAN;*/
211
212 t = p2*exp(alpha*M_LN2+g+p1-c*log(ch));
213 b = t/ch;
214 a = 0.5*t - b*c;
215 s1 = (210+ a*(140+a*(105+a*(84+a*(70+60*a))))) * i420;
216 s2 = (420+ a*(735+a*(966+a*(1141+1278*a)))) * i2520;
217 s3 = (210+ a*(462+a*(707+932*a))) * i2520;
218 s4 = (252+ a*(672+1182*a) + c*(294+a*(889+1740*a))) * i5040;
219 s5 = (84+2264*a + c*(1175+606*a)) * i2520;
220
221 ch += t*(1+0.5*t*s1-b*c*(s1-b*(s2-b*(s3-b*(s4-b*(s5-b*s6))))));
222 if(fabs(q - ch) < EPS2*ch)
223 goto END;
224 if(fabs(q - ch) > 0.1*ch) {/* diverging? -- also forces ch > 0 */
225 if(ch < q) ch = 0.9 * q; else ch = 1.1 * q;
226 }
227 }
228 /* no convergence in MAXIT iterations -- but we add Newton now... */
229 #ifdef DEBUG_q
230 MATHLIB_WARNING3("qgamma(%g) not converged in %d iterations; rel.ch=%g\n",
231 p, MAXIT, ch/fabs(q - ch));
232 #endif
233 /* was
234 * ML_ERROR(ME_PRECISION, "qgamma");
235 * does nothing in R !*/
236
237 END:
238 /* PR# 2214 : From: Morten Welinder <terra@diku.dk>, Fri, 25 Oct 2002 16:50
239 -------- To: R-bugs@biostat.ku.dk Subject: qgamma precision
240
241 * With a final Newton step, double accuracy, e.g. for (p= 7e-4; nu= 0.9)
242 *
243 * Improved (MM): - only if rel.Err > EPS_N (= 1e-15);
244 * - also for lower_tail = FALSE or log_p = TRUE
245 * - optionally *iterate* Newton
246 */
247 x = 0.5*scale*ch;
248 if(max_it_Newton) {
249 /* always use log scale */
250 if (!log_p) {
251 p = log(p);
252 log_p = TRUE;
253 }
254 if(x == 0) {
255 const double _1_p = 1. + 1e-7;
256 const double _1_m = 1. - 1e-7;
257 x = DBL_MIN;
258 p_ = pgamma(x, alpha, scale, lower_tail, log_p);
259 if(( lower_tail && p_ > p * _1_p) ||
260 (!lower_tail && p_ < p * _1_m))
261 return(0.);
262 /* else: continue, using x = DBL_MIN instead of 0 */
263 }
264 else
265 p_ = pgamma(x, alpha, scale, lower_tail, log_p);
266 if(p_ == ML_NEGINF) return 0; /* PR#14710 */
267 for(i = 1; i <= max_it_Newton; i++) {
268 p1 = p_ - p;
269 #ifdef DEBUG_qgamma
270 if(i == 1) REprintf("\n it=%d: p=%g, x = %g, p.=%g; p1=d{p}=%g\n",
271 i, p, x, p_, p1);
272 if(i >= 2) REprintf(" x{it= %d} = %g, p.=%g, p1=d{p}=%g\n",
273 i, x, p_, p1);
274 #endif
275 if(fabs(p1) < fabs(EPS_N * p))
276 break;
277 /* else */
278 if((g = dgamma(x, alpha, scale, log_p)) == R_D__0) {
279 #ifdef DEBUG_q
280 if(i == 1) REprintf("no final Newton step because dgamma(*)== 0!\n");
281 #endif
282 break;
283 }
284 /* else :
285 * delta x = f(x)/f'(x);
286 * if(log_p) f(x) := log P(x) - p; f'(x) = d/dx log P(x) = P' / P
287 * ==> f(x)/f'(x) = f*P / P' = f*exp(p_) / P' (since p_ = log P(x))
288 */
289 t = log_p ? p1*exp(p_ - g) : p1/g ;/* = "delta x" */
290 t = lower_tail ? x - t : x + t;
291 p_ = pgamma (t, alpha, scale, lower_tail, log_p);
292 if (fabs(p_ - p) > fabs(p1) ||
293 (i > 1 && fabs(p_ - p) == fabs(p1)) /* <- against flip-flop */) {
294 /* no improvement */
295 #ifdef DEBUG_qgamma
296 if(i == 1 && max_it_Newton > 1)
297 REprintf("no Newton step done since delta{p} >= last delta\n");
298 #endif
299 break;
300 } /* else : */
301 #ifdef Harmful_notably_if_max_it_Newton_is_1
302 /* control step length: this could have started at
303 the initial approximation */
304 if(t > 1.1*x) t = 1.1*x;
305 else if(t < 0.9*x) t = 0.9*x;
306 #endif
307 x = t;
308 }
309 }
310
311 return x;
312 }
313