Incomplete Gamma function in scipy

Question

I would like to compute what wolfram alpha calls the incomplete gamma function (see here):

`gamma[0, 0.1]`

The wolfram alpha output is 1.822. The only thing scipy gives me that resembles this is scipy.special.gammainc, but it has a different definition than how wolfram alpha defines their incomplete gamma function.

Not surprisingly

import scipy
scipy.special.gammainc(0, 0.1)

gives me nan. Does scipy support what I'm looking for?

Answer 1

According to http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.special.gammainc.html, the first argument must be positive, whereas you have zero; that's why you're getting NaN.

That said, suppose we try to compute Gamma[0.01,0.1] instead. In this case WolframAlpha returns 1.80324:

According to http://mathworld.wolfram.com/IncompleteGammaFunction.html, this is the Upper Incomplete Gamma Function, whereas what Scipy outputs is a scaled version of what WolframAlpha calls the Lower Incomplete Gamma Function. By using the identity in Equation 10, one can see that in cases where a>0, you can use the following:

from scipy.special import gammainc
from scipy.special import gamma
gamma(0.01)*(1 - gammainc(0.01,0.1))

which returns 1.8032413569025461 in agreement with WolframAlpha.

In short, Gamma[a,x] in WolframAlpha corresponds to gamma(a)*(1-gammainc(a,x)) in Scipy, provided that a>0.

Answer 2

Unfortunately scipys gammaincc does not support a value of a=0. But you could define your own inc_gamma:

from scipy.special import gamma, gammaincc, exp1
def inc_gamma(a, x):
    return exp1(x) if a == 0 else gamma(a)*gammaincc(a, x)

Then inc_gamma(0, 0.1) will give you 1.8229239584193906.