The gammaincinv function
The inverse of the incomplete gamma function
Definition
Given two real values $y\in[0,1]$ and $a>0$, gammaincinv(y,a,tail) gives the
value $x$ such that
\[P(x,a) = y,\quad\text{ or }\quad Q(x,a) = y\]
depending on the value of tail ("lower" or "upper").
Computation
The computation is carried out via standard Newton method, solving, for $x$
\[ y - P(x,a) = 0 \quad\text{ or }\quad y - Q(x,a) = 0. \]
The Newton method uses as exit criterion a control on the residual and on the
maximum number of iterations.
For numerical stability, the initial guess and which gammainc function is
inverted depends on the values of $y$ and $a$. We use the following property:
\[\text{gammaincinv}(y,a,\text{`upper'})=\text{gammaincinv}(1-y,a,\text{`lower'})\]
which follows immediately by the property
\[P(x,a) + Q(x,a) = 1.\]
We start computing the trivial
values, i.e. $a=1,y=0$ and $y=1$. Then we rename the variables as follows (in order to obtain a more classical notation): if
tail = "lower", then $p=y,q=1-y$; while if tail = "upper", $q=y,p=1-y$. We will
invert with $p$ as $y$ in the cases in which we invert the lower incomplete gamma function, and with $q$ as $y$ in the cases in which we invert the upper one.
Now, the initial guesses (and the relative version of the incomplete gamma function to invert) are the following:
-
$p<\dfrac{(0.2(1+a))^a}{\Gamma(1+a)}$. We define
\[r = (p\,\Gamma(1+a))^{\frac{1}{a}}\]
and the parameters
\begin{align*}
c_2 &= \dfrac{1}{a+1}\\
c_3 &= \dfrac{3a+5}{2(a+1)^2(a+2)}\\
c_4 &= \dfrac{8a^2+33a+31}{3(a+1)^3(a+2)(a+3)}\\
c_5 &= \dfrac{125a^4+1179a^3+3971a^2+5661a+2888}{24(1+a)^4(a+2)^2(a+3)(a+4)}\\
\end{align*}
Then the initial guess is
\[x_0 = r + c_2\,r^2+c_3\,r^3+c_4\,r^4+c_5\,r^5\]
and we invert the lower incomplete gamma function.
-
$q<\dfrac{e^{-a/2}}{\Gamma(a+1)}, a\in(0,10)$. The initial guess is
\[x_0 = - \log (q) - \log(\Gamma(a))\]
and we invert the upper incomplete gamma function.
-
$a\in(0,1)$. The initial guess is
\[x_0 = (p\Gamma(a+1))^{\frac{1}{a}}\]
and we invert the lower incomplete gamma function.
-
Other cases. The initial guess is
\[x_0 = a\,t^3\]
where
\[t = 1-d-\sqrt{d}\,\text{norminv}(q),\quad d=\dfrac{1}{9a}\]
(where norminv is the inverse of the normal distribution function) and we
invert the upper incomplete gamma function.
The first three initial guesses can be found in the article "Efficient and
accurate algorithms for the computation and inversion of the incomplete gamma
function rations" by Gil, Segura ad Temme; while the last one can be found in
the documentation of the function igami.c of the cephes library.
You can find my work on my repository on bitbucket
here, bookmark "gammainc".
Comments
Working on this function I found some problems in gammainc. To solve them, I
just changed the flag for the condition 5 and added separate functions to
treat Inf cases.
