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