betai.c File Reference

Routines for computing the incomplete beta function. More...

#include "mdb.h"

Go to the source code of this file.

Macros
#define	BETAI_ACCURACY   1e-10
#define	MAXIMUM_ITERATIONS   10000

Functions
double	betaIncSum (double x, double a, double b)
double	betaInc1 (double x, double a, double b)
double	betaComp (double a, double b)
	Compute the complete beta function component.
double	lnBetaComp (double a, double b)
double	betaInc (double a, double b, double x)
	Compute the incomplete beta function.

Detailed Description

Routines for computing the incomplete beta function.

Copyright

• (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
• (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.

License

This file is distributed under the terms of the Software License Agreement found in the file LICENSE included with this distribution.

Author

M. Borland, R. Soliday

Definition in file betai.c.

Macro Definition Documentation

BETAI_ACCURACY

#define BETAI_ACCURACY   1e-10

Definition at line 27 of file betai.c.

MAXIMUM_ITERATIONS

#define MAXIMUM_ITERATIONS   10000

Definition at line 89 of file betai.c.

Function Documentation

betaComp()

double betaComp	(	double	a,
		double	b )

Compute the complete beta function component.

Uses the gamma function to compute the value of Beta(a,b) = Γ(a)*Γ(b)/Γ(a+b).

Parameters

a	First parameter.
b	Second parameter.

Returns

The value of Beta(a,b).

Definition at line 38 of file betai.c.

                                    {
#if defined(vxWorks)
  fprintf(stderr, "lgamma function not implemented on vxWorks\n");
  exit(1);
#else
  return exp(lgamma(a) + lgamma(b) - lgamma(a + b));
#endif
}

betaInc()

double betaInc	(	double	a,
		double	b,
		double	x )

Compute the incomplete beta function.

Calculates the incomplete beta function I_x(a,b) for given parameters a, b and x. If necessary, parameters may be swapped internally to improve convergence.

Parameters

a	First parameter (must be > 0).
b	Second parameter (must be > 0).
x	The integration limit (0 ≤ x ≤ 1).

Returns

The value of I_x(a,b), or -1.0 if x is out of range.

Definition at line 67 of file betai.c.

                                             {
  double xLimit, sum, lnBab = 0;
  short swapped = 0;
 
  if (x < 0 || x > 1)
    return -1.0;
  if (a + b != 2) {
    xLimit = (a + 1) / (a + b - 2);
    if (x > xLimit) {
      swapped = 1;
      x = 1 - x;
      SWAP_DOUBLE(a, b);
    }
  }
  lnBab = lnBetaComp(a, b);
  /*    sum = pow(x, a)*pow(1-x, b)*betaIncSum(a, b, x)/(Bab*a); */
  sum = exp(a * log(x) + b * log(1 - x) - lnBab) * betaIncSum(a, b, x) / a;
  if (!swapped)
    return sum;
  return 1 - sum;
}

betaInc1()

double betaInc1	(	double	x,
		double	a,
		double	b )

Definition at line 139 of file betai.c.

                                              {
  double xp, sum, term;
  long i;
  xp = x;
  i = sum = 0;
  do {
    sum += (term = betaComp(a + 1, i + 1.0) / betaComp(a + b, i + 1.0) * xp);
    xp = x * xp;
    i++;
  } while (term > BETAI_ACCURACY);
  return (sum + 1) * pow(x, a) * pow(1 - x, b) / (a * betaComp(a, b));
}

betaIncSum()

double betaIncSum	(	double	x,
		double	a,
		double	b )

Definition at line 91 of file betai.c.

                                                {
  double f1, f2;
  double A0, B0, A_1, B_1, A1, B1, d;
  double m, aM1, aP1, aPb, mT2;
 
  /* evaluate the first step (d1) outside the loop to get started */
  A_1 = B_1 = A0 = 1;
  B0 = 1 - (a + b) / (a + 1) * x;
  aPb = a + b;
  aP1 = a + 1;
  aM1 = a - 1;
  m = 1;
  do {
    mT2 = m * 2;
    d = m * (b - m) * x / ((aM1 + mT2) * (a + mT2));
    A1 = A0 + d * A_1;
    B1 = B0 + d * B_1;
    f1 = A1 / B1;
    A_1 = A0;
    B_1 = B0;
    A0 = A1;
    B0 = B1;
    d = -(a + m) * (aPb + m) * x / ((a + mT2) * (aP1 + mT2));
    A1 = A0 + d * A_1;
    B1 = B0 + d * B_1;
    f2 = A1 / B1;
    A_1 = A0;
    B_1 = B0;
    A0 = A1;
    B0 = B1;
    if (B1) {
      /* rescale the recurrence to avoid over/underflow */
      A_1 /= B1;
      B_1 /= B1;
      A0 /= B1;
      B0 = 1;
    }
    m++;
  } while (fabs(f1 - f2) > BETAI_ACCURACY && m < MAXIMUM_ITERATIONS);
  /*
    if (m>=MAXIMUM_ITERATIONS)
        fprintf(stderr, "warning: too many iterations for incomplete beta function sum (%e, %e, %e) (betaInc)\n",
                a, b, x);
*/
  return f2;
}

lnBetaComp()

double lnBetaComp	(	double	a,
		double	b )

Definition at line 47 of file betai.c.

                                      {
#if defined(vxWorks)
  fprintf(stderr, "lgamma function not implemented on vxWorks\n");
  exit(1);
#else
  return lgamma(a) + lgamma(b) - lgamma(a + b);
#endif
}