incomplete gamma function in python?

Question

the scipy.special.gammainc can not take negative values for the first argument. Are there any other implementations that could in python? I can do a manual integration for sure but I'd like to know if there are good alternatives that already exist.

Correct result: 1 - Gamma[-1,1] = 0.85

Use Scipy: scipy.special.gammainc(-1, 1) = 0

Thanks.

score 9 · Accepted Answer · answered May 10 '12 at 22:27

9

I typically reach for mpmath whenever I need special functions and I'm not too concerned about performance. (Although its performance in many cases is pretty good anyway.)

For example:

>>> import mpmath
>>> mpmath.gammainc(-1,1)
mpf('0.14849550677592205')
>>> 1-mpmath.gammainc(-1,1)
mpf('0.85150449322407795')
>>> mpmath.mp.dps = 50 # arbitrary precision!
>>> 1-mpmath.gammainc(-1,1)
mpf('0.85150449322407795208164000529866078158523616237514084')

answered May 10 '12 at 22:27

DSM

342,061
65
592
494

thanks, I did a brief search online but failed to find a standalone solution. Guess this is one of the few occasions that I install an extra package for a single functionality ;-) – nye17 May 11 '12 at 01:30
4

what does the 'mpf' stand for? – LuizAngioletti Jan 04 '14 at 14:12

score 2 · Answer 2 · answered Feb 28 '13 at 10:24

I just had the same issue and ended up using the recurrence relations for the function when a<0. http://en.wikipedia.org/wiki/Incomplete_gamma_function#Properties

Note also that the scipy functions gammainc and gammaincc give the regularized forms Gamma(a,x)/Gamma(a)

score 1 · Answer 3 · answered Oct 28 '21 at 22:14

Still an issue in 2021, and they still haven't improved this in scipy. Especially it is frustrating that scipy does not even provide unregularised versions of the upper and lower incomplete Gamma functions. I also ended up using mpmath, which uses its own data type (here mpf for mpmath floating - which supports arbitrary precision). In order to cook up something quick for the upper and lower incomplete Gamma function that works with numpy arrays, and that behaves like one would expect from evaluating those integrals I came up with the following:

import numpy as np
from mpmath import gammainc

"""
In both functinos below a is a float and z is a numpy.array.
"""

def gammainc_up(a,z):
    return np.asarray([gammainc(a, zi, regularized=False)
                       for zi in z]).astype(float)

def gammainc_low(a,z):
    return np.asarray([gamainc(a, 0, zi, regularized=False)
                       for zi in z]).astype(float)

Note again, this is for the un-regularised functions (Eq. 8.2.1 and 8.2.2 in the DLMF), the regularised functions (Eq. 8.2.3 and 8.2.4) can be obtined in mpmath by setting the keyword regularized=True.

incomplete gamma function in python?

3 Answers3

Linked