1 /* incbet.c
2 *
3 * Incomplete beta integral
4 *
5 *
6 * SYNOPSIS:
7 *
8 * double a, b, x, y, incbet();
9 *
10 * y = incbet(a, b, x);
11 *
12 *
13 * DESCRIPTION:
14 *
15 * Returns incomplete beta integral of the arguments, evaluated
16 * from zero to x. The function is defined as
17 *
18 * x
19 * - -
20 * | (a+b) | | a-1 b-1
21 * ----------- | t (1-t) dt.
22 * - - | |
23 * | (a) | (b) -
24 * 0
25 *
26 * The domain of definition is 0 <= x <= 1. In this
27 * implementation a and b are restricted to positive values.
28 * The integral from x to 1 may be obtained by the symmetry
29 * relation
30 *
31 * 1 - incbet(a, b, x) = incbet(b, a, 1-x).
32 *
33 * The integral is evaluated by a continued fraction expansion
34 * or, when b*x is small, by a power series.
35 *
36 * ACCURACY:
37 *
38 * Tested at uniformly distributed random points (a,b,x) with a and b
39 * in "domain" and x between 0 and 1.
40 * Relative error
41 * arithmetic domain # trials peak rms
42 * IEEE 0,5 10000 6.9e-15 4.5e-16
43 * IEEE 0,85 250000 2.2e-13 1.7e-14
44 * IEEE 0,1000 30000 5.3e-12 6.3e-13
45 * IEEE 0,10000 250000 9.3e-11 7.1e-12
46 * IEEE 0,100000 10000 8.7e-10 4.8e-11
47 * Outputs smaller than the IEEE gradual underflow threshold
48 * were excluded from these statistics.
49 *
50 * ERROR MESSAGES:
51 * message condition value returned
52 * incbet domain x<0, x>1 0.0
53 * incbet underflow 0.0
54 */
55
56
57 /*
58 Cephes Math Library, Release 2.8: June, 2000
59 Copyright 1984, 1995, 2000 by Stephen L. Moshier
60 */
61
62 #include "mconf.h"
63
64 #ifdef DEC
65 #define MAXGAM 34.84425627277176174
66 #else
67 #define MAXGAM 171.624376956302725
68 #endif
69
70 static double incbcf (double, double, double);
71 static double incbd (double, double, double);
72 static double pseries (double, double, double);
73
74 static double big = 4.503599627370496e15;
75 static double biginv = 2.22044604925031308085e-16;
76
77
incbet(double aa,double bb,double xx)78 double incbet (double aa, double bb, double xx)
79 {
80 double a, b, t, x, xc, w, y;
81 int flag;
82
83 if (aa <= 0.0 || bb <= 0.0) {
84 goto domerr;
85 }
86
87 if (xx <= 0.0 || xx >= 1.0) {
88 if (xx == 0.0) {
89 return 0.0;
90 }
91 if (xx == 1.0) {
92 return 1.0;
93 }
94 domerr:
95 mtherr("incbet", CEPHES_DOMAIN);
96 return 0.0;
97 }
98
99 flag = 0;
100
101 if (bb * xx <= 1.0 && xx <= 0.95) {
102 t = pseries(aa, bb, xx);
103 goto done;
104 }
105
106 w = 1.0 - xx;
107
108 /* Reverse a and b if x is greater than the mean. */
109 if (xx > (aa/(aa+bb))) {
110 flag = 1;
111 a = bb;
112 b = aa;
113 xc = xx;
114 x = w;
115 } else {
116 a = aa;
117 b = bb;
118 xc = w;
119 x = xx;
120 }
121
122 if (flag == 1 && (b * x) <= 1.0 && x <= 0.95) {
123 t = pseries(a, b, x);
124 goto done;
125 }
126
127 /* Choose expansion for better convergence. */
128 y = x * (a+b-2.0) - (a-1.0);
129
130 if (y < 0.0) {
131 w = incbcf(a, b, x);
132 } else {
133 w = incbd(a, b, x) / xc;
134 }
135
136 /* Multiply w by the factor
137 a b _ _ _
138 x (1-x) | (a+b) / (a | (a) | (b)) .
139 */
140
141 y = a * log(x);
142 t = b * log(xc);
143
144 if ((a+b) < MAXGAM && fabs(y) < MAXLOG && fabs(t) < MAXLOG) {
145 t = pow(xc,b);
146 t *= pow(x,a);
147 t /= a;
148 t *= w;
149 t *= cephes_gamma(a+b) / (cephes_gamma(a) * cephes_gamma(b));
150 goto done;
151 }
152
153 /* Resort to logarithms. */
154 y += t + lgam(a+b) - lgam(a) - lgam(b);
155 y += log(w/a);
156
157 if (y < MINLOG) {
158 t = 0.0;
159 } else {
160 t = exp(y);
161 }
162
163 done:
164
165 if (flag == 1) {
166 if (t <= MACHEP) {
167 t = 1.0 - MACHEP;
168 } else {
169 t = 1.0 - t;
170 }
171 }
172
173 return t;
174 }
175
176 /* Continued fraction expansion #1 for incomplete beta integral
177 */
178
incbcf(double a,double b,double x)179 static double incbcf (double a, double b, double x)
180 {
181 double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
182 double k1, k2, k3, k4, k5, k6, k7, k8;
183 double r, t, ans, thresh;
184 int n;
185
186 k1 = a;
187 k2 = a + b;
188 k3 = a;
189 k4 = a + 1.0;
190 k5 = 1.0;
191 k6 = b - 1.0;
192 k7 = k4;
193 k8 = a + 2.0;
194
195 pkm2 = 0.0;
196 qkm2 = 1.0;
197 pkm1 = 1.0;
198 qkm1 = 1.0;
199 ans = 1.0;
200 r = 1.0;
201 n = 0;
202 thresh = 3.0 * MACHEP;
203
204 do {
205 xk = -(x * k1 * k2)/(k3 * k4);
206 pk = pkm1 + pkm2 * xk;
207 qk = qkm1 + qkm2 * xk;
208 pkm2 = pkm1;
209 pkm1 = pk;
210 qkm2 = qkm1;
211 qkm1 = qk;
212
213 xk = (x * k5 * k6)/(k7 * k8);
214 pk = pkm1 + pkm2 * xk;
215 qk = qkm1 + qkm2 * xk;
216 pkm2 = pkm1;
217 pkm1 = pk;
218 qkm2 = qkm1;
219 qkm1 = qk;
220
221 if (qk != 0) {
222 r = pk/qk;
223 }
224
225 if (r != 0) {
226 t = fabs((ans - r)/r);
227 ans = r;
228 } else {
229 t = 1.0;
230 }
231
232 if (t < thresh) {
233 goto cdone;
234 }
235
236 k1 += 1.0;
237 k2 += 1.0;
238 k3 += 2.0;
239 k4 += 2.0;
240 k5 += 1.0;
241 k6 -= 1.0;
242 k7 += 2.0;
243 k8 += 2.0;
244
245 if ((fabs(qk) + fabs(pk)) > big) {
246 pkm2 *= biginv;
247 pkm1 *= biginv;
248 qkm2 *= biginv;
249 qkm1 *= biginv;
250 }
251
252 if ((fabs(qk) < biginv) || (fabs(pk) < biginv)) {
253 pkm2 *= big;
254 pkm1 *= big;
255 qkm2 *= big;
256 qkm1 *= big;
257 }
258 } while (++n < 300);
259
260 cdone:
261
262 return ans;
263 }
264
265 /* Continued fraction expansion #2 for incomplete beta integral
266 */
267
incbd(double a,double b,double x)268 static double incbd (double a, double b, double x)
269 {
270 double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
271 double k1, k2, k3, k4, k5, k6, k7, k8;
272 double r, t, ans, z, thresh;
273 int n;
274
275 k1 = a;
276 k2 = b - 1.0;
277 k3 = a;
278 k4 = a + 1.0;
279 k5 = 1.0;
280 k6 = a + b;
281 k7 = a + 1.0;;
282 k8 = a + 2.0;
283
284 pkm2 = 0.0;
285 qkm2 = 1.0;
286 pkm1 = 1.0;
287 qkm1 = 1.0;
288 z = x / (1.0-x);
289 ans = 1.0;
290 r = 1.0;
291 n = 0;
292 thresh = 3.0 * MACHEP;
293
294 do {
295 xk = -(z * k1 * k2)/(k3 * k4);
296 pk = pkm1 + pkm2 * xk;
297 qk = qkm1 + qkm2 * xk;
298 pkm2 = pkm1;
299 pkm1 = pk;
300 qkm2 = qkm1;
301 qkm1 = qk;
302
303 xk = (z * k5 * k6)/(k7 * k8);
304 pk = pkm1 + pkm2 * xk;
305 qk = qkm1 + qkm2 * xk;
306 pkm2 = pkm1;
307 pkm1 = pk;
308 qkm2 = qkm1;
309 qkm1 = qk;
310
311 if (qk != 0) {
312 r = pk/qk;
313 }
314
315 if (r != 0) {
316 t = fabs((ans - r)/r);
317 ans = r;
318 } else {
319 t = 1.0;
320 }
321
322 if (t < thresh) {
323 goto cdone;
324 }
325
326 k1 += 1.0;
327 k2 -= 1.0;
328 k3 += 2.0;
329 k4 += 2.0;
330 k5 += 1.0;
331 k6 += 1.0;
332 k7 += 2.0;
333 k8 += 2.0;
334
335 if ((fabs(qk) + fabs(pk)) > big) {
336 pkm2 *= biginv;
337 pkm1 *= biginv;
338 qkm2 *= biginv;
339 qkm1 *= biginv;
340 }
341
342 if ((fabs(qk) < biginv) || (fabs(pk) < biginv)) {
343 pkm2 *= big;
344 pkm1 *= big;
345 qkm2 *= big;
346 qkm1 *= big;
347 }
348 } while (++n < 300);
349
350 cdone:
351
352 return ans;
353 }
354
355 /* Power series for incomplete beta integral.
356 Use when b*x is small and x not too close to 1. */
357
pseries(double a,double b,double x)358 static double pseries (double a, double b, double x)
359 {
360 double s, t, u, v, n, t1, z, ai;
361
362 ai = 1.0 / a;
363 u = (1.0 - b) * x;
364 v = u / (a + 1.0);
365 t1 = v;
366 t = u;
367 n = 2.0;
368 s = 0.0;
369 z = MACHEP * ai;
370
371 while (fabs(v) > z) {
372 u = (n - b) * x / n;
373 t *= u;
374 v = t / (a + n);
375 s += v;
376 n += 1.0;
377 }
378
379 s += t1;
380 s += ai;
381
382 u = a * log(x);
383
384 if ((a+b) < MAXGAM && fabs(u) < MAXLOG) {
385 t = cephes_gamma(a+b) / (cephes_gamma(a) * cephes_gamma(b));
386 s = s * t * pow(x, a);
387 } else {
388 t = lgam(a+b) - lgam(a) - lgam(b) + u + log(s);
389 if (t < MINLOG) {
390 s = 0.0;
391 } else {
392 s = exp(t);
393 }
394 }
395
396 return s;
397 }
398