/* 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 <math.h>
#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);
}