How to calculate quantiles for a gamma distribution?

I would like to compute quantiles for a gamma distribution. I found a purported example here given as

$$\text{quantile}(a, b, p) = \frac{\gamma^{-1}(a, \Gamma(a) p )}{b}$$

where $$\gamma^{-1}$$ is the inverse of the lower incomplete gamma function and $$\Gamma$$ is the gamma function.

Trying this out in SciPy I found that for many values I got nan, suggesting something numerical going wrong.

from scipy.special import (gammaincinv, gamma as gamma_function)

def gamma_quantile_function(a, b, p):

return gammaincinv(a, gamma_function(a) * p) / b

gamma_quantile_function(3.8756542398707046, 5349.756221, 0.99) # returns nan


Is there just a problem with my implementation, or is the math wrong?

• I'm not familiar with python, but R (?qgamma) uses the algorithm proposed by Best and Roberts (1975) doi.org/10.2307/2347113 Jan 29 at 19:46
• Your method seems reasonable, so it looks like a Python code issue. In R both qgamma(p=0.99,shape=3.8756542398707046,rate=5349.756221) and library(zipfR); Igamma.inv(3.8756542398707046,gamma(3.8756542398707046)*0.99)/5349.756221 give 0.001840512 Jan 29 at 19:50

There are two Gamma incomplete functions: the regularized one and the non-regularized one. The Python function gammaincinv is the inverse of the regularized one, so you don't have to include the factor $$\Gamma(a)$$ in the second argument:

from scipy.special import gammaincinv

def gamma_quantile_function(a, b, p):

return gammaincinv(a, p) / b

gamma_quantile_function(3.8756542398707046, 5349.756221, 0.99)


This returns 0.001840512168532855, which coincides with the output of the R command qgamma(0.99, 3.8756542398707046, 5349.756221).

• It seems I suffered from some inattentional blindness while reading the docs. Thank you for noticing and pointing out the regularization. Jan 29 at 20:04