GammaCHI: a package for the inversion and computation of the gamma and chi-square cumulative distribution functions (central and noncentral)

0.21 MB·English

by  A. Gil

GammaCHI: a package for the inversion and
computation of the gamma and chi-square cumulative
distribution functions (central and noncentral)
5
Amparo Gila, Javier Segurab, Nico M. Temmec
1
0
a
2 Depto. de Matem´atica Aplicada y Ciencias de la Comput. Universidad de Cantabria.
39005-Santander, Spain. e-mail: [email protected]
n
b
a Depto. de Matem´aticas, Estad´ıstica y Comput. Universidad de Cantabria.
J 39005-Santander, Spain
7 c IAA, 1391 VD 18, Abcoude, The Netherlands1

]
S
M
. Abstract
s
c
[ A Fortran 90 module GammaCHI for computing and inverting the gamma

1 and chi-square cumulative distribution functions (central and noncentral) is
v presented. The main novelty of this package are the reliable and accurate
8
7 inversion routines for the noncentral cumulative distribution functions. Ad-
5 ditionally, the package also provides routines for computing the gamma func-
1
tion, the error function and other functions related to the gamma function.
0
. The module includes the routines cdfgamC, invcdfgamC, cdfgamNC,
1
0 invcdfgamNC, errorfunction, inverfc, gamma, loggam, gamstar and
5
quotgamm for the computation of the central gamma distribution function
1
: (and its complementary function), the inversion of the central gamma dis-
v
i tribution function, the computation of the noncentral gamma distribution
X
function (and its complementary function), the inversion of the noncentral
r
a gamma distribution function, the computation of the error function and its
complementary function, the inversion of the complementary error function,
the computation of: the gamma function, the logarithm of the gamma func-
tion, the regulated gamma function and the ratio of two gamma functions,
respectively.
PROGRAM SUMMARY
1
Former address: CWI, 1098 XG Amsterdam, The Netherlands
Preprint submitted to Computer Physics Communications January 8, 2015
Manuscript Title:
GammaCHI: a package for the inversion and computation of the gamma and chi-
square cumulative distribution functions (central and noncentral)
Authors:
Amparo Gil, Javier Segura, Nico M. Temme
Program Title:
Module GammaCHI
Journal Reference:
Catalogue identifier:
Licensing provisions:
Programming language:
Fortran 90
Computer:
Any supporting a FORTRAN compiler.
Operating system:
Any supporting a FORTRAN compiler.
RAM:
a few MB
Number of processors used:
Keywords:
Gamma cumulative distribution function; chi-square cumulative distribution func-
tion; inversion of cumulative distribution functions; error function; complemen-
tary error function; gamma function; logarithm of the gamma function; regulated
gamma function; quotient of gamma functions.
Classification:
4.7 Other functions.
External routines/libraries:
None.
Subprograms used:
Nature of problem:
The computation and inversion of gamma and chi-square cumulative distribu-
tion functions (central and noncentral) as well as the computation of the error
and gamma functions is needed in many problems of applied and mathematical
physics.
Solution method:
The algorithms use different methods of computation depending on the range of
parameters: asymptotic expansions, quadrature methods, etc.
Restrictions:
In the inversion of the central gamma/chi-square distribution functions, very
2
small input function values P (x,y), Q (x,y) (lower than 10−150) are not admis-
µ µ
sible.
The admissible input parameter ranges for computing the noncentral cumu-
lative gamma distribution functions P (x,y), Q (x,y) in standard IEEE double
µ µ
precision arithmetic are 0 x 10000, 0 y 10000, 0.5 µ 10000 and
≤ ≤ ≤ ≤ ≤ ≤
the related parameter ranges for the noncentral chi-square cumulative distribution
function. In the inversion of the noncentral gamma/chi-square distribution func-
tions, very small input function values P (x,y), Q (x,y) (lower than 10−25 and
µ µ
10−35, respectively) are not admissible.
Running time:
It varies depending on the function and the parameter range.
1. Introduction
The Fortran 90 module GammaCHI provides reliable and fast routines
for the inversion and computation of the gamma and chi-square distribution
functions. These functions appear in many problems of applied probability
including, of course, a large number of problems in Physics. The module
also includes routines for the computation of the gamma function, the error
function and its complementary function, the inverse of the complementary
error function, the logarithm of the gamma function, the regulated gamma
function and the ratio of two gamma functions.
The main novelty of the package GammaCHI are the inversion routines
for the noncentral cumulative distribution functions. These algorithms are
based on the asymptotic methods presented in [1] although improvements in
the estimation of the initial values for small values of the parameter µ are
included in the present algorithms. Among other applications, the compu-
tation and the inversion of the noncentral cumulative distribution functions
included in GammaCHI appear in the analysis of signal detection in differ-
ent physical scenarios such as optics [2], radiometry [3] or quantum detection
[4], [5]. Also, the inversion routines can be used to generate normal variates
as well as central and noncentral gamma (or chi-square) variates, which can
be useful, for example, in Monte Carlo simulations.
2. Theoretical background
2.1. Central gamma distribution function
The incomplete gamma functions are defined by
3
x ∞
γ(a,x) = ta−1e−tdt, Γ(a,x) = ta−1e−tdt. (1)
Z Z
0 x
Let us define
1 1
P(a,x) = γ(a,x), Q(a,x) = Γ(a,x), (2)
Γ(a) Γ(a)
where we assume that a and x are real positive numbers.
The functions P(a,x), Q(a,x) are the central gamma distribution func-
tion and its complementary function, respectively. These distribution func-
tionswhichappearinmanyproblemsofappliedprobabilityinclude, aspartic-
ularcases, thestandardchi-squareprobabilityfunctionsP(χ2 ν)andQ(χ2 ν)
| |
with parameters a = ν/2 and x = χ2/2.
The central gamma distribution functions satisfy the complementary re-
lation
P(a,x)+Q(a,x) = 1. (3)
In hypothesis testing, and using the chi-square distribution, it is usual to
consider the problem as Q(χ2 ν) = α; α is the probability that a variable
|
distributed according tothechi-square distributionwithν degrees offreedom
exceeds the value χ2. Two kinds of problems are usually considered: a)
computing the confidence level α given an experimentally determined chi-
square value χ2; b) computing percentage points for certain values of α and
for different degrees of freedom ν (the standard tables which are common
in statistics text-books usually include this information). The first problem
involves the direct computation of the cumulative distribution function while
the second one is an inversion problem.
The module GammaCHI includes routines both for the direct compu-
tation of the distribution functions (central gamma/chi-square) and their
inversion. Regarding the direct computation, the routine first computes
min P(a,x),Q(a,x) ,andtheotheronebyusing(3). ComputingtheQ(a,x)
{ }
function simply as 1 P(a,x) when P(a,x) is close to 1 can lead to serious
−
cancellation problems, which are avoided in our algorithms.
On the other hand, the inversion routine solves the equations
P(a,x) = p, Q(a,x) = q, 0 < p, q < 1, (4)
for a given value of a.
4
The algorithm for computing the gamma distribution function and its
inversion is described in [6]. The computation of the distribution functions
is based on the use of Taylor expansions, continued fractions or uniform
asymptotic expansions, depending on the parameter values. For the inver-
sion, asymptoticexpansions incombinationwithhigh-orderNewtonmethods
are used.
The version included in the package GammaCHI for the inversion in-
cludes some improvements in its performance for small values of p and q.
For such small values, some of the elementary bounds for incomplete gamma
functions given in [7] provide accurate enough approximations for estimating
sufficiently accurate starting values for the Newton iterative process.
2.2. Noncentral gamma distribution function
The noncentral gamma cumulative distribution functions can be defined
as
∞
xk
P (x,y) = e−x P(µ+k,y),
µ
k!
X
k=0 (5)
∞
xk
Q (x,y) = e−x Q(µ+k,y),
µ
k!
X
k=0
where P(µ,y) and Q(µ,y) are the central gamma cumulative distribution
function and its complementary function, respectively; x is called the non-
centrality parameter of the distribution functions.
As in the case of the central distribution functions, the noncentral func-
tions satisfy the relation
P (x,y)+Q (x,y) = 1. (6)
µ µ
On the other hand, if χ2(λ) is a random variable with a noncentral chi-
n
square distribution with n > 0 degrees of freedom and noncentrality param-
eter λ 0, then the lower and upper tail probabilities for the distribution
≥
function of χ2(λ) are given by
n
∞ (λ/2)k
Prob(χ2(λ) < t) = e−λ/2 P(n/2+k, t/2),
n k!
X
k=0 (7)
∞ (λ/2)k
Prob(χ2(λ) > t) = e−λ/2 Q(n/2+k, t/2),
n k!
X
k=0
5
from which, and comparing with (5), it follows the simple relation between
distribution functions: if the random variable Y has a noncentral gamma
distribution with parameters a and noncentrality parameter λ, then X = 2Y
has a noncentral chi-square distribution with parameter n = 2a and with
noncentrality parameter 2λ.
Algorithms for computing the functions P (x,y) and Q (x,y) for a large
µ µ
range of the parameters µ (µ 1), x, y are described in [8]. The methods
≥
include series expansions in terms of the incomplete gamma functions, recur-
rence relations, continued fractions, asymptotic expansions, and numerical
quadrature. In the module GammaCHI the lower µ range of computation
of the algorithm for computing the noncentral gamma cumulative distribu-
tionfunctionisextended toinclude theinterval 1 µ < 1. Forextending the
2 ≤
algorithms to this interval, we use a single step of the following three-term
homogeneous recurrence relation (Eq.14 of [8])
y I (2√xy)
µ
y (1+c )y +c y = 0, c = . (8)
µ+1 µ µ µ µ−1 µ
− rxI (2√xy)
µ−1
Both Q (x,y) and P (x,y) satisfy (8). In this expression, ratios of Bessel
µ µ
functions appear. These ratios are efficiently computed using a continued
fraction representation.
Regarding the inversion of the noncentral gamma/chi-square distribution
functions, two different kinds of inversion problems can be considered:
Problem 1: Find x from the equation
Q (x,y) = q, (9)
µ
with fixed y, µ and q.
Problem 2: Find y from the equation
P (x,y) = p, (10)
µ
with fixed x, µ and p.
The inversion of Q (x,y) with respect to x corresponds to the problem
µ
of inverting the distribution function with respect to the noncentrality pa-
rameter given the upper tail probability. On the other hand, the inversion
of P (x,y) with respect to y with fixed x corresponds to the problem of
µ
6
computing the p-quantiles of the distribution function. For noncentral chi-
square distributions, the inversion with respect to y allows the computation
of a variate with a given number of degrees of freedom and a given noncen-
trality parameter.
In the inversion with respect to x, it is important to note that not for all
the possible values of the pair (y ,q ) there is a solution of the equation
0 1
Q (x,y ) = q . (11)
µ 0 1
In particular, the values of q such that q < q , where q = Q (0,y ) =
1 1 0 0 µ 0
Q (y ), do not make sense because Q (x,y) is increasing as a function of x.
µ 0 µ
Both in the inversion with respect to x and y, a combination of asymp-
totic expansions and secant methods can be used, as discussed in [1]. A
Fortran 90 implementation is now available in the module GammaCHI. As
for the direct computation of function distribution values, in this implemen-
tation the range of computation for the inversion of the gamma cumulative
distribution function is extended to include the interval 1 µ < 1. The
2 ≤
computation in this interval is made by using the double asymptotic prop-
erty of the expansions derived in [1]: the formulas obtained in [1] can all be
used for small values of µ if ξ = 2√xy is large. More details can be found in
the appendix. For small q this is suplemented with estimations for the case
µ = 1/2, which can be written in terms of error functions; we use the result
of the inversion for the case µ = 1/2 to estimate a starting value when µ is
close to 1/2. It is easy to check that
1
Q (x,y) = erfc(√x+√y)+erfc(√y √x)
1/2
2 −
(cid:0) (cid:1)
and for this function, computing derivatives is simple and one can apply the
fourth order method described in [9].
2.3. Error function and its complementary function
The error function erf(x) [10] is defined by the integral
2 x
erf(x) = e−t2 dt. (12)
√π Z
0
Its complement with respect to 1 is denoted by
2 ∞
erfc(x) = 1 erf(x) = e−t2 dt. (13)
− √π Z
x
7
The complementary error function is very closely related to the normal
distribution functions, which are defined by
1 x 1 ∞
P(x) = e−t2/2dt, Q(x) = e−t2/2dt (14)
√2π Z √2π Z
−∞ x
with the property P(x)+Q(x) = 1.
From the definition of the complementary error function it follows that
1 1
P(x) = erfc( x/√2), Q(x) = erfc(x/√2). (15)
2 − 2
The algorithm for the error functions is based on the use of rational ap-
proximations [11] andit includes thepossibility of computing thescaled com-
plementary error function ex2erfc(x). This is important to avoid underflow
problems in numerical computations of the complementary error function
(the dominant exponential behavior of the function is e−x2). In general, it is
quite useful (when possible) todefine scaled valuesof mathematical functions
with the dominant exponential behavior factored out [12, 12.1.3].
§
Our module includes also the computation of the inverse of the comple-
mentary error function y = inverfcx. Using the relation (15), this routine
can be also useful in the generation of normal variates.
For values of x close to 1, the following expansion is used
inverfcx = t+ 1t3 + 7 t5 + 127t7 +..., 0 < x < 2. (16)
3 30 630
with t = 1√π(1 x). For these and more coefficients, see [13].
2 −
2.4. Gamma function
Three standard definitions of the gamma function are usually considered:
∞
Γ(z) = e−ttz−1dt, Re z > 0 (Euler). (17)
Z
0
n!nz
Γ(z) = lim , z = 0, 1, 2,.... (Euler). (18)
n→∞ z(z +1) (z +n) 6 − −
···
∞
1 z
= zeγz 1+ e−z/n (Weierstrass), (19)
Γ(z) n
nY=1(cid:16) (cid:17)
with γ = 0.57721..., Euler’s constant. The equivalence of these definitions is
proved in [14].
8
Important relations are
Γ(z +1) = zΓ(z), (recursion), (20)
and
1 sinπz
= , (21)
Γ(1+z)Γ(1 z) πz
−
which is called the reflection formula. This formula is important for range
reduction and it is easily proved by using (19).
For computing the gamma function, we use recursion for small values of
x and the relation given in (22) for larger values. The special cases x = n
and x = n+ 1 are treated separately.
2
2.4.1. Logarithm of the gamma function
The following representation of logΓ(z) is very important for deriving
expansions for large values of z :
| |
logΓ(z) = z 1 logz z + 1 log(2π)+S(z). (22)
− 2 − 2
(cid:16) (cid:17)
For large z , S(z) gives a small correction with respect to the remaining
| |
terms of the right-hand side. Neglecting S(z) gives the well-known Stirling
formula
1/2
2π
Γ(z) zze−z , z . (23)
∼ (cid:18) z (cid:19) → ∞
S(z)canbewrittenintermsofanexpansion involving Bernoulli numbers:
N−1
B
2n+2
S(z) = +E (z). (24)
(2n+1)(2n+2)z2n+1 N
Xn=0
The Bernoulli numbers B are given by B = B (0). The first few are
n n n
1 1 1
B = 1, B = , B = , B = 0, B = . More information on
0 1 −2 2 6 3 4 −30
Bernoulli polynomials and Bernoulli numbers can be found in [15, 24.2(i)].
§
A bound for E in (24) can be obtained as
N
| |
B K(z)
E 2N+2 z −2N−1, (25)
N
| | ≤ (2N +1)(2N +2)| |
9
where
z2
K(z) = sup .
(cid:12)z2 +u2(cid:12)
u≥0(cid:12) (cid:12)
(cid:12) (cid:12)
If argz < π/4, then K(z) = 1. (cid:12)For real(cid:12)positive z, E (z) is less in
N
| |
absolute value than the first term neglected in (24) and it has the same sign.
From (25) it follows that E (z) = (z−2N−1) for Re z + . Inserting the
N
O → ∞
valuesofthefirstBernoullinumberswearriveattherepresentation(Stirling’s
series):
1 1
logΓ(z)= z logz z + log(2π) +
− 2 − 2
(cid:16) (cid:17) (26)
1 1 + 1 1 + (z−9).
12z − 360z3 1260z5 − 1680z7 O
The algorithm for the logarithm of the gamma function makes use of the
representation (22). The computation of S(z) is based on the use of the sum
(22), Chebyshev series and a minimax rational approximation. Coefficients
of rational approximations are taken from the the tables in [16, 6.6]. On
§
some intervals we use Chebyshev expansions, with coefficients derived by
using high precision Maple codes.
2.5. The function Γ∗(x)
This function (the regulated gamma function) is defined by
Γ(x)
Γ∗(x) = , x > 0. (27)
2π/xxxe−x
p
When x is large, this function can be very important in algorithms where
the function Γ(x) is involved because Γ∗(x) = 1 + (1/x), as can be seen
O
from its asymptotic expansion (Stirling series):
Γ∗(x) 1+ 1 x−1 + 1 x−2 +..., x . (28)
∼ 12 288 → ∞
The algorithm for computing Γ∗(x), uses the relation Γ∗(x) = exp(S(x))
forx > 3,whereS(x)isgivenin(24). ThecomputationofS(x)wasdescribed
in the previous section.
10