2
$\begingroup$

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.

$\endgroup$

marked as duplicate by whuber Mar 25 at 14:11

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • 1
    $\begingroup$ I recently posted an efficient, accurate, well-tested implementation of this incomplete Beta function at stats.stackexchange.com/a/394983/919. It is arguably the "standard" way you are asking about. $\endgroup$ – whuber Mar 25 at 14:11
1
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ Numerical integration will work fine for large parameters (albeit inefficiently), but if either one of them is less than $1,$ it may fail badly because of the singularity of the density near $0$ or $1.$ $\endgroup$ – whuber Mar 25 at 14:13
  • $\begingroup$ @whuber, I’m not surprised that my pragmatic engineering approach is different than a true statistician :) $\endgroup$ – gunes Mar 25 at 15:00
  • 1
    $\begingroup$ I learned my numerical computing skills by making mistakes. You remember those! I spent a few weeks one summer debugging intricate code only to discover that I had been just a little sloppy in my implementation of a Bessel function, and propagation of a small numerical error in a few special cases was ruining everything. It was a good lesson. $\endgroup$ – whuber Mar 25 at 17:05

Not the answer you're looking for? Browse other questions tagged or ask your own question.