The Beautiful Beta Functions: In Raw Python!

Plots of Complete Beta Functions

For almost a week I have been trying very very hard to bend the computer to code for the incomplete Beta function in Python. Let me show you guys how I finally achieved it using only the most basic of Python functionality, namely, importing only math routine. The reason why I have been busting myself for this is that it (along with gamma functions) is among the foundations of many probability distributions. Once you get through these special functions, the rest of statistical computing will be only about translating from mathematics to code. And this is my sole motivation for this. Hence, I will only discuss the maths of it only superficially. I have used some useful resources to learn and code beta function (please see below: Further reading/Resources).

Complete Beta Function:

Complete Beta function can be defined in terms of gamma function as such:

$\boldsymbol {B(a,b)=B(b,a)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}}$

Proof:
for $\boldsymbol{m}$, given $\boldsymbol{m}$ is an integer, one can write $\boldsymbol {\int_0^{\infty} p^{m-1}e^{-p}dp=(m-1)!}$ for two of such functions, one can write $\boldsymbol {P(m,n)=(m-1)!(n-1)! = \int_0^{\infty}p^{m-1}e{-p}dp\times \int_0^{\infty}q^{n-1}e^{-q}dq}$ $\boldsymbol {= 2\int_0^{\infty}x^{2m-1}e^{-x^2}dx\times 2\int_0^{\infty}y^{2n-1}e^{-y^2}dy.}$ Pass the polar coordinates; $\boldsymbol { x= r \cos \theta,\qquad y=r\sin\theta,\qquad dx\,dy\rightarrow rdr\,d\theta}$ and this becomes $\boldsymbol {P(m,n)= 4\int_0^{\infty}r^{2(m+n-1)}e^{-r^2}rdr\times \int_0^{\pi/2}\cos^{2m-1}\theta\sin^{2n-1}\theta \,d\theta}$ $\boldsymbol {=2(m+n-1)!\int_0^{\pi/2}\cos^{2m-1}\theta\sin^{2n-1}\theta \,d\theta}$ remember $\boldsymbol{P\left (m,n \right ) = \left (m - 1\right )! \left (n - 1 \right )! }$ from above. Hence, $\boldsymbol {2\int_0^{\pi/2}\cos^{2m-1}\theta\sin^{2n-1}\theta \,d\theta=\frac{(m-1)!(n-1)!}{(m+n-1)!}.}$ change variable $\boldsymbol {t=\sin^2\theta}$ so $\boldsymbol {dt=2\sin\theta\cos\theta\,d\theta}$ and $\boldsymbol {\cos^2\theta = 1-t }$ giving, $\boldsymbol{ 2\int_0^{\pi/2}\cos^{2m-1}\theta\sin^{2n-1}\theta \,d\theta=\int_0^{\pi/2}\underbrace{\cos^{2m-2}\theta}_{(1-t)^{m-1}}\times\underbrace{\sin^{2n-2}\theta}_{t^{n-1}}\times\underbrace{2\sin\theta\cos\theta\,d\theta}_{dt}}$ $\boldsymbol{=\int_0^1 t^{n-1}(1-t)^{m-1}dt}$ if we define Beta function in terms of Gamma function like $\boldsymbol {\displaystyle \int_0^{\infty}p^{m-1}e^{-p}dp=\Gamma(m);}$ then from what we have seen above, $\boldsymbol {B(m,n)=\int_0^{1}t^{m-1}(1-t)^{n-1}dt}$ $\boldsymbol {=\frac{\Gamma(m)\Gamma(n)}{\Gamma(m+n)}.}$

To code the complete beta function we can use the functonality of Python math library (the only library that I ever use!). This, I think is not too much of a dependency. Sometime later, I would write a new code for log-gamma and gamma functions. Till then this short but cool code for complete beta would suffice. Lets use both math.gamma() and math.lgamma() functions (implementations of gamma and log-gamma functions respectively in math lib of python)

import math

def beta_1(a,b):

'''uses gamma function or inbuilt math.gamma() to compute values of beta function'''

beta = math.gamma(a)*math.gamma(b)/math.gamma(a+b)
return beta



We can also compute complete beta function using log-gamma function. In most of resources, including ‘Numerical Recipes in C’ I have seen log-gamma being used to compute values of beta function.

import math

def beta_2(a,b):

'''uses log-gamma function or inbuilt math.lgamma() to compute values of beta function'''

beta = math.exp(math.lgamma(a) + math.lgamma(b) - math.lgamma(a+b))
return beta



Incomplete Beta Function:

Incomplete Beta function was what took me so long to code. incomplete beta function can be defined as:

$\boldsymbol { \mathrm{B}(z;a,b)}$

and the regularized beta function is closely related to both incomplete and complete beta function. We will see later in my future posts how this regularized function keeps popping up in probability theory.

$\boldsymbol {I_z(a,b) = \dfrac{\mathrm{B}(z;a,b)}{\mathrm{B}(a,b)}}$


import math

def incompbeta(a, b, x):

''' incompbeta(a,b,x) evaluates incomplete beta function, here a, b > 0 and 0 <= x <= 1. This function requires contfractbeta(a,b,x, ITMAX = 200)
(Code translated from: Numerical Recipes in C.)'''

if (x == 0):
return 0;
elif (x == 1):
return 1;
else:
lbeta = math.lgamma(a+b) - math.lgamma(a) - math.lgamma(b) + a * math.log(x) + b * math.log(1-x)
if (x < (a+1) / (a+b+2)):
return math.exp(lbeta) * contfractbeta(a, b, x) / a;
else:
return 1 - math.exp(lbeta) * contfractbeta(b, a, 1-x) / b;


As you guys might have noted above that the function incompbeta(a,b,x) requires betacontfract(a,b,x, ITMAX), here is the code. Follow the resources below to read in detail about continued fraction method of computing incomplete beta.

import math

def contfractbeta(a,b,x, ITMAX = 200):

""" contfractbeta() evaluates the continued fraction form of the incomplete Beta function; incompbeta().
(Code translated from: Numerical Recipes in C.)"""

EPS = 3.0e-7
bm = az = am = 1.0
qab = a+b
qap = a+1.0
qam = a-1.0
bz = 1.0-qab*x/qap

for i in range(ITMAX+1):
em = float(i+1)
tem = em + em
d = em*(b-em)*x/((qam+tem)*(a+tem))
ap = az + d*am
bp = bz+d*bm
d = -(a+em)*(qab+em)*x/((qap+tem)*(a+tem))
app = ap+d*az
bpp = bp+d*bz
aold = az
am = ap/bpp
bm = bp/bpp
az = app/bpp
bz = 1.0
if (abs(az-aold)<(EPS*abs(az))):
return az

print 'a or b too large or given ITMAX too small for computing incomplete beta function.'


Python Functions Implemented in this Post

beta_1(a,b): uses gamma function or inbuilt math.gamma() to compute values of beta function

beta_2(a,b): uses log-gamma function or inbuilt math.lgamma() to compute values of beta function

incompbeta(a, b, x): incompbeta(a,b,x) evaluates incomplete beta function, for a, b > 0 and 0 <= x <= 1. This function requires contfractbeta(a,b,x, ITMAX = 200)

contfractbeta(a,b,x, ITMAX = 200): evaluates the continued fraction form of the incomplete Beta function; incompbeta().