## Gamma Function in Scipy.Special

I was playing with the Gamma Distribution Function in python. As typically, I used a combination of ~~bumpy~~ numpy arrays and the scipy.special functions for Gamma and incomplete gamma functions. At one point I realized that the Gamma distribution function (bottom panel on the figure below) was not reaching unity for large x values and for some parameters (k=0.63, theta=0.05). In the given case, the deviation was about 0.3 in probability space, which is significant.

Note: The Gamma Distribution Function (its cdf or pdf) is not the same thing as the Gamma function!

The solution I came up with is to use the module mpmath. When I do, everything looks fine.

I googled only briefly, but did not find any description of this error in scipy.special. Also I am not sure, if using mpmath is the best/fastest/most elegant solution. However, at least it seems to work.

I get the same strange plots if I define the cdf of the gamma distribution as:

gamma_cdf = lambda x, alpha, beta: sp.special.gammainc(alpha, beta * x) / sp.special.gamma(alpha)

or, if you like the other formulation better:

gamma2 = lambda x, k, beta: gamma_cdf(x, k, 1 / beta)

which is the obvious thing to do. However in scipy.special, gammainc is defined as:

“1 / gamma(a) * integral(exp(-t) * t**(a-1), t=0..x)”

So, dividing by sp.special.gamma(alpha) is the straw that broke the gamma’s back. This will bring you up to 1:

gamma_cdf = lambda x, alpha, beta: sp.special.gammainc(alpha, beta * x)

By the way: are “bumpy arrays” more uneven than numpy ones? SCNR (I am also sorry for the url)

well, gamma gamma, thanks for pointing out what breaks ones neck.. 😉

Things may get even hairier when you use more complex distributions in scipy.stats.

One can get quite interesting insights into some distributions – after the headaches of transforming parameters from textbooks or wikipedia into loc, scale, etc. of the stats functions to achieve the same results, have gone away.

I think scipy.stats.gamma is one of those candidates.

By the way, why don’t you use that one?

thanks for pointing towards the functions integrated into scipy.stats — I guess I wanted to see things more directly… 😉