3
$\begingroup$

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?

$\endgroup$
2
  • 1
    $\begingroup$ I'm not familiar with python, but R (?qgamma) uses the algorithm proposed by Best and Roberts (1975) doi.org/10.2307/2347113 $\endgroup$
    – utobi
    Jan 29 at 19:46
  • $\begingroup$ 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 $\endgroup$
    – Henry
    Jan 29 at 19:50

1 Answer 1

5
$\begingroup$

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.