scipy.stats.chi2.ppf() without scipy

    top = pow(x, s)*math.exp(-x)*pow(x,k)
OverflowError: (34, 'Numerical result out of range')

import scipy.stats
import math
import numpy
 
def bisection(function, k, p, a, b, tol):
	# http://code.activestate.com/recipes/578417-bisection-method-in-python/
	assert (function(a,k)-p)*(function(b,k)-p) < 0 # a, b must bracket root
	c = (a+b)/2.0
	while (b-a)/2.0 > tol:
		if (function(c, k)-p) == 0:
			return c
		elif (function(a,k)-p)*(function(c,k)-p) < 0:
			b = c
		else :
			a = c
		c = (a+b)/2.0
	return c
 
def lower_gamma(s,x):
	# lower incomplete gamma function
	# https://en.wikipedia.org/wiki/Incomplete_gamma_function#Evaluation_formulae
	g  = 0
	last_g = 1.0
	done = False
	tol = 1.0e-6
	k=0
	while not done:
		top = pow(x, s)*math.exp(-x)*pow(x,k)
		bot = numpy.prod( [float(s+j) for j in range(k+1) ] )
		dg = float(top)/float(bot)
		if dg == float("Inf"):
			break
		g += dg
		k += 1
		if k>100: # get at least 100 terms in the sum
			if g==0:
				break
			delta = abs(dg/g)
			if delta == float("Inf"):
				break
			if delta < tol:
				done = True
		last_g = g
	return g
 
def chi2_cdf(x, k):
	# chi-squared cumulative density function
	# cdf(x; k) = lower_gamma(k/2, x/2) / gamma(k/2)
	return lower_gamma(k/2.0, x/2.0) / math.gamma(k/2.0)
 
def chi2_ppf(p, k):
	# chi-squared Percent point function (inverse of cdf percentiles).
	# look for x such that
	# p = chi2_cdf( x=chi2_ppf(p, k), k)
	tol = 1e-8
	lolim = 0
	hilim = k
	while (chi2_cdf(lolim,k)-p)*(chi2_cdf(hilim,k)-p) > 0: 
		hilim *= 1.5
	return bisection( chi2_cdf, k, p, lolim, hilim, tol)
 
print "scipy cdf: ",scipy.stats.chi2.cdf(55, 33)	
print "own   cdf: ",chi2_cdf(55, 33)
 
 
print "scipy ppf ", scipy.stats.chi2.ppf(0.4, 33)
print "  own ppf ", chi2_ppf(0.4, 33)
 
# test that we really found the inverse
print scipy.stats.chi2.cdf(scipy.stats.chi2.ppf(0.4, 33), 33)
print chi2_cdf( chi2_ppf(0.4, 33), 33 )
 
# try to check the scipy function against our own function
# for some random input of (p, k)
for n in range(100):
	k = numpy.random.randint(20, 200)
	p = numpy.random.random()
	print k, p,
	a=scipy.stats.chi2.ppf(p, k)
	b=chi2_ppf(p, k)
	ok = numpy.isclose(a, b)
	if ok:
		print  ok
	else:
		print  ok, a, b
	assert ok

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

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

   while not done: if x != 0: log_top = math.log(x) * (s + k) - x else: log_top = 0 log_bot = sum(math.log(s+j) for j in range(k+1)) dg = math.exp(log_top - log_bot)

   while not done: if x != 0: log_top = math.log(x) * (s + k) - x else: log_top = 0 log_bot = sum(math.log(s+j) for j in range(k+1)) dg = math.exp(log_top - log_bot)

   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 https://emergent.unpythonic.net/files/sandbox/chi2_ppf_jepler.py

   boost discusses their implementation: http://www.boost.org/doc/libs/1_35_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_gamma/igamma.html#math_toolkit.special.sf_gamma.igamma.implementation using several different methods depending on the values of (what they call) x and a, may be worth a gander.

   Reply

   1. 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.

      Reply

2. 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():

   def chi2_pdf(x,k): bin = x[1]-x[0] fx = (x**((k/2)-1)*np.exp(-x/2))/(2**(k/2)) fx = fx/(bin*np.sum(fx)) return fx, bin   def chi2_cdf(x,k): pdf,bin = chi2_pdf(x,k) cdf = np.zeros(len(x)) for i in range(len(x)): cdf[i] = bin*np.sum(pdf[:i]) return cdf   def chi2_inv(p,x,k): cdf = chi2_cdf(x,k) ix = np.argmin(np.abs(cdf-p)) return x[ix]

   def chi2_pdf(x,k): bin = x[1]-x[0] fx = (x**((k/2)-1)*np.exp(-x/2))/(2**(k/2)) fx = fx/(bin*np.sum(fx)) return fx, bin def chi2_cdf(x,k): pdf,bin = chi2_pdf(x,k) cdf = np.zeros(len(x)) for i in range(len(x)): cdf[i] = bin*np.sum(pdf[:i]) return cdf def chi2_inv(p,x,k): cdf = chi2_cdf(x,k) ix = np.argmin(np.abs(cdf-p)) return x[ix]

   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...

   Reply

3. 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?

   Reply

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