Efficiently Computing The Beta CDF

Question

I am using numba to JIT compile some looped python functions as part of a larger application. Ideally, everything will run in numba's "no python" mode, such that the loop can be parallelised.

One of these functions involves use of the Beta distribution CDF (the regularised incomplete beta function):

$$I_z(z| \alpha, \beta) = \frac{B(z | \alpha, \beta)}{B(\alpha, \beta)}$$

where

$$B(\alpha, \beta) =\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)}$$

and

$$B(z|\alpha, \beta) = \int_0^{z} u ^{\alpha -1} (1 -u)^{\beta-1} du$$

To run in nopython mode, I would like to implement the above using only functions from python's standard math library and numpy.

Implementing $B(\alpha, \beta)$ is no problem as I can use math.gamma from the python standard library.

The problem is how to evaluate the integral for $B(z|\alpha, \beta)$, which -- based on my admittedly limited knowledge of the problem -- must be performed numerically.

Is there a "standard" way that this is implemented inside most statistical libraries, eg. scipy? This question appears to ask a similar question but remains unanswered.

Please note that $\alpha$ and $\beta$ could be varied quite substantially, so a lookup table with interpolation is not an option.

Solution for those searching for similar: I ended up following this linked answer from the comment below by whuber. Doing some further reading lead me to this implementation on Github, along with the author's explanation of the algorithm.

Answer 1

There are several ways to do numeric integration, you just need to setup the resolution. One way is to use naive implementation:

v = np.linspace(0,z,1000)
dv = v[1]-v[0]
np.sum(v**(a-1) * (1-v)**(b-1) * dv)

which gives 0.055928172212507204 when z=0.5,a=3,b=0.4 for example.

Or, use numpy.trapz, a slighly more accurate way of integration due to using trapezoids, instead of rectangles (which gives 0.05583334509998778):

np.trapz(v**(a-1) * (1-v)**(b-1), v)

and you can compare them with scipy.integrate:

import scipy.integrate as integrate
integrate.quad(lambda x: x**(a-1) * (1-x)**(b-1), 0, z)[0]

which gives 0.055833303967108956.