/* incbet.c
*
* Incomplete beta integral
*
*
* SYNOPSIS:
*
* double a, b, x, y, incbet();
*
* y = incbet( 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
* or, when b*x is small, by a power series.
*
* ACCURACY:
*
* Tested at uniformly distributed random points (a,b,x) with a and b
* in "domain" and x between 0 and 1.
* Relative error
* arithmetic domain # trials peak rms
* IEEE 0,5 10000 6.9e-15 4.5e-16
* IEEE 0,85 250000 2.2e-13 1.7e-14
* IEEE 0,1000 30000 5.3e-12 6.3e-13
* IEEE 0,10000 250000 9.3e-11 7.1e-12
* IEEE 0,100000 10000 8.7e-10 4.8e-11
* Outputs smaller than the IEEE gradual underflow threshold
* were excluded from these statistics.
*
* ERROR MESSAGES:
* message condition value returned
* incbet domain x<0, x>1 0.0
* incbet underflow 0.0
*/
/*
Cephes Math Library, Release 2.3: March, 1995
Copyright 1984, 1995 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef DEC
#define MAXGAM 34.84425627277176174
#else
#define MAXGAM 171.624376956302725
#endif
extern double MACHEP, MINLOG, MAXLOG;
static double incbcf (double, double, double);
static double incbd (double, double, double);
static double pseries (double, double, double);
static double big = 4.503599627370496e15;
static double biginv = 2.22044604925031308085e-16;
double
incbet (double aa, double bb, double xx)
{
double a, b, t, x, xc, w, y;
int flag;
if (aa <= 0.0 || bb <= 0.0)
goto domerr;
if ((xx <= 0.0) || (xx >= 1.0)) {
if (xx == 0.0)
return (0.0);
if (xx == 1.0)
return (1.0);
domerr:
mtherr ("incbet", MATHERR_DOMAIN);
return (0.0);
}
flag = 0;
if ((bb * xx) <= 1.0 && xx <= 0.95) {
t = pseries (aa, bb, xx);
goto done;
}
w = 1.0 - xx;
/* Reverse a and b if x is greater than the mean. */
if (xx > (aa / (aa + bb))) {
flag = 1;
a = bb;
b = aa;
xc = xx;
x = w;
} else {
a = aa;
b = bb;
xc = w;
x = xx;
}
if (flag == 1 && (b * x) <= 1.0 && x <= 0.95) {
t = pseries (a, b, x);
goto done;
}
/* Choose expansion for better convergence. */
y = x * (a + b - 2.0) - (a - 1.0);
if (y < 0.0)
w = incbcf (a, b, x);
else
w = incbd (a, b, x) / xc;
/* Multiply w by the factor
a b _ _ _
x (1-x) | (a+b) / ( a | (a) | (b) ) . */
y = a * log (x);
t = b * log (xc);
if ((a + b) < MAXGAM && fabs (y) < MAXLOG && fabs (t) < MAXLOG) {
t = pow (xc, b);
t *= pow (x, a);
t /= a;
t *= w;
t *= gammafn (a + b) / (gammafn (a) * gammafn (b));
goto done;
}
/* Resort to logarithms. */
y += t + lgam (a + b) - lgam (a) - lgam (b);
y += log (w / a);
if (y < MINLOG)
t = 0.0;
else
t = exp (y);
done:
if (flag == 1) {
if (t <= MACHEP)
t = 1.0 - MACHEP;
else
t = 1.0 - t;
}
return (t);
}
/* Continued fraction expansion #1
* for incomplete beta integral
*/
static double
incbcf (double a, double b, double x)
{
double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
double k1, k2, k3, k4, k5, kk6, kk7, k8;
double r, t, ans, thresh;
int n;
k1 = a;
k2 = a + b;
k3 = a;
k4 = a + 1.0;
k5 = 1.0;
kk6 = b - 1.0;
kk7 = k4;
k8 = a + 2.0;
pkm2 = 0.0;
qkm2 = 1.0;
pkm1 = 1.0;
qkm1 = 1.0;
ans = 1.0;
r = 1.0;
n = 0;
thresh = 3.0 * MACHEP;
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 * kk6) / (kk7 * 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 = fabs ((ans - r) / r);
ans = r;
} else
t = 1.0;
if (t < thresh)
goto cdone;
k1 += 1.0;
k2 += 1.0;
k3 += 2.0;
k4 += 2.0;
k5 += 1.0;
kk6 -= 1.0;
kk7 += 2.0;
k8 += 2.0;
if ((fabs (qk) + fabs (pk)) > big) {
pkm2 *= biginv;
pkm1 *= biginv;
qkm2 *= biginv;
qkm1 *= biginv;
}
if ((fabs (qk) < biginv) || (fabs (pk) < biginv)) {
pkm2 *= big;
pkm1 *= big;
qkm2 *= big;
qkm1 *= big;
}
}
while (++n < 300);
cdone:
return (ans);
}
/* Continued fraction expansion #2
* for incomplete beta integral
*/
static double
incbd (double a, double b, double x)
{
double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
double k1, k2, k3, k4, k5, kk6, kk7, k8;
double r, t, ans, z, thresh;
int n;
k1 = a;
k2 = b - 1.0;
k3 = a;
k4 = a + 1.0;
k5 = 1.0;
kk6 = a + b;
kk7 = 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 = 1.0;
n = 0;
thresh = 3.0 * MACHEP;
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 * kk6) / (kk7 * 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 = fabs ((ans - r) / r);
ans = r;
} else
t = 1.0;
if (t < thresh)
goto cdone;
k1 += 1.0;
k2 -= 1.0;
k3 += 2.0;
k4 += 2.0;
k5 += 1.0;
kk6 += 1.0;
kk7 += 2.0;
k8 += 2.0;
if ((fabs (qk) + fabs (pk)) > big) {
pkm2 *= biginv;
pkm1 *= biginv;
qkm2 *= biginv;
qkm1 *= biginv;
}
if ((fabs (qk) < biginv) || (fabs (pk) < biginv)) {
pkm2 *= big;
pkm1 *= big;
qkm2 *= big;
qkm1 *= big;
}
}
while (++n < 300);
cdone:
return (ans);
}
/* Power series for incomplete beta integral.
Use when b*x is small and x not too close to 1. */
static double
pseries (double a, double b, double x)
{
double s, t, u, v, n, t1, z, ai;
ai = 1.0 / a;
u = (1.0 - b) * x;
v = u / (a + 1.0);
t1 = v;
t = u;
n = 2.0;
s = 0.0;
z = MACHEP * ai;
while (fabs (v) > z) {
u = (n - b) * x / n;
t *= u;
v = t / (a + n);
s += v;
n += 1.0;
}
s += t1;
s += ai;
u = a * log (x);
if ((a + b) < MAXGAM && fabs (u) < MAXLOG) {
t = gammafn (a + b) / (gammafn (a) * gammafn (b));
s = s * t * pow (x, a);
} else {
t = lgam (a + b) - lgam (a) - lgam (b) + u + log (s);
if (t < MINLOG)
s = 0.0;
else
s = exp (t);
}
return (s);
}
syntax highlighted by Code2HTML, v. 0.9.1