3
* Inverse of imcomplete beta integral
9
* double a, b, x, y, incbi();
11
* x = incbi( a, b, y );
17
* Given y, the function finds x such that
19
* incbet( a, b, x ) = y .
21
* The routine performs interval halving or Newton iterations to find the
22
* root of incbet(a,b,x) - y = 0.
29
* arithmetic domain domain # trials peak rms
30
* IEEE 0,1 .5,10000 50000 5.8e-12 1.3e-13
31
* IEEE 0,1 .25,100 100000 1.8e-13 3.9e-15
32
* IEEE 0,1 0,5 50000 1.1e-12 5.5e-15
33
* VAX 0,1 .5,100 25000 3.5e-14 1.1e-15
34
* With a and b constrained to half-integer or integer values:
35
* IEEE 0,1 .5,10000 50000 5.8e-12 1.1e-13
36
* IEEE 0,1 .5,100 100000 1.7e-14 7.9e-16
37
* With a = .5, b constrained to half-integer or integer values:
38
* IEEE 0,1 .5,10000 10000 8.3e-11 1.0e-11
43
Cephes Math Library Release 2.4: March,1996
44
Copyright 1984, 1996 by Stephen L. Moshier
49
extern double MACHEP, MAXNUM, MAXLOG, MINLOG;
51
double ndtri(), exp(), fabs(), log(), sqrt(), lgam(), incbet();
54
double incbi( aa, bb, yy0 )
57
double a, b, y0, d, y, x, x0, x1, lgm, yp, di, dithresh, yl, yh, xt;
58
int i, rflg, dir, nflg;
72
if( aa <= 1.0 || bb <= 1.0 )
80
y = incbet( a, b, x );
87
/* approximation to inverse function */
107
lgm = (yp * yp - 3.0)/6.0;
108
x = 2.0/( 1.0/(2.0*a-1.0) + 1.0/(2.0*b-1.0) );
109
d = yp * sqrt( x + lgm ) / x
110
- ( 1.0/(2.0*b-1.0) - 1.0/(2.0*a-1.0) )
111
* (lgm + 5.0/6.0 - 2.0/(3.0*x));
118
x = a/( a + b * exp(d) );
119
y = incbet( a, b, x );
124
/* Resort to interval halving if not close enough. */
129
for( i=0; i<100; i++ )
133
x = x0 + di * (x1 - x0);
139
x = x0 + di * (x1 - x0);
143
y = incbet( a, b, x );
144
yp = (x1 - x0)/(x1 + x0);
145
if( fabs(yp) < dithresh )
148
if( fabs(yp) < dithresh )
161
di = 1.0 - (1.0 - di) * (1.0 - di);
165
di = (y0 - y)/(yh - yl);
184
y = incbet( a, b, x );
195
if( rflg == 1 && x1 < MACHEP )
211
di = (y - y0)/(yh - yl);
215
mtherr( "incbi", PLOSS );
224
mtherr( "incbi", UNDERFLOW );
234
lgm = lgam(a+b) - lgam(a) - lgam(b);
238
/* Compute the function at this point. */
261
if( x == 1.0 || x == 0.0 )
263
/* Compute the derivative of the function at this point. */
264
d = (a - 1.0) * log(x) + (b - 1.0) * log(1.0-x) + lgm;
270
/* Compute the step to the next approximation of x. */
275
y = (x - x0) / (x1 - x0);
276
xt = x0 + 0.5 * y * (x - x0);
282
y = (x1 - x) / (x1 - x0);
283
xt = x1 - 0.5 * y * (x1 - x);
288
if( fabs(d/x) < 128.0 * MACHEP )
291
/* Did not converge. */
292
dithresh = 256.0 * MACHEP;