[SciPy-User] discrepancy in chi square cumulative distribution and gamma functions?

[James Jensen]

Thanks, Chad.

My reason for looking into this was that I am working with cdf values that
are very close to zero or one, such that chi2.cdf rounds them to 0.0 or
1.0. I hope that by working in log space I could distinguish between them.
So to get the log of the cdf, I was considering using something along the
lines of

ln(lower_incomplete_gamma(k/2, x/2)) - ln(gamma(k/2))

and scipy.special has a gammaln function.

I then realized that scipy.stats.chi2 has a logcdf function. But it doesn't
appear to have the additional precision I'm looking for. Is it possible
that it's implemented first computing the cdf, and then taking the
logarithm of it? If so, the precision would be lost prior to taking the
logarithm.


On Wed, Mar 25, 2015 at 4:43 PM, James Jensen <jdjensen <at> eng.ucsd.edu>
wrote:
> From what I understand, the cdf of a chi-square distribution with k
degrees
> of freedom and a test statistic x is
>
> lower_incomplete_gamma(k/2, x/2) / gamma(k/2)
>
> So you would expect the following two lines of code to yield the same
> result:
>
> cdf1 = scipy.stats.chi2.cdf(x, k)
> cdf2 = scipy.special.gammainc(0.5 * k, 0.5 * x) / scipy.special.gamma(0.5
*
> k)
>

It appears to be simply a naming/expectation issue. What is being
calculated with scipy.special.gammainc is not what e.g. Wikipedia
calls the lower incomplete gamma function but instead is what it calls
the regularized gamma function, see [1]. MATLAB's naming is consistent
with Scipy's (not with Wikipedia's), see [2].

In fact, both scipy.stats.chi2.cdf and scipy.special.gammainc call the
exact same underlying C function, it's just that scipy.stats.chi2.cdf
divides the arguments by 2 first.

You'll notice that the discrepancy-ratios you give are just equal to
scipy.special.gamma(0.5 * k), k = 1, 2, ..., 20.

[1]
http://en.wikipedia.org/wiki/Incomplete_gamma_function#Regularized_Gamma_functions_and_Poisson_random_variables
[2] http://www.mathworks.com/help/matlab/ref/gammainc.html

Cheers,
Chad
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20150326/7fca0ea9/attachment.html>