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.3: March, 1995
59 * Copyright 1984, 1995 by Stephen L. Moshier
60 */
61
62 #include "mconf.h"
63
64 #define MAXGAM 171.624376956302725
65
66 extern double MACHEP, MINLOG, MAXLOG;
67
68 static double big = 4.503599627370496e15;
69 static double biginv = 2.22044604925031308085e-16;
70
71 static double incbcf(double a, double b, double x);
72 static double incbd(double a, double b, double x);
73 static double pseries(double a, double b, double x);
74
incbet(aa,bb,xx)75 double incbet(aa, bb, xx)
76 double aa, bb, xx;
77 {
78 double a, b, t, x, xc, w, y;
79 int flag;
80
81 if (aa <= 0.0 || bb <= 0.0)
82 goto domerr;
83
84 if ((xx <= 0.0) || (xx >= 1.0)) {
85 if (xx == 0.0)
86 return (0.0);
87 if (xx == 1.0)
88 return (1.0);
89 domerr:
90 sf_error("incbet", SF_ERROR_DOMAIN, NULL);
91 return (NPY_NAN);
92 }
93
94 flag = 0;
95 if ((bb * xx) <= 1.0 && xx <= 0.95) {
96 t = pseries(aa, bb, xx);
97 goto done;
98 }
99
100 w = 1.0 - xx;
101
102 /* Reverse a and b if x is greater than the mean. */
103 if (xx > (aa / (aa + bb))) {
104 flag = 1;
105 a = bb;
106 b = aa;
107 xc = xx;
108 x = w;
109 }
110 else {
111 a = aa;
112 b = bb;
113 xc = w;
114 x = xx;
115 }
116
117 if (flag == 1 && (b * x) <= 1.0 && x <= 0.95) {
118 t = pseries(a, b, x);
119 goto done;
120 }
121
122 /* Choose expansion for better convergence. */
123 y = x * (a + b - 2.0) - (a - 1.0);
124 if (y < 0.0)
125 w = incbcf(a, b, x);
126 else
127 w = incbd(a, b, x) / xc;
128
129 /* Multiply w by the factor
130 * a b _ _ _
131 * x (1-x) | (a+b) / ( a | (a) | (b) ) . */
132
133 y = a * log(x);
134 t = b * log(xc);
135 if ((a + b) < MAXGAM && fabs(y) < MAXLOG && fabs(t) < MAXLOG) {
136 t = pow(xc, b);
137 t *= pow(x, a);
138 t /= a;
139 t *= w;
140 t *= 1.0 / beta(a, b);
141 goto done;
142 }
143 /* Resort to logarithms. */
144 y += t - lbeta(a,b);
145 y += log(w / a);
146 if (y < MINLOG)
147 t = 0.0;
148 else
149 t = exp(y);
150
151 done:
152
153 if (flag == 1) {
154 if (t <= MACHEP)
155 t = 1.0 - MACHEP;
156 else
157 t = 1.0 - t;
158 }
159 return (t);
160 }
161
162 /* Continued fraction expansion #1
163 * for incomplete beta integral
164 */
165
incbcf(a,b,x)166 static double incbcf(a, b, x)
167 double a, b, x;
168 {
169 double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
170 double k1, k2, k3, k4, k5, k6, k7, k8;
171 double r, t, ans, thresh;
172 int n;
173
174 k1 = a;
175 k2 = a + b;
176 k3 = a;
177 k4 = a + 1.0;
178 k5 = 1.0;
179 k6 = b - 1.0;
180 k7 = k4;
181 k8 = a + 2.0;
182
183 pkm2 = 0.0;
184 qkm2 = 1.0;
185 pkm1 = 1.0;
186 qkm1 = 1.0;
187 ans = 1.0;
188 r = 1.0;
189 n = 0;
190 thresh = 3.0 * MACHEP;
191 do {
192
193 xk = -(x * k1 * k2) / (k3 * k4);
194 pk = pkm1 + pkm2 * xk;
195 qk = qkm1 + qkm2 * xk;
196 pkm2 = pkm1;
197 pkm1 = pk;
198 qkm2 = qkm1;
199 qkm1 = qk;
200
201 xk = (x * k5 * k6) / (k7 * k8);
202 pk = pkm1 + pkm2 * xk;
203 qk = qkm1 + qkm2 * xk;
204 pkm2 = pkm1;
205 pkm1 = pk;
206 qkm2 = qkm1;
207 qkm1 = qk;
208
209 if (qk != 0)
210 r = pk / qk;
211 if (r != 0) {
212 t = fabs((ans - r) / r);
213 ans = r;
214 }
215 else
216 t = 1.0;
217
218 if (t < thresh)
219 goto cdone;
220
221 k1 += 1.0;
222 k2 += 1.0;
223 k3 += 2.0;
224 k4 += 2.0;
225 k5 += 1.0;
226 k6 -= 1.0;
227 k7 += 2.0;
228 k8 += 2.0;
229
230 if ((fabs(qk) + fabs(pk)) > big) {
231 pkm2 *= biginv;
232 pkm1 *= biginv;
233 qkm2 *= biginv;
234 qkm1 *= biginv;
235 }
236 if ((fabs(qk) < biginv) || (fabs(pk) < biginv)) {
237 pkm2 *= big;
238 pkm1 *= big;
239 qkm2 *= big;
240 qkm1 *= big;
241 }
242 }
243 while (++n < 300);
244
245 cdone:
246 return (ans);
247 }
248
249
250 /* Continued fraction expansion #2
251 * for incomplete beta integral
252 */
253
incbd(a,b,x)254 static double incbd(a, b, x)
255 double a, b, x;
256 {
257 double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
258 double k1, k2, k3, k4, k5, k6, k7, k8;
259 double r, t, ans, z, thresh;
260 int n;
261
262 k1 = a;
263 k2 = b - 1.0;
264 k3 = a;
265 k4 = a + 1.0;
266 k5 = 1.0;
267 k6 = a + b;
268 k7 = a + 1.0;;
269 k8 = a + 2.0;
270
271 pkm2 = 0.0;
272 qkm2 = 1.0;
273 pkm1 = 1.0;
274 qkm1 = 1.0;
275 z = x / (1.0 - x);
276 ans = 1.0;
277 r = 1.0;
278 n = 0;
279 thresh = 3.0 * MACHEP;
280 do {
281
282 xk = -(z * k1 * k2) / (k3 * k4);
283 pk = pkm1 + pkm2 * xk;
284 qk = qkm1 + qkm2 * xk;
285 pkm2 = pkm1;
286 pkm1 = pk;
287 qkm2 = qkm1;
288 qkm1 = qk;
289
290 xk = (z * k5 * k6) / (k7 * k8);
291 pk = pkm1 + pkm2 * xk;
292 qk = qkm1 + qkm2 * xk;
293 pkm2 = pkm1;
294 pkm1 = pk;
295 qkm2 = qkm1;
296 qkm1 = qk;
297
298 if (qk != 0)
299 r = pk / qk;
300 if (r != 0) {
301 t = fabs((ans - r) / r);
302 ans = r;
303 }
304 else
305 t = 1.0;
306
307 if (t < thresh)
308 goto cdone;
309
310 k1 += 1.0;
311 k2 -= 1.0;
312 k3 += 2.0;
313 k4 += 2.0;
314 k5 += 1.0;
315 k6 += 1.0;
316 k7 += 2.0;
317 k8 += 2.0;
318
319 if ((fabs(qk) + fabs(pk)) > big) {
320 pkm2 *= biginv;
321 pkm1 *= biginv;
322 qkm2 *= biginv;
323 qkm1 *= biginv;
324 }
325 if ((fabs(qk) < biginv) || (fabs(pk) < biginv)) {
326 pkm2 *= big;
327 pkm1 *= big;
328 qkm2 *= big;
329 qkm1 *= big;
330 }
331 }
332 while (++n < 300);
333 cdone:
334 return (ans);
335 }
336
337 /* Power series for incomplete beta integral.
338 * Use when b*x is small and x not too close to 1. */
339
pseries(a,b,x)340 static double pseries(a, b, x)
341 double a, b, x;
342 {
343 double s, t, u, v, n, t1, z, ai;
344
345 ai = 1.0 / a;
346 u = (1.0 - b) * x;
347 v = u / (a + 1.0);
348 t1 = v;
349 t = u;
350 n = 2.0;
351 s = 0.0;
352 z = MACHEP * ai;
353 while (fabs(v) > z) {
354 u = (n - b) * x / n;
355 t *= u;
356 v = t / (a + n);
357 s += v;
358 n += 1.0;
359 }
360 s += t1;
361 s += ai;
362
363 u = a * log(x);
364 if ((a + b) < MAXGAM && fabs(u) < MAXLOG) {
365 t = 1.0 / beta(a, b);
366 s = s * t * pow(x, a);
367 }
368 else {
369 t = -lbeta(a,b) + u + log(s);
370 if (t < MINLOG)
371 s = 0.0;
372 else
373 s = exp(t);
374 }
375 return (s);
376 }
377