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:
for , given is an integer, one can write for two of such functions, one can write Pass the polar coordinates; and this becomes remember from above. Hence, change variable so and giving, if we define Beta function in terms of Gamma function like then from what we have seen above,
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:
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.
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().
- Lanczos, C. 1964 SIAM Journal on Numerical Analysis, ser. B, vol. 1, pp. 86-96
- Numerical Recipes: The Art of Scientific Computing, Third Edition (2007)ISBN-10: 0521880688)