Approximating the
Incomplete Gamma Function
Daniel D. Warner
July 25, 1999
The Incomplete Gamma Function and its complement are defined by
, and
.
From the definition of the Gamma Function we have the identity
.
The Incomplete Gamma Function is of particular interest in approximating the Cummulative Poisson Probability Function and the Cummualtive
Probability Function.
and
.
The Incomplete Gamma Function is also useful in approximating the Error Function and the Exponential Integrals
and
.
Except for the Exponential Integrals, these applications are restricted to
and
. Although the Incomplete Gamma Function can be defined for nonpositive
and negative
, we will not consider that situation in this study.
There is a rich collection of material on approximations for the Incomplete Gamma Function, starting with the material in Abramowitz and Stegun and ending with the sophisticated algorithm published by Walter Gautschi in 1979. Both Gautschi and the authors of the Numerical Recipes book employ a series expansion and a continued fraction expansion. In the Numerical Recipes approach the
is calculated using a series when
, and
is calculated using a continued fraction when
. Gautschi's algorithm takes a similar approach but there are more regions and the shape of the regions is more subtle. In part, the added complexity is caused by the fact that he is developing an approximation that willl work for nonpositive values of
. However, some of the complexity comes from avoiding unnecessary loss of precision through cancillation. That will be a concern in the approximation developed here.
>
IG := (a,x) -> (1/GAMMA(a)) * int(exp(-t)*t^(a-1),t = x .. infinity);
Ig := (a,x) -> (1/GAMMA(a)) * int(exp(-t)*t^(a-1),t = 0..x);
> plot({IG(0.5,x),IG(1.0,x),IG(3.0,x),IG(10.0,x)},x=0..15);
> plot({Ig(0.5,x),Ig(1.0,x),Ig(3.0,x),Ig(10.0,x)},x=0..15);
The following continued fraction, due to Legendre, will converge for positive
and for arbitrary real
.
The approach in the Numerical Recipes is to compute this continued fraction using a general continued fraction algorithm. However, Gautschi introduces two transformations which reduce the calculation to the calculation of a series, whose convergence is easier to verify.
This yields the approximation
,
where
, and
,
where
, and
,
and where
=
.
The following calculations serve to verify the transformations.
>
t[0] := 1:
s[0] := 1:
z[0] := t[0]/(x+1-a);
a[1] := 1*(a-1)/((x+2*1-a)^2-1):
s[1] := 1/(1+a[1]*s[0]):
t[1] := t[0]*(s[1]-1):
z[1] := normal((t[0]+t[1])/(x+1-a));
a[2] := 2*(a-2)/((x+2*2-a)^2-1):
s[2] := 1/(1+a[2]*s[1]):
t[2] := t[1]*(s[2]-1):
z[2] := normal((t[0] + t[1] + t[2])/(x+1-a));
>
z[0] := normal(1/(x+(1-a)/(1)));
z[1] := normal(1/(x+(1-a)/(1+1/(x+(2-a)/(1)))));
z[2] := normal(1/(x+(1-a)/(1+1/(x+(2-a)/(1+2/(x+(3-a)/(1)))))));
And the following calculations verify the actual algorithm.
>
p0 := 0:
q0 := (x-1-a)*(x+1-a):
r0 := 4*(x+1-a):
s0 := -a+1:
u0 := 0:
t0 := 1:
z0 := t0/(x+1-a);
p1 := p0 + s0:
q1 := simplify(q0 + r0):
r1 := r0 + 8:
s1 := s0 + 2:
v := p1*(1+u0):
u1 := simplify(v/(q1-v)):
t1 := u1*t0:
z1 := normal((t0+t1)/(x+1-a));
p2 := p1 + s1:
q2 := q1 + r1:
r2 := r1 + 8:
s2 := s1 + 2:
v := p2*(1+u1):
u2 := v/(q2-v):
t2 := u2*t1:
z2 := normal((t0+t1+t2)/(x+1-a));
Now consider
,
and
, then
will be positive and
. This means that the series will start off as an alternating series whose terms have monotone decreasing magnitudes. On the other hand, once
is larger than
then
can grow, but since
, we can show that
remains between 1 and 2.
>
IGCF := proc( aa,xx )
local x, a, tol, p, q, r, s, u, v, t, w, z, k;
x := evalf( xx );
a := evalf( aa );
tol := evalf(10^(-15));
p := 0;
q := (x-1.0-a)*(x+1.0-a);
r := 4.0*(x+1.0-a);
s := 1.0-a;
u := 0;
t := 1;
z := t;
for k from 1 to 100 do
w := t;
p := p + s;
q := q + r;
r := r + 8.0;
s := s + 2.0;
v := p*(1.0+u);
u := v/(q-v);
t := u*t;
z := z+t;
if ((0<w) and (0<u) and (w < tol*(1-u))) then
break
fi;
od;
RETURN( evalf( z/(x+1-a) ) );
end:
>
IgS := proc( aa, xx )
local a, x, tol, t, z, s, r, k;
x := eval( xx );
a := eval( aa );
tol := evalf(10^(-15));
t := 1;
z := t;
for k from 1 to 100 do
s := t;
r := x/(a+k);
t := t*r;
z := z+t;
if (s < tol*(1-r)) then
break
fi;
od;
RETURN( evalf( z ) );
end:
>
IG1 := (a,x) -> `if`(x<(a+1), evalf(1-exp(-x+a*log(x)-lnGAMMA(a+1))*IgS(a,x)), evalf(exp(-x+a*log(x)-lnGAMMA(a))*IGCF(a,x))):
Ig1 := (a,x) -> `if`(x<(a+1), evalf(exp(-x+a*log(x)-lnGAMMA(a+1))*IgS(a,x)), evalf(1-exp(-x+a*log(x)-lnGAMMA(a))*IGCF(a,x))):
> plot({IG1(0.5,x),IG1(1.0,x),IG1(3.0,x),IG1(10.0,x)},x=0..15);
> plot({Ig1(0.5,x),Ig1(1.0,x),Ig1(3.0,x),Ig1(10.0,x)},x=0..15);
> Digits := 20;
> plot3d((IG - IG1), 0..10, 0..15, title="Plot of the Absolute Errors",axes=BOX);
> plot3d((IG - IG1)/IG, 0..10, 0..15, title="Plot of the Relative Errors",axes=BOX);
>
readlib(C);
C(IGCF,ansi);
#include <math.h>
double IGCF(double aa,double xx)
{
double x;
double a;
double tol;
double p;
double q;
double r;
double s;
double u;
double v;
double t;
int w;
double z;
int k;
{
x = 0.1E1*xx;
a = 0.1E1*aa;
tol = 0.1E-14;
p = 0.0;
q = (x-0.1E1-a)*(x+0.1E1-a);
r = 0.4E1*x+0.4E1-0.4E1*a;
s = 0.1E1-a;
u = 0.0;
t = 1.0;
z = t;
for(k = 1;k <= 100.0;k++)
{
w = t;
p += s;
q += r;
r += 0.8E1;
s += 0.2E1;
v = p*(0.1E1+u);
u = v/(q-v);
t = u*t;
z += t;
if( 0.0 < w && 0.0 < u && w < tol*(1.0-u) )
break;
}
return(0.1E1*z/(x+1.0-a));
}
}
> C(IgS,ansi);
#include <math.h>
double IgS(double aa,double xx)
{
double a;
double x;
double tol;
double t;
double z;
int s;
double r;
int k;
{
x = (* xx);
a = (* aa);
tol = 0.1E-14;
t = 1.0;
z = t;
for(k = 1;k <= 100.0;k++)
{
s = t;
r = x/(a+k);
t = t*r;
z += t;
if( s < tol*(1.0-r) )
break;
}
return(0.1E1*z);
}
}
>
References:
Abramowitz, M. and Stegun, I.A.,
Handbook of Mathematical Functions
,
Applied Mathematics Series, Vol. 55, National Bureau of Standards, 1964.
reprinted by Dover Publications, 1968.
Walter Gautschi,
"A Computational Procedure for Incomplete Gamma Functions",
ACM Transactions on Mathematical Software
,
Vol. 5, No. 4, December 1979, pp. 466-481.
Luke, Y. L.,
The Special Functions and Their Approximations
,
Volumes I and II, Academic Press, 1969.
William H. Press, Brian P. Flannery, Saul A. Teukolsky, and Wiliam T. Vetterling,
Numericl Recipes, The Art of Scientific Computing
,
Cambridge University Press. 1986.