/* incbetf.c * * Incomplete beta integral * * * SYNOPSIS: * * float a, b, x, y, incbetf(); * * y = incbetf( a, b, x ); * * * DESCRIPTION: * * Returns incomplete beta integral of the arguments, evaluated * from zero to x. The function is defined as * * x * - - * | (a+b) | | a-1 b-1 * ----------- | t (1-t) dt. * - - | | * | (a) | (b) - * 0 * * The domain of definition is 0 <= x <= 1. In this * implementation a and b are restricted to positive values. * The integral from x to 1 may be obtained by the symmetry * relation * * 1 - incbet( a, b, x ) = incbet( b, a, 1-x ). * * The integral is evaluated by a continued fraction expansion. * If a < 1, the function calls itself recursively after a * transformation to increase a to a+1. * * ACCURACY: * * Tested at random points (a,b,x) with a and b in the indicated * interval and x between 0 and 1. * * arithmetic domain # trials peak rms * Relative error: * IEEE 0,30 10000 3.7e-5 5.1e-6 * IEEE 0,100 10000 1.7e-4 2.5e-5 * The useful domain for relative error is limited by underflow * of the single precision exponential function. * Absolute error: * IEEE 0,30 100000 2.2e-5 9.6e-7 * IEEE 0,100 10000 6.5e-5 3.7e-6 * * Larger errors may occur for extreme ratios of a and b. * * ERROR MESSAGES: * message condition value returned * incbetf domain x<0, x>1 0.0 */ /* Cephes Math Library, Release 2.2: July, 1992 Copyright 1984, 1987, 1992 by Stephen L. Moshier Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */ #include #ifdef ANSIC float lgamf(float), expf(float), logf(float); static float incbdf(float, float, float); static float incbcff(float, float, float); float incbpsf(float, float, float); #else float lgamf(), expf(), logf(); float incbpsf(); static float incbcff(), incbdf(); #endif #define fabsf(x) ( (x) < 0 ? -(x) : (x) ) /* BIG = 1/MACHEPF */ #define BIG 16777216. extern float MACHEPF, MAXLOGF; #define MINLOGF (-MAXLOGF) float incbetf( float aaa, float bbb, float xxx ) { float aa, bb, xx, ans, a, b, t, x, onemx; int flag; aa = aaa; bb = bbb; xx = xxx; if( (xx <= 0.0) || ( xx >= 1.0) ) { if( xx == 0.0 ) return(0.0); if( xx == 1.0 ) return( 1.0 ); mtherr( "incbetf", DOMAIN ); return( 0.0 ); } onemx = 1.0 - xx; /* transformation for small aa */ if( aa <= 1.0 ) { ans = incbetf( aa+1.0, bb, xx ); t = aa*logf(xx) + bb*logf( 1.0-xx ) + lgamf(aa+bb) - lgamf(aa+1.0) - lgamf(bb); if( t > MINLOGF ) ans += expf(t); return( ans ); } /* see if x is greater than the mean */ if( xx > (aa/(aa+bb)) ) { flag = 1; a = bb; b = aa; t = xx; x = onemx; } else { flag = 0; a = aa; b = bb; t = onemx; x = xx; } /* transformation for small aa */ /* if( a <= 1.0 ) { ans = a*logf(x) + b*logf( onemx ) + lgamf(a+b) - lgamf(a+1.0) - lgamf(b); t = incbetf( a+1.0, b, x ); if( ans > MINLOGF ) t += expf(ans); goto bdone; } */ /* Choose expansion for optimal convergence */ if( b > 10.0 ) { if( fabsf(b*x/a) < 0.3 ) { t = incbpsf( a, b, x ); goto bdone; } } ans = x * (a+b-2.0)/(a-1.0); if( ans < 1.0 ) { ans = incbcff( a, b, x ); t = b * logf( t ); } else { ans = incbdf( a, b, x ); t = (b-1.0) * logf(t); } t += a*logf(x) + lgamf(a+b) - lgamf(a) - lgamf(b); t += logf( ans/a ); if( t < MINLOGF ) { t = 0.0; if( flag == 0 ) { mtherr( "incbetf", UNDERFLOW ); } } else { t = expf(t); } bdone: if( flag ) t = 1.0 - t; return( t ); } /* Continued fraction expansion #1 * for incomplete beta integral */ static float incbcff( float aa, float bb, float xx ) { float a, b, x, xk, pk, pkm1, pkm2, qk, qkm1, qkm2; float k1, k2, k3, k4, k5, k6, k7, k8; float r, t, ans; static float big = BIG; int n; a = aa; b = bb; x = xx; k1 = a; k2 = a + b; k3 = a; k4 = a + 1.0; k5 = 1.0; k6 = b - 1.0; k7 = k4; k8 = a + 2.0; pkm2 = 0.0; qkm2 = 1.0; pkm1 = 1.0; qkm1 = 1.0; ans = 1.0; r = 0.0; n = 0; do { xk = -( x * k1 * k2 )/( k3 * k4 ); pk = pkm1 + pkm2 * xk; qk = qkm1 + qkm2 * xk; pkm2 = pkm1; pkm1 = pk; qkm2 = qkm1; qkm1 = qk; xk = ( x * k5 * k6 )/( k7 * k8 ); pk = pkm1 + pkm2 * xk; qk = qkm1 + qkm2 * xk; pkm2 = pkm1; pkm1 = pk; qkm2 = qkm1; qkm1 = qk; if( qk != 0 ) r = pk/qk; if( r != 0 ) { t = fabsf( (ans - r)/r ); ans = r; } else t = 1.0; if( t < MACHEPF ) goto cdone; k1 += 1.0; k2 += 1.0; k3 += 2.0; k4 += 2.0; k5 += 1.0; k6 -= 1.0; k7 += 2.0; k8 += 2.0; if( (fabsf(qk) + fabsf(pk)) > big ) { pkm2 *= MACHEPF; pkm1 *= MACHEPF; qkm2 *= MACHEPF; qkm1 *= MACHEPF; } if( (fabsf(qk) < MACHEPF) || (fabsf(pk) < MACHEPF) ) { pkm2 *= big; pkm1 *= big; qkm2 *= big; qkm1 *= big; } } while( ++n < 100 ); cdone: return(ans); } /* Continued fraction expansion #2 * for incomplete beta integral */ static float incbdf( float aa, float bb, float xx ) { float a, b, x, xk, pk, pkm1, pkm2, qk, qkm1, qkm2; float k1, k2, k3, k4, k5, k6, k7, k8; float r, t, ans, z; static float big = BIG; int n; a = aa; b = bb; x = xx; k1 = a; k2 = b - 1.0; k3 = a; k4 = a + 1.0; k5 = 1.0; k6 = a + b; k7 = a + 1.0;; k8 = a + 2.0; pkm2 = 0.0; qkm2 = 1.0; pkm1 = 1.0; qkm1 = 1.0; z = x / (1.0-x); ans = 1.0; r = 0.0; n = 0; do { xk = -( z * k1 * k2 )/( k3 * k4 ); pk = pkm1 + pkm2 * xk; qk = qkm1 + qkm2 * xk; pkm2 = pkm1; pkm1 = pk; qkm2 = qkm1; qkm1 = qk; xk = ( z * k5 * k6 )/( k7 * k8 ); pk = pkm1 + pkm2 * xk; qk = qkm1 + qkm2 * xk; pkm2 = pkm1; pkm1 = pk; qkm2 = qkm1; qkm1 = qk; if( qk != 0 ) r = pk/qk; if( r != 0 ) { t = fabsf( (ans - r)/r ); ans = r; } else t = 1.0; if( t < MACHEPF ) goto cdone; k1 += 1.0; k2 -= 1.0; k3 += 2.0; k4 += 2.0; k5 += 1.0; k6 += 1.0; k7 += 2.0; k8 += 2.0; if( (fabsf(qk) + fabsf(pk)) > big ) { pkm2 *= MACHEPF; pkm1 *= MACHEPF; qkm2 *= MACHEPF; qkm1 *= MACHEPF; } if( (fabsf(qk) < MACHEPF) || (fabsf(pk) < MACHEPF) ) { pkm2 *= big; pkm1 *= big; qkm2 *= big; qkm1 *= big; } } while( ++n < 100 ); cdone: return(ans); } /* power series */ float incbpsf( float aa, float bb, float xx ) { float a, b, x, t, u, y, s; a = aa; b = bb; x = xx; y = a * logf(x) + (b-1.0)*logf(1.0-x) - logf(a); y -= lgamf(a) + lgamf(b); y += lgamf(a+b); t = x / (1.0 - x); s = 0.0; u = 1.0; do { b -= 1.0; if( b == 0.0 ) break; a += 1.0; u *= t*b/a; s += u; } while( fabsf(u) > MACHEPF ); if( y < MINLOGF ) { mtherr( "incbetf", UNDERFLOW ); s = 0.0; } else s = expf(y) * (1.0 + s); /*printf( "incbpsf: %.4e\n", s );*/ return(s); }