1 /*
2 * R : A Computer Language for Statistical Data Analysis
3 * Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka
4 * Copyright (C) 1998--2007 The R Core Team
5 * based on code (C) 1979 and later Royal Statistical Society
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
22 * Reference:
23 * Cran, G. W., K. J. Martin and G. E. Thomas (1977).
24 * Remark AS R19 and Algorithm AS 109,
25 * Applied Statistics, 26(1), 111-114.
26 * Remark AS R83 (v.39, 309-310) and the correction (v.40(1) p.236)
27 * have been incorporated in this version.
28 */
29
30
31 #include "nmath.h"
32 #include "dpq.h"
33
34 /* set the exponent of accu to -2r-2 for r digits of accuracy */
35 /*---- NEW ---- -- still fails for p = 1e11, q=.5*/
36
37 #define fpu 3e-308
38 /* acu_min: Minimal value for accuracy 'acu' which will depend on (a,p);
39 acu_min >= fpu ! */
40 #define acu_min 1e-300
41 #define lower fpu
42 #define upper 1-2.22e-16
43
44 #define const1 2.30753
45 #define const2 0.27061
46 #define const3 0.99229
47 #define const4 0.04481
48
49
qbeta(double alpha,double p,double q,int lower_tail,int log_p)50 double qbeta(double alpha, double p, double q, int lower_tail, int log_p)
51 {
52 int swap_tail, i_pb, i_inn;
53 double a, adj, logbeta, g, h, pp, p_, prev, qq, r, s, t, tx, w, y, yprev;
54 double acu;
55 volatile double xinbta;
56
57 /* test for admissibility of parameters */
58
59 #ifdef IEEE_754
60 if (ISNAN(p) || ISNAN(q) || ISNAN(alpha))
61 return p + q + alpha;
62 #endif
63 if(p < 0. || q < 0.) ML_ERR_return_NAN;
64
65 R_Q_P01_boundaries(alpha, 0, 1);
66
67 p_ = R_DT_qIv(alpha);/* lower_tail prob (in any case) */
68
69 if(log_p && (p_ == 0. || p_ == 1.))
70 return p_; /* better than NaN or infinite loop;
71 FIXME: suboptimal, since -Inf < alpha ! */
72
73 /* initialize */
74 logbeta = lbeta(p, q);
75
76 /* change tail if necessary; afterwards 0 < a <= 1/2 */
77 if (p_ <= 0.5) {
78 a = p_; pp = p; qq = q; swap_tail = 0;
79 } else { /* change tail, swap p <-> q :*/
80 a = (!lower_tail && !log_p)? alpha : 1 - p_;
81 pp = q; qq = p; swap_tail = 1;
82 }
83
84 /* calculate the initial approximation */
85
86 /* y := {fast approximation of} qnorm(1 - a) :*/
87 r = sqrt(-2 * log(a));
88 y = r - (const1 + const2 * r) / (1. + (const3 + const4 * r) * r);
89 if (pp > 1 && qq > 1) {
90 r = (y * y - 3.) / 6.;
91 s = 1. / (pp + pp - 1.);
92 t = 1. / (qq + qq - 1.);
93 h = 2. / (s + t);
94 w = y * sqrt(h + r) / h - (t - s) * (r + 5. / 6. - 2. / (3. * h));
95 xinbta = pp / (pp + qq * exp(w + w));
96 } else {
97 r = qq + qq;
98 t = 1. / (9. * qq);
99 t = r * pow(1. - t + y * sqrt(t), 3.0);
100 if (t <= 0.)
101 xinbta = 1. - exp((log1p(-a)+ log(qq) + logbeta) / qq);
102 else {
103 t = (4. * pp + r - 2.) / t;
104 if (t <= 1.)
105 xinbta = exp((log(a * pp) + logbeta) / pp);
106 else
107 xinbta = 1. - 2. / (t + 1.);
108 }
109 }
110
111 /* solve for x by a modified newton-raphson method, */
112 /* using the function pbeta_raw */
113
114 r = 1 - pp;
115 t = 1 - qq;
116 yprev = 0.;
117 adj = 1;
118 /* Sometimes the approximation is negative! */
119 if (xinbta < lower)
120 xinbta = 0.5;
121 else if (xinbta > upper)
122 xinbta = 0.5;
123
124 /* Desired accuracy should depend on (a,p)
125 * This is from Remark .. on AS 109, adapted.
126 * However, it's not clear if this is "optimal" for IEEE double prec.
127
128 * acu = fmax2(acu_min, pow(10., -25. - 5./(pp * pp) - 1./(a * a)));
129
130 * NEW: 'acu' accuracy NOT for squared adjustment, but simple;
131 * ---- i.e., "new acu" = sqrt(old acu)
132
133 */
134 acu = fmax2(acu_min, pow(10., -13 - 2.5/(pp * pp) - 0.5/(a * a)));
135 tx = prev = 0.; /* keep -Wall happy */
136
137 for (i_pb=0; i_pb < 1000; i_pb++) {
138 y = pbeta_raw(xinbta, pp, qq, /*lower_tail = */ TRUE, FALSE);
139 #ifdef IEEE_754
140 if(!R_FINITE(y))
141 #else
142 if (errno)
143 #endif
144 ML_ERR_return_NAN;
145
146 y = (y - a) *
147 exp(logbeta + r * log(xinbta) + t * log1p(-xinbta));
148 if (y * yprev <= 0.)
149 prev = fmax2(fabs(adj),fpu);
150 g = 1;
151 for (i_inn=0; i_inn < 1000;i_inn++) {
152 adj = g * y;
153 if (fabs(adj) < prev) {
154 tx = xinbta - adj; /* trial new x */
155 if (tx >= 0. && tx <= 1) {
156 if (prev <= acu) goto L_converged;
157 if (fabs(y) <= acu) goto L_converged;
158 if (tx != 0. && tx != 1)
159 break;
160 }
161 }
162 g /= 3;
163 }
164 if (fabs(tx - xinbta) < 1e-15*xinbta) goto L_converged;
165 xinbta = tx;
166 yprev = y;
167 }
168 /*-- NOT converged: Iteration count --*/
169 ML_ERROR(ME_PRECISION, "qbeta");
170
171 L_converged:
172 return swap_tail ? 1 - xinbta : xinbta;
173 }
174