CodeCogs

The regularized incomplete Beta integral

You're viewing an older version of this page (#493). View the current version.

Interface

#include <codecogs/maths/special/gamma/beta_reg.h>

using namespace Maths::Special::Gamma;

double fraction_exp(double x, double a, double b, bool inverse=false)
double betaLower_reg_pow(double x, double a, double b)
double beta_reg(double x, double a, double b, bool upper = false)

Overview

These function compute the regularized incomplete beta integral of the arguments, evaluated from zero to x. This is occasionally also called the beta function ratio, given that it returns the probability that a random variable drawn from a beta distribution with parameters a and b will be less that or equal to x.

This ratio (or regularization) function is calculated by dividing the incomplete beta integral by the complete beta integral, i.e.

If the incomplete beta integral is defined by:

$B_z(a,b) = \int_0^x t^{a-1} (1-t)^{b-1} dt$

(1)

and the complete beta integral is defined by:

$\frac{\Gamma(a) \Gamma(b)}{\Gamma(a+b)}$

(2)

where $\Gamma$ is the Gamma function.

Therefore the regularized beta integral is:

$\frac{\Gamma(a+b)}{\Gamma(a) \Gamma(b)} * B_x(a,b)$

(3)

where $B_x(a,b)$ .

All the functions are only valid in the range 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, i.e.

1 - beta_reg( x, a, b )  =  beta_reg( 1-x, a, b ).

Example:

The following example evaluates the upper regularized incomplete beta integral for 10 points equally spaced in the interval from 0 to 1.

#include <stdio.h>
#include <codecogs/maths/special/gamma/beta_reg.h>
void main()
{
  for (double x=0; x<1; x+=0.1)
  {
    double y = Maths::Special::Gamma::beta_reg(x, 2, 1, true);
    printf("beta_reg(2, %.1lf, 1, true) = %lf\n", x,y);
  }
}

Output:

beta_reg(2, 0.0, 1, true) = 1.000000
beta_reg(2, 0.1, 1, true) = 0.990000
beta_reg(2, 0.2, 1, true) = 0.960000
beta_reg(2, 0.3, 1, true) = 0.910000
beta_reg(2, 0.4, 1, true) = 0.840000
beta_reg(2, 0.5, 1, true) = 0.750000
beta_reg(2, 0.6, 1, true) = 0.640000
beta_reg(2, 0.7, 1, true) = 0.510000
beta_reg(2, 0.8, 1, true) = 0.360000
beta_reg(2, 0.9, 1, true) = 0.190000
beta_reg(2, 1.0, 1, true) = 0.000000

FUNCTION

fraction_exp

Computes using continued fraction expansion that main term need to compute the incomplete Beta integral at the extremes.

{Methodolgy} The overall incomplete Beta Integral can be computed using fraction expansion, such that

$w \cdot x^a (1-x)^b \frac{ \Gamma (a+b)}{a \Gamma a \Gamma b}$

(4)

where w is the fraction term evaluated here as

$w= \frav{\sum p_n}{\sum q_n}$

(5)

and

$p_{n+1} = p_{n-1}(1+qk_n} + p_{n-2} pk_n$

(6)

$q_{n+1} = q_{n-1}(1+qk_n} + q_{n-2} pk_n$

(7)

where the nth term of the numerator is

$pk_n=\frac{x (a+n)(a+n+n)}{(a+1+2n)(a+2n)}$

(8)

and the denominator is

$qk_n=\frac{x (n+1)(b-1-n)}{(a+1+2n)(a+2+2n) }$

(9)

with initial values of

$p_0 = 1, q_0=1, p_{-1}=0, q_{-1}=1$

(10)

Parameters

upper limit of integration. (0<x<1)

first shape argument, must be positive.

second shape argument, must be positive.

inverse

computes the fraction expansion $x/(1-x)$

Interactive Calculator

inverse

Result –

FUNCTION

betaLower_reg_pow

Power series for the regularized incomplete beta integral, $B_x(a,b)$ , Use when $b x$ is little and x is not too close to 1.

Parameters

upper limit of integration (0<x<1)

first shape argument, must be positive

second shape argument, must be positive

Author

Interactive Calculator

Result –

FUNCTION

beta_reg

Computes the regularize beta integral using a variety of method depending on the input values:

When b*x<=1 and x<0.95 then the function is calculated using the power series expansion defined in betaLower_reg_pow.
If x is greater than the mean of the beta function, i.e.
$x >\frac{a}{a+b}$
(11)
Then the values of a and b are reversed and the calculation proceeds as before only with x=1-x.
A second attempt is then made to use the power series expansion within the range b*x<=1 and x<0.95.
Failing this a fractional expansion is used within the tails of this function.

Parameters

the value at which to evaluate the function, must be in range 0..1

the 1st number of degrees of freedom, must be strictly positive

the 2nd number of degrees of freedom, must be strictly positive

upper

Default value = false

Author

Interactive Calculator

upper

Result –