scipy.stats.chi2.ppf() without scipy

Allantools has a dependence on scipy because it uses scipy.stats.chi2.ppf() (inverse of the chi-squared cumulative distribution function) in the code for confidence intervals.

When testing scipy takes quite a long time to install and the whole package seems a bit overkill for just this one function.

So I tried implementing it using just numpy. It sort of works, but not for corner cases where p is very close to 1.0 or k is a big number. I get for example:

Maybe there's something trick when numerically evaluating lower_gamma() using the series expansion that I am missing...

4 replies on “scipy.stats.chi2.ppf() without scipy”

Try computing dg using log-math. not sure if this will render in a comment or not:

however, this has some failures of the 'isclose' test e.g., at
190 0.9999 False 271.179788314 271.221549805
(relative error .00015)

full test program adapted from yours at

boost discusses their implementation: using several different methods depending on the values of (what they call) x and a, may be worth a gander.

Thanks! I will take a look at the links.

HTML tags 'pre lang="python"' around the code makes it render using the WP-Syntax plugin.

You could maybe get away with exploiting the fact that the gamma function in the denominator of the pdf of a chi-squared random variable is just a normalizing constant. So, e.g., these three functions get me pretty close to the output of scipy.stats.chi2.ppf():

The "binsize" of x matters, both for speed and accuracy. And if you don't know what k is ahead of time, the fact that you have to define x is kind of a pain, though you can probably also exploit the fact that the expected value of a chi-squared random variable is k and the variance is 2k to define x over an appropriate range once k is known.

Here's hoping I used the "pre" tag correctly...

I have a feeling that the cdf(x,k) function is very flat at high x as we get closer to p=1.0. Conversely ppf(p,k) is very steep as we get closer to p=1.0

One improvement could be to change the bisection() function that looks for the right value of x to instead stop searching with a tolerance check on p (instead of x). Maybe?

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.