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