1 /* cdf/betainv.c
2 *
3 * Copyright (C) 2004 Free Software Foundation, Inc.
4 * Copyright (C) 2006, 2007 Brian Gough
5 * Written by Jason H. Stover.
6 * Modified for GSL by Brian Gough
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 3 of the License, or (at
11 * 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, write to the Free Software
20 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.
21 */
22
23 /*
24 * Invert the Beta distribution.
25 *
26 * References:
27 *
28 * Roger W. Abernathy and Robert P. Smith. "Applying Series Expansion
29 * to the Inverse Beta Distribution to Find Percentiles of the
30 * F-Distribution," ACM Transactions on Mathematical Software, volume
31 * 19, number 4, December 1993, pages 474-480.
32 *
33 * G.W. Hill and A.W. Davis. "Generalized asymptotic expansions of a
34 * Cornish-Fisher type," Annals of Mathematical Statistics, volume 39,
35 * number 8, August 1968, pages 1264-1273.
36 */
37
38 #include "gsl__config.h"
39 #include <math.h>
40 #include "gsl_math.h"
41 #include "gsl_errno.h"
42 #include "gsl_sf_gamma.h"
43 #include "gsl_cdf.h"
44 #include "gsl_randist.h"
45
46 #include "gsl_cdf__error.h"
47
48 static double
bisect(double x,double P,double a,double b,double xtol,double Ptol)49 bisect (double x, double P, double a, double b, double xtol, double Ptol)
50 {
51 double x0 = 0, x1 = 1, Px;
52
53 while (fabs(x1 - x0) > xtol) {
54 Px = gsl_cdf_beta_P (x, a, b);
55 if (fabs(Px - P) < Ptol) {
56 /* return as soon as approximation is good enough, including on
57 the first iteration */
58 return x;
59 } else if (Px < P) {
60 x0 = x;
61 } else if (Px > P) {
62 x1 = x;
63 }
64 x = 0.5 * (x0 + x1);
65 }
66 return x;
67 }
68
69
70 double
gsl_cdf_beta_Pinv(const double P,const double a,const double b)71 gsl_cdf_beta_Pinv (const double P, const double a, const double b)
72 {
73 double x, mean;
74
75 if (P < 0.0 || P > 1.0)
76 {
77 CDF_ERROR ("P must be in range 0 < P < 1", GSL_EDOM);
78 }
79
80 if (a < 0.0)
81 {
82 CDF_ERROR ("a < 0", GSL_EDOM);
83 }
84
85 if (b < 0.0)
86 {
87 CDF_ERROR ("b < 0", GSL_EDOM);
88 }
89
90 if (P == 0.0)
91 {
92 return 0.0;
93 }
94
95 if (P == 1.0)
96 {
97 return 1.0;
98 }
99
100 if (P > 0.5)
101 {
102 return gsl_cdf_beta_Qinv (1 - P, a, b);
103 }
104
105 mean = a / (a + b);
106
107 if (P < 0.1)
108 {
109 /* small x */
110
111 double lg_ab = gsl_sf_lngamma (a + b);
112 double lg_a = gsl_sf_lngamma (a);
113 double lg_b = gsl_sf_lngamma (b);
114
115 double lx = (log (a) + lg_a + lg_b - lg_ab + log (P)) / a;
116 if (lx <= 0) {
117 x = exp (lx); /* first approximation */
118 x *= pow (1 - x, -(b - 1) / a); /* second approximation */
119 } else {
120 x = mean;
121 }
122
123 if (x > mean)
124 x = mean;
125 }
126 else
127 {
128 /* Use expected value as first guess */
129 x = mean;
130 }
131
132 /* Do bisection to get closer */
133 x = bisect (x, P, a, b, 0.01, 0.01);
134
135 {
136 double lambda, dP, phi;
137 unsigned int n = 0;
138
139 start:
140 dP = P - gsl_cdf_beta_P (x, a, b);
141 phi = gsl_ran_beta_pdf (x, a, b);
142
143 if (dP == 0.0 || n++ > 64)
144 goto end;
145
146 lambda = dP / GSL_MAX (2 * fabs (dP / x), phi);
147
148 {
149 double step0 = lambda;
150 double step1 = -((a - 1) / x - (b - 1) / (1 - x)) * lambda * lambda / 2;
151
152 double step = step0;
153
154 if (fabs (step1) < fabs (step0))
155 {
156 step += step1;
157 }
158 else
159 {
160 /* scale back step to a reasonable size when too large */
161 step *= 2 * fabs (step0 / step1);
162 };
163
164 if (x + step > 0 && x + step < 1)
165 {
166 x += step;
167 }
168 else
169 {
170 x = sqrt (x) * sqrt (mean); /* try a new starting point */
171 }
172
173 if (fabs (step0) > 1e-10 * x)
174 goto start;
175 }
176
177 end:
178
179 if (fabs(dP) > GSL_SQRT_DBL_EPSILON * P)
180 {
181 GSL_ERROR_VAL("inverse failed to converge", GSL_EFAILED, GSL_NAN);
182 }
183
184 return x;
185 }
186 }
187
188 double
gsl_cdf_beta_Qinv(const double Q,const double a,const double b)189 gsl_cdf_beta_Qinv (const double Q, const double a, const double b)
190 {
191
192 if (Q < 0.0 || Q > 1.0)
193 {
194 CDF_ERROR ("Q must be inside range 0 < Q < 1", GSL_EDOM);
195 }
196
197 if (a < 0.0)
198 {
199 CDF_ERROR ("a < 0", GSL_EDOM);
200 }
201
202 if (b < 0.0)
203 {
204 CDF_ERROR ("b < 0", GSL_EDOM);
205 }
206
207 if (Q == 0.0)
208 {
209 return 1.0;
210 }
211
212 if (Q == 1.0)
213 {
214 return 0.0;
215 }
216
217 if (Q > 0.5)
218 {
219 return gsl_cdf_beta_Pinv (1 - Q, a, b);
220 }
221 else
222 {
223 return 1 - gsl_cdf_beta_Pinv (Q, b, a);
224 };
225 }
226