INTEGRAL REPRESENTATIONS FOR CONVOLUTIONS OF
NON–CENTRAL MULTIVARIATE GAMMA DISTRIBUTIONS
T. ROYEN
Fachhochschule Bingen, University of applied sciences,
Berlinstrasse 109, D–55411 Bingen, Germany
E–mail: royen@fh-bingen.de
Abstract. Three types of integral representations for the cumulative distribution
functions of convolutions of Γp(αk,Σk,∆k)–distributions with non–centrality matri-
ces ∆k are given by integration of products of simple complex functions over the
p–cube (−π, π]p. In particular, the joint distribution of the diagonal elements of a
generalized quadratic form XAX ′ with n independentNp(µk,Σ)–distributed columns
in Xp×n and a fixed A ≥ 0 is obtained. For a single Γp(α,Σ,∆)–cdf (p− 1)–variate
integrals over (−π, π]p−1 are derived. The integrals are numerically more favourable
than integrals obtained from the Fourier– or Laplace inversion formula.
Key words and phrases: Convolutions of multivariate distributions, generalized
quadratic forms of normal random vectors, multivariate chi–square distribution, Mul-
tivariate gamma distribution
AMS 2000 subject classifications: 62H10, 62E15
1. Introduction
The following notations are used:
(n) stands for
n1+...+np=n
n1, . . . , np ∈ N0 and
without any indices means
(n). The notation D ≥ 0
is also used for non–symmetrical matrices Dp×p with only non–negative eigenvalues.
The spectral norm of a p× p–matrix B is denoted by ‖B‖, I or Ip is always an identity
matrix and Cp is the p–cube (−π, π]p.
The Laplace transform (L.t.) of a p–variate non–central Γp(α,Σ,∆)–density with
α > 0, Σ > 0 and a non–centrality matrix ∆ ≥ 0 was originally obtained from the L.t.
of a non-central Wp(2α,Σ,∆)–Wishart distribution (with an additional scale factor 2)
and is given by
f̂(t1, . . . , tp;α,Σ,∆) = |Ip +ΣT |−αetr(−ΣT (I +ΣT )−1∆), (1)
T = diag(t1, . . . , tp), t1, . . . , tp ≥ 0.
This function f̂ is generally the L.t. of the density of a real measure on (0,∞)p which
is not always a probability measure. The term ”Γp(α,Σ,∆)–distribution” is used here
in this general sense. The exact set of values α, leading to a probability density (pdf)
f(x1, . . . , xp;α,Σ,∆), depends on Σ and presumedly on ∆. To obtain a pdf, all positive
integers 2α (degrees of freedom) are admissible and all 2α > p − 1. Moreover, in the
central case all non–integer values 2α > p− 2 ≥ 0 are allowed. For p− 2 < 2α < p− 1
http://arxiv.org/abs/0704.0539v1
see Royen (1997). Furthermore all α > 0 are admissible if |I + ΣT |−1 is infinitely
divisible. Two characterizations of infinite divisibility of a Γp(α,Σ)–distribution are
found in Griffiths (1984) and Bapat (1989). Further conditions for admissible non–
integer 2α < p− 2 are given in Royen (1997), (2006).
Three integral representations by integration over Cp are provided by theorem 2 in
section 4 for the functions
F (x1, . . . , xp;α1, . . . , αn,Σ1, . . . ,Σn,∆1, . . . ,∆n) (2)
. . .
f(ξ1, . . . , ξp;α1, . . . , αn,Σ1, . . . ,Σn,∆1, . . . ,∆n)dξ1 . . . dξp,
where f has the L.t.
|Ip +ΣkT |−αketr(−ΣkT (I +ΣkT )−1∆k), (3)
α1, . . . , αn > 0,Σ1, . . . ,Σn > 0,∆1, . . . ,∆n ≥ 0.
Thus, F is not always the cumulative distribution function (cdf) of a probability measure.
In particular letXp×n be aNp×n(Mp×n,Σp×p⊗In)–random matrix and An×n ≥ 0 of
rank q with T ′AT = Λ = diag(λ1, . . . , λn), λ1 ≥ . . . ≥ λn. Then the joint distribution of
the diagonal elements of the generalized quadratic form 1
XAX ′ equals the distribution
of the diagonal of 1
Y ΛY ′ with a Np×n(MT,Σp×p⊗In)–distributed Y = XT . This is the
distribution of a sum of q independent Γp(
, λkΣ,∆k =
−1)–random vectors,
where µ∗k is the k–th column of M
∗ = MT . This joint distribution of p quadratic
forms of normal random vectors is comprised within theorem 2 as a special case with
, Σk = λkΣ, k = 1, . . . , q. For methods under more general assumptions see also
Blacher (2003). For a survey of univariate quadratic forms of normal random variables
see chapter 4 in Mathai and Provost (1992). For several quadratic forms of skew elliptical
distributions see B.Q. Fang (2005).
In Royen (1991), (1992) three different types of series expansions for the χ2p(2α,Σ)–
cdf were derived from three different representations of the χ2p(2α,Σ)–L.t. which are
extended to the general Γp(α,Σ,∆)–L.t. in section 3 in a similar way as in Royen
(1995).
Some series expansions, closely related to the first two types, are already found in
Khatri, Krishnaiah and Sen (1977). The third type was introduced because of its superior
convergence properties. The simple method to transform many series expansions into
integrals over Cp is explained in more detail in section 2 and summarized in theorem 1.
The idea is as follows:
If A(z1, . . . , zp) and B(z1, . . . , zp) are analytical functions whose power series have
the coefficients a(m1, . . . ,mp) and b(n1, . . . , np) and which are absolutely convergent for
max |zj | < rA and max |zj | < rB respectively, where r−1B < rA, then
(2π)−p
A(y1, . . . , yp)B(y
1 , . . . , y
p )dϕ1 . . . dϕp (4)
a(n1, . . . , np)b(n1, . . . , np)
holds with yj = re
iϕj , −π < ϕj ≤ π, j = 1, . . . , p and r−1B < r < rA.
The integrals in (4) might be more economical than the series if the generating
functions A and B are simple available functions and if the series are slowly convergent
with very intricate coefficients. For non–central multivariate gamma distributions series
expansions are practically not feasible.
The integral representations in theorem 2 of section 4 are of the type in (4). As
long as no elementary density formulas are availale it should be a reasonable way to
obtain the joint cdf by integration of elementary terms only over Cp and not over Rp as
by the Fourier or Laplace inversion formula. A single Γp(α,Σ,∆)–cdf is represented by
a (p− 1)–variate integral over Cp−1 in section 5.
A totally different
–variate integral representation of the Γp(α,Σ)–cdf has
been given recently by Royen (2006), which is based on m–factorial decompositions∑
p×p = D −BB′, where D is a real or complex diagonal matrix minimizing the rank
m of Σ−1 −D. Approximations to a Γp(α,Σ)–cdf are obtained by m–factorial approx-
imations to Σ with a low value of m. These approximations are improved further by
successive correction terms.
2. The method
Theorem 1 in this section can be generalized in many ways, e.g. for Fourier trans-
forms, but the version below is sufficient for the purpose of the underlying paper.
Let f̂(t1, . . . , tp), t1, . . . , tp ≥ 0, be a given L.t. of an unknown function f(x1, . . . , xp)
with f = 0 for minxj < 0. It is assumed that there are univariate L.t. ĝj0(t) of
some probability densities gj0(x) on (0,∞) and further functions hj(t) with |hj(t)| ≤ 1,
uniformly for t ≥ 0, which enable a representation
f̂(t1, . . . , tp) =
ĝj0(tj)
B (h1(t1), . . . , hp(tp)) (5)
with an analytical function B(z1, . . . , zp) whose power series expansion
b(n1, . . . , np)
j (6)
is absolutely convergent for |z1|, . . . , |zp| < rB with a certain value rB > 1.
Furthermore, the products ĝj0(t)(hj(t))
n are supposed to be the L.t. of continuous
functions gjn(x), x > 0, which satisfy the conditions
|gjn(x)| ≤ nck(x) with a constant c and
(7)∫ ∞
k(x)e−txdx < ∞ for all t > 0 .
Hence, the generating functions (generators)
gj(x, y) =
gjn(x)y
n, j = 1, . . . , p, (8)
are defined for all x > 0 and |y| < 1, and they have the L.t.
ĝj(t, y) =
ĝj0(t)
1− yhj(t)
, t ≥ 0. (9)
Theorem 1. Under the assumptions from (6) and (7) f̂ in (5) is the L.t. of
f(x1, . . . , xp) = (2π)
B(y−11 , . . . , y
gj(xj , yj)dϕj (10)
with yj = re
iϕj , −π < ϕ ≤ π, r−1B < r < 1, gj from (8).
Proof. The integral in (10) is evaluated by
(2π)−p
b(m1, . . . ,mp)
gjnj (xj)y
 dϕ1 . . . dϕp
b(n1, . . . , np)
j=1 gjnj (xj)
and this series has the L.t. from (5).
Some further remarks: With
Gj(xj , yj) =
gj(ξ, yj)dξ (11)
instead of the gj in (10), the corresponding representation arises for
F (x1, . . . , xp) =
. . .
f(ξ1, . . . , ξp)dξ1 . . . dξp. (12)
If the series in (8) are absolutely convergent for all y ∈ C then additionally
hj(t) = 0 (13)
is supposed to hold. Then the rhs of (9) is the L.t. of gj(x, y) for any fixed y and all
sufficiently large t.
In some cases the functions gj0 and their L.t. ĝj0 are known from univariate marginal
distributions apart from some scale factors. If the functions uj = hj(t) are explicitly
invertible then
B(u1, . . . , up) =
h−11 (u1), . . . , h
p (up)
j=1 ĝj0
h−1j (uj)
) (14)
can sometimes be found easily from the given f̂ .
3. Three representations for the Γp(α,Σ,∆)–Laplace transform
and the related generators
With any v > 0 we define
zj = (1 + v
−1tj)
−1, tj ≥ 0, uj = 1− zj = v−1tjzj , ωj = zj − uj ,
Z = diag(z1, . . . , zp), U = diag(u1, . . . , up), Ω = diag(ω1, . . . , ωp).
The scale factor v is introduced to obtain ‖B‖ < 1 for the matrices B defined in (20)
below and to effect the convergence of some series expansions. For a more general scaling
see remarks following theorem 2 in section 4.
From the relations
v−1T = UZ−1, Ip = Z + U, Ω = Z − U, (16)
it follows for the matrices I +ΣT in the L.t. (1):
I +ΣT = I + vΣUZ−1 = (Z + vΣU)Z−1 (17)
Z + vΣU =
I + (vΣ− I)U, (18a)
vΣ(I + (v−1Σ−1 − I)Z), (18b)
(I + vΣ)(I + (2(I + vΣ)−1 − I)Ω), (18c)
and therefore
|I +ΣT |−α = cα|Z|α|I +BY |−α (19)
Y = U, B = vΣ− I, c = 1, (20a)
Y = Z, B = (vΣ)−1 − I, c = |I +B|, (20b)
Y = Ω, B = 2(I + vΣ)−1 − I, c = |I +B|. (20c)
It should be noticed that ‖B‖ < 1 in (20c) for every v > 0 and Σ > 0.
Now, using (16), by a straightforward calculation the L.t. in (1) can be represented
f̂(t1, . . . , tp;α,Σ,∆) =
|Z|α|I +BU |−αetr(−(I +B)U(I +BU)−1∆), (21a)
|I +B|αetr(−∆)|Z|α|I +BZ|−αetr(Z(I +BZ)−1(I +B)∆), (21b)
|I +B|αetr(− 1
∆(I −B)) (21c)
·|Z|α|I +BΩ|−αetr(1
Ω(I +BΩ)−1(I +B)∆(I −B)),
with the corresponding matrices B from (20) and Z,U,Ω from (15).
For the former series expansions the following relations were used:
Laplace transform f̂(t): f(x): F (x) =
f(ξ)dξ:
zαun vg
α+n(vx) G
α+n(vx) (22a)
zα+n vgα+n(vx) Gα+n(vx) (22b)
zαωn vhα,n(vx) Hα,n(vx) (22c)
where z = (1 + v−1t)−1, gα+n(x) = e
−xxα−1+n/Γ(α+ n),
α+n(x) =
gα+n(x) =
α− 1 + n
L(α−1)n (x)gα(x)
with the generalized Laguerre polynomials L
(α−1)
n and
hα,n(x) = (−1)n
α− 1 + n
L(α−1)n (2x)gα(x).
The last identity is verified by L.t.
The following bounds are derived from (22.14.13) in Abramowitz and Stegun (1965):
∣∣∣g(n)α+n(x)
∣∣∣ ≤
ex/2gα(x), α ≥ 1
2nα−1ex/2gα(x), 0 < α < 1
|hα,n(x)| ≤
xα−1/Γ(α), α ≥ 1
2nxα−1/Γ(α+ 1), 0 < α < 1
, (24)
matching with the conditions in (7).
The following generators (generating functions) with the Γ(α+n)–cdf Gα+n(x) are
required for the formulas in theorem 2:
Fα(x, y) =
n=0 G
α+n(x)y
n = 1
, |y| < 1, (25a)
n=0 Gα+n(x)y
n = Gα(x, y), y ∈ C, (25b)
n=0 Hα,n(x)y
n = 1
x, 2y
, |y| < 1 (25c)
The identities (a) and (c) are verified by the L.t. of fα(x, y) =
Fα(x, y). A short
calculation shows
Gα(x, y) =
Gα(x) − y1−αe(y−1)x Gα(xy)
, y 6= 1, α > 0
Gα−1(x) − y1−αe(y−1)x Gα−1(xy)
, α ≥ 1, G0 := 1
xgα(x) + (1 + x− α) Gα(x), y = 1
gα(x, y) =
Gα(x, y) =
gα(x) + y
1−αe(y−1)x Gα(xy), α > 0
y1−αe(y−1)x Gα−1(xy), α ≥ 1
The functions Fα(x, y) are especially simple for α ∈ N since Gα(z) = 1−e−z
j=0 z
j/j!,
α ∈ N.
Besides,
Gk+1/2(z) = erf(z
1/2)− e−z
zj−1/2
Γ(j + 1/2)
, k ∈ N0.
The following simple lemma is used for the proof of theorem 2.
Lemma 1. If B is a symmetrical p× p–matrix with ‖B‖ < 1 and
Y = diag(y1, . . . , yp) then the power series expansion
|I +BY |−α =
b(n1, . . . , np)
is absolutely convergent for max |yj | < rB = ‖B‖−1.
This follows from
(n) |b(n1, . . . , np)| = O(ϑn) with any ϑ > ‖B‖, which has been
already shown in (2.1.16) . . . (2.1.18) in Royen (1991) (with the notation −C instead of
4. The integral representations
In theorem 2 below the functions F (x1, . . . , xp;α1, . . . , αn,Σ1, . . . ,Σn,∆1, . . . ,∆n)
from (2) are represented by three different integrals over Cp = (−π, π]p. Together with
the generators Fα from (25), α =
k=1 αk, the following matrices are used with a scale
factor v to enforce ‖Bk‖ < 1:
Bk = vΣk − I, Dk = ∆k(I +Bk), Fα from (25a), (27a)
Bk = (vΣk)
−1 − I, Dk = (I +Bk)∆k, Fα from (25b), (27b)
Bk = 2(I + vΣk)
−1 − I, Dk = 12 (I +Bk)∆k(I −Bk), Fα from (25c). (27c)
Furthermore, we define λmax = max ‖Σk‖ , λ−1min = max ‖Σ
k ‖, yj = reiϕj ,
−π < ϕj ≤ π, Y = diag(y1, . . . , yp),
K = K(y1, . . . , yp) =
etr(±(Y +Bk)−1Dk)|I +BkY −1|−αk ,
where the negative sign occurs only with Bk, Dk from (27a), and
Fαdϕ =
j=1 Fα(vxj , yj)dϕj .
Theorem 2. With the above notations the functions F from (2) are respresentable
by each of the following three integrals:
(2π)−p
KFαdϕ, (28)
Fα from (25a), Bk, Dk from (27a), ‖Bk‖ < 1 if v < 2λ−1max, max ‖Bk‖ < r < 1,
etr(−∆k)|I +Bk|αk
(2π)−p
KFαdϕ, (29)
Fα from (25b), Bk, Dk from (27b), ‖Bk‖ < 1 if v > 12λ
min, max ‖Bk‖ < r,
∆k(I −Bk)
|I +Bk|αk
(2π)−p
KFαdϕ, (30)
Fα from (25c), Bk, Dk from (27c), v > 0, max ‖Bk‖ < r < 1.
Proof. Because of lemma 1 the assumptions of theorem 1 are satisfied with
ĝj0(t) = z
j = (1 + v
−1tj)
−α and hj(tj) corresponding to zj or uj = v
−1tjzj = 1 − zj
or ωj = zj − uj respectively. The functions ĝj0(t)(hj(t))n are the L.t. of the functions
in the second column of (22) from which type (a) and (c) have the bounds in (23), (24),
satisfying the condition (7) for theorem 1. The series
n=0 Gα+n(x)y
n = Gα(x, y) in
(25b) is absolutely convergent for every y ∈ C. Thus, all r > max ‖Bk‖ are admissible
in (29). In (30) we have max ‖Bk‖ < 1 for every v > 0. Hence, theorem 1 together with
the respresentations of the L.t. in (21) implies (28), (29) and (30).
The univariate case of (29) provides
F (x;α1, . . . , αn, σ
1 , . . . , σ
1 , . . . , δ
n) = (31)(
σ−2αkk e
Gα(vx, e
δ2k/(1 + vσ
iϕ − 1))
1 + (v−1σ−2k − 1)e−iϕ
with 2v > maxσ−2k , r = 1, Gα from (26). With p = 1 similar formulas arise from (28)
or (30).
The cdf of a quadratic form 1
x′Ax with T ′AT = diag(λ1, . . . , λn) ≥ 0 of rank q
and a N (µ, σ2In)–random vector x is a special case of (31) with αk = 12 , σ
k = λkσ
2 and
non–centrality parameters δ2k =
µ∗2k /σ
2, k = 1, . . . , q, µ∗ = T ′µ.
Some further remarks: In (29) also ‖Bk‖ > 1 is allowed since every r = ‖Y ‖ >
max ‖Bk‖ is admissible, which entails max ‖BkY −1‖ < 1.
With ϑ = λmax/λmin it follows with special values of v:
max ‖Bk‖ ≤
in (28) with v = 2(λmin + λmax)
max ‖Bk‖ ≤
in (29) with v =
(λ−1min + λ
max),
max ‖Bk‖ ≤
ϑ− 1√
in (30) with v = (λminλmax)
−1/2.
More generally, the scale factor v = w2 can be replaced by a scale matrix
W 2 = diag(w21 , . . . , w
p) > 0. Then with Tw = W
−1TW−1, Σw = WΣW , ∆w = W∆W
the L.t. (1) equals
|I +ΣwTw|−α etr(−ΣwTw(I +ΣwTw)−1∆w). (32)
Consequently, besides the substitutions vΣk → WΣkW , ∆k → W∆kW−1, the matrices
I + Bk in theorem 2 must be replaced by WΣkW , (WΣkW )
−1 and 2(I + WΣkW )
respectively, and the generators Fα(vxj , yj) by Fα(w
jxj , yj).
In particular for a single Γp(α,Σ,∆)–distribution this more general scaling can be
used to minimize ‖B‖ or for a ”natural scaling” i.e. to standardize I+B to a correlation
matrix. However, ‖B‖ < 1 must be taken into account in (28), whereas this condition
is satisfied in (30) for every scaling. It was shown in Royen (1991) that natural scaling
can always be accomplished also in I +B = 2(I +WΣW )−1 by a unique W 2.
5. Representations of the Γp(α,Σ,∆) distribution function by (p− 1)–variate
integrals
For a single Γp(α,Σ,∆)–cdf it is always possible to perform the integration over a
single variable ϕj within the integrals from theorem 2.
We use the following functions
Gα(x, y) = e−y
Gα+n(x)
α+n(x)
(−y)n
x, y ∈ C, Gα+n, G(n)α+n from (22), and
G∗α(x, y) = eyGα(x, y).
For positive half integers α = 1/2 + k these functions can also be computed by the
erf–function and a sum of k terms which are essentially given by the modified Bessel
functions Ij−1/2(2(xy)
1/2), j = 1, . . . , k, (see e.g. Royen (1995) or (2006)).
Now let be W 2 = diag(w21 , . . . , w
p) a general scale matrix,
Y = diag(y1, . . . , yp), yj = re
iϕj , −π < ϕj ≤ π,
Bpp bp
b′p bpp
WΣW − I, (34a)
(WΣW )−1 − I, (34b)
2(I +WΣW )−1 − I, (34c)
Dpp dp
dp dpp
W∆ΣW, (35a)
W−1Σ−1∆W−1, (35b)
2(I +WΣW )−1W∆W−1(I − (I +WΣW )−1), (35c)
y0 = y0(y1, . . . , yp−1) = b
p(Ypp +Bpp)
−1bp − bpp (36)
q = q(y1, . . . , yp−1) = (b
p(Ypp +Bpp)
−1,−1)D
(Ypp +Bpp)
Kα = Kα(y1, . . . , yp−1) = etr(±(Ypp +Bpp)−1Dpp)|I +BppY −1pp |−α,
where the negative sign is only taken for Bpp from (34a).
Theorem 3. With the above notations the Γp(α,Σ,∆)–cdf F (x1, . . . , xp;α,Σ,∆)
is given by each of the following three integrals:
(2π)p−1
w2pxp
1− y0
1− y0
1− yj
w2jxj ,
yj − 1
dϕj , (38)
B from (34a), D from (35a), ‖B‖ < r < 1,
etr(−W∆W−1)
|WΣW |α ·
(2π)p−1
(1− y0)−α G∗α
(1− y0)w2pxp,
1− y0
jxj , yj)dϕj ,
B from (34b), D from (35b), ‖B‖ < r,
2αpetr(− 1
W∆W−1(I −B))
|I +WΣW |α ·
(2π)p−1
1− y0
(1− y0)−αGα
1− y0
1 + y0
w2pxp,
1− y20
1 + yj
w2jxj ,
yj + 1
dϕj ,
B from (34c), D from (35c), ‖B‖ < r < 1.
For the proof of theorem 3 the following two lemmas are required.
Lemma 2. With Y = diag(y1, . . . , yp), yj = re
iϕj , ‖B‖ < r, B,D, y0, q from (34),
(35), (36), (37) the following decomposition is obtained
etr((Y +B)−1D)|Y +B|−α
= etr
(Ypp +Bpp)
−1Dpp
|Ypp +Bpp|−α exp
yp − y0
(yp − y0)−α.
Proof. From frequently used formulas for p× p–matrices, (see e.g. complements
and problems 2.4, 2.7 in chapter 1b of Rao (1973)) it follows for
A = Y +B =
App bp
b′p yp + bpp
|A| = |App|(yp + bpp − b′pA−1pp bp) = |Ypp +Bpp|(yp − y0),
A−1 =
A−1pp +
yp−y0
A−1pp bpb
pp − 1yp−y0A
pp bp
yp−y0
yp−y0
trace(A−1D)
= trace
A−1pp Dpp +
yp − y0
A−1pp bpb
pp Dpp −A−1pp bpdp
yp − y0
(dpp − b′pA−1pp dp)
= trace(A−1pp Dpp) +
yp − y0
, which implies (41).
Lemma 3. Let be q any number, Sr = {y ∈ C
∣∣|y| = r}, y0 any number with
|y0| < r, then with Fα from (25), Gα,G∗α from (33), and the negative sign in ±q only for
(42a)
y − y0
Fα(x, y)(y − y0)−αyα−1dy
Fα from (25a), r < 1, (42a)
(1− y0)−α G∗α
(1 − y0)x, q1−y0
, Fα from (25b), (42b)
(1− y0)−α Gα
x, 2q
, Fα from (25c), r < 1, (42c)
Proof. It is sufficient to verify (42) for the corresponding derivatives fα =
At first, (42a) is shown:
With Fα from (25a) and the binomial series for (1− y0/y)−(α+n) we obtain
fα(x, y)(y − y0)−(α+n)yα−1dy
α+m(x)y
α+ n+ k − 1
y−n−1dy.
With z = (1 + t)−1, u = tz, the last integral has the L.t.
(uy)m
α+ n+ k − 1
y−n−1dy
m=n+k
α+ n+ k − 1
yk0 = z
αun(1− uy0)−(α+n).
Multiplication by (−q)n/n! and summation over n leads to the L.t.
(1− uy0)α
1− uy0
(1 + (1− y0)t)α
(1− y0)t
1 + (1− y0)t
and this is the L.t. of ∂
To verify (42b) we obtain with Fα from (25b):
fα(x, y)(y − y0)−(α+n)yα−1dy
gα+m(x)y
Γ(α+ n+ k)
Γ(α+ n)k!
y−n−1dy
m=n+k
gα+m(x)
Γ(α + n+ k)
Γ(α+ n)k!
yk0 = gα+n(x)e
= (1− y0)−(α+n)(1− y0)gα+n((1 − y0)x).
Multiplication by qn/n! and summation provides (1− y0)−α ∂∂x G
(1− y0)x, q1−y0
(42c) can be shown by L.t. in a similar way as (42a).
Proof of theorem 3. Without loss of generality yp is selected from the variables
yj = re
iϕj in Y = diag(y1, . . . , yp) with any fixed r > ‖B‖. If yp is replaced by a variable
y with any |y| then the equation
|Y +B| = |Ypp +Bpp|(yp + bpp − b′p(Ypp +Bpp)−1bp) = 0
has always a unique solution
y = y0 = b
p(Ypp +Bpp)
−1bp − bpp
with |y0| < r since ‖Bpp‖ ≤ ‖B‖.
Hence, with lemma 2 and lemma 3, theorem 3 is obtained by integration over ϕp in
the integrals of theorem 2 with n = 1.
References
Abramowitz, M. and Stegun, I.A. (1968). Handbook of Mathematical Functions, Dover, New York.
Bapat, R.B. (1989). Infinite divisibility of multivariate gamma distributions and M–matrices, Sankhyā,
Series A 51, 73–78.
Blacher, R. (2003). Multivariate quadratic forms of random vectors, Journal of Multivariate Analysis
87, 2–23.
Fang, B.Q. (2005). Noncentral quadratic forms of the skew elliptical variables, Journal of Multivariate
Analysis 95, 410–430.
Griffiths, R.C. (1984). Characterization of infinitely divisible multivariate gamma distributions, Journal
of Multivariate Analysis 15, 13–20.
Khatri, C.G., Krishnaiah, P.R. and Sen, P.K. (1977). A note on the joint distribution of correlated
quadratic forms, Journal of Statistical Planning and Inference 1, 299–307.
Krishnamoorthy, A.S. and Parthasarathy, M. (1951). A multivariate gamma type distribution, Annals
of Mathematical Statistics 22, 549–557 (correction: ibid. (1960), 31, p. 229).
Mathai, A.M. and Provost, S.B. (1992). Quadratic forms in random variables: Theory and applications,
Marcel Dekker, New York.
Rao, C.R. (1973). Linear Statistical Inference and its Applications, 2nd edition, Wiley, New York.
Royen, T. (1991). Expansions for the multivariate chi–square distribution, Journal of Multivariate
Analysis 38, 213–232.
Royen, T. (1992). On representation and computation of multivariate gamma distributions, in: Data
Analysis and Statistical Inference - Festschrift in Honour of Friedhelm Eicker, 201–216, Verlag
Josef Eul, Bergisch Gladbach, Köln.
Royen, T. (1995). On some central and non–central multivariate chi–square distributions, Statistica
Sinica 5, 373–397.
Royen, T. (1997). Multivariate gamma distributions (Update), Encyclopedia of Statistical Sciences,
Update Volume 1, 419–425, Wiley, New York.
Royen, T. (2006). Integral representations and approximations for multivariate gamma distributions,
Annals of the Institute of Statistical Mathematics, DOI 10.1007/s10463-006-0057-5.
	Introduction 
	The method 
	Three representations for the p(,,)–Laplace transform and the related generators
	The integral representations
	Representations of the p(,,) distribution function by (p-1)–variate integrals

  Three types of integral representations for the cumulative distribution
functions of convolutions of non-central p-variate gamma distributions are
given by integration of elementary complex functions over the p-cube Cp =
(-pi,pi]x...x(-pi,pi]. In particular, the joint distribution of the diagonal
elements of a generalized quadratic form XAX' with n independent normally
distributed column vectors in X is obtained. For a single p-variate gamma
distribution function (p-1)-variate integrals over Cp-1 are derived. The
integrals are numerically more favourable than integrals obtained from the
Fourier or laplace inversion formula.

Introduction
The following notations are used:
(n) stands for
n1+...+np=n
n1, . . . , np ∈ N0 and
without any indices means
(n). The notation D ≥ 0
is also used for non–symmetrical matrices Dp×p with only non–negative eigenvalues.
The spectral norm of a p× p–matrix B is denoted by ‖B‖, I or Ip is always an identity
matrix and Cp is the p–cube (−π, π]p.
The Laplace transform (L.t.) of a p–variate non–central Γp(α,Σ,∆)–density with
α > 0, Σ > 0 and a non–centrality matrix ∆ ≥ 0 was originally obtained from the L.t.
of a non-central Wp(2α,Σ,∆)–Wishart distribution (with an additional scale factor 2)
and is given by
f̂(t1, . . . , tp;α,Σ,∆) = |Ip +ΣT |−αetr(−ΣT (I +ΣT )−1∆), (1)
T = diag(t1, . . . , tp), t1, . . . , tp ≥ 0.
This function f̂ is generally the L.t. of the density of a real measure on (0,∞)p which
is not always a probability measure. The term ”Γp(α,Σ,∆)–distribution” is used here
in this general sense. The exact set of values α, leading to a probability density (pdf)
f(x1, . . . , xp;α,Σ,∆), depends on Σ and presumedly on ∆. To obtain a pdf, all positive
integers 2α (degrees of freedom) are admissible and all 2α > p − 1. Moreover, in the
central case all non–integer values 2α > p− 2 ≥ 0 are allowed. For p− 2 < 2α < p− 1
http://arxiv.org/abs/0704.0539v1
see Royen (1997). Furthermore all α > 0 are admissible if |I + ΣT |−1 is infinitely
divisible. Two characterizations of infinite divisibility of a Γp(α,Σ)–distribution are
found in Griffiths (1984) and Bapat (1989). Further conditions for admissible non–
integer 2α < p− 2 are given in Royen (1997), (2006).
Three integral representations by integration over Cp are provided by theorem 2 in
section 4 for the functions
F (x1, . . . , xp;α1, . . . , αn,Σ1, . . . ,Σn,∆1, . . . ,∆n) (2)
. . .
f(ξ1, . . . , ξp;α1, . . . , αn,Σ1, . . . ,Σn,∆1, . . . ,∆n)dξ1 . . . dξp,
where f has the L.t.
|Ip +ΣkT |−αketr(−ΣkT (I +ΣkT )−1∆k), (3)
α1, . . . , αn > 0,Σ1, . . . ,Σn > 0,∆1, . . . ,∆n ≥ 0.
Thus, F is not always the cumulative distribution function (cdf) of a probability measure.
In particular letXp×n be aNp×n(Mp×n,Σp×p⊗In)–random matrix and An×n ≥ 0 of
rank q with T ′AT = Λ = diag(λ1, . . . , λn), λ1 ≥ . . . ≥ λn. Then the joint distribution of
the diagonal elements of the generalized quadratic form 1
XAX ′ equals the distribution
of the diagonal of 1
Y ΛY ′ with a Np×n(MT,Σp×p⊗In)–distributed Y = XT . This is the
distribution of a sum of q independent Γp(
, λkΣ,∆k =
−1)–random vectors,
where µ∗k is the k–th column of M
∗ = MT . This joint distribution of p quadratic
forms of normal random vectors is comprised within theorem 2 as a special case with
, Σk = λkΣ, k = 1, . . . , q. For methods under more general assumptions see also
Blacher (2003). For a survey of univariate quadratic forms of normal random variables
see chapter 4 in Mathai and Provost (1992). For several quadratic forms of skew elliptical
distributions see B.Q. Fang (2005).
In Royen (1991), (1992) three different types of series expansions for the χ2p(2α,Σ)–
cdf were derived from three different representations of the χ2p(2α,Σ)–L.t. which are
extended to the general Γp(α,Σ,∆)–L.t. in section 3 in a similar way as in Royen
(1995).
Some series expansions, closely related to the first two types, are already found in
Khatri, Krishnaiah and Sen (1977). The third type was introduced because of its superior
convergence properties. The simple method to transform many series expansions into
integrals over Cp is explained in more detail in section 2 and summarized in theorem 1.
The idea is as follows:
If A(z1, . . . , zp) and B(z1, . . . , zp) are analytical functions whose power series have
the coefficients a(m1, . . . ,mp) and b(n1, . . . , np) and which are absolutely convergent for
max |zj | < rA and max |zj | < rB respectively, where r−1B < rA, then
(2π)−p
A(y1, . . . , yp)B(y
1 , . . . , y
p )dϕ1 . . . dϕp (4)
a(n1, . . . , np)b(n1, . . . , np)
holds with yj = re
iϕj , −π < ϕj ≤ π, j = 1, . . . , p and r−1B < r < rA.
The integrals in (4) might be more economical than the series if the generating
functions A and B are simple available functions and if the series are slowly convergent
with very intricate coefficients. For non–central multivariate gamma distributions series
expansions are practically not feasible.
The integral representations in theorem 2 of section 4 are of the type in (4). As
long as no elementary density formulas are availale it should be a reasonable way to
obtain the joint cdf by integration of elementary terms only over Cp and not over Rp as
by the Fourier or Laplace inversion formula. A single Γp(α,Σ,∆)–cdf is represented by
a (p− 1)–variate integral over Cp−1 in section 5.
A totally different
–variate integral representation of the Γp(α,Σ)–cdf has
been given recently by Royen (2006), which is based on m–factorial decompositions∑
p×p = D −BB′, where D is a real or complex diagonal matrix minimizing the rank
m of Σ−1 −D. Approximations to a Γp(α,Σ)–cdf are obtained by m–factorial approx-
imations to Σ with a low value of m. These approximations are improved further by
successive correction terms.
2. The method
Theorem 1 in this section can be generalized in many ways, e.g. for Fourier trans-
forms, but the version below is sufficient for the purpose of the underlying paper.
Let f̂(t1, . . . , tp), t1, . . . , tp ≥ 0, be a given L.t. of an unknown function f(x1, . . . , xp)
with f = 0 for minxj < 0. It is assumed that there are univariate L.t. ĝj0(t) of
some probability densities gj0(x) on (0,∞) and further functions hj(t) with |hj(t)| ≤ 1,
uniformly for t ≥ 0, which enable a representation
f̂(t1, . . . , tp) =
ĝj0(tj)
B (h1(t1), . . . , hp(tp)) (5)
with an analytical function B(z1, . . . , zp) whose power series expansion
b(n1, . . . , np)
j (6)
is absolutely convergent for |z1|, . . . , |zp| < rB with a certain value rB > 1.
Furthermore, the products ĝj0(t)(hj(t))
n are supposed to be the L.t. of continuous
functions gjn(x), x > 0, which satisfy the conditions
|gjn(x)| ≤ nck(x) with a constant c and
(7)∫ ∞
k(x)e−txdx < ∞ for all t > 0 .
Hence, the generating functions (generators)
gj(x, y) =
gjn(x)y
n, j = 1, . . . , p, (8)
are defined for all x > 0 and |y| < 1, and they have the L.t.
ĝj(t, y) =
ĝj0(t)
1− yhj(t)
, t ≥ 0. (9)
Theorem 1. Under the assumptions from (6) and (7) f̂ in (5) is the L.t. of
f(x1, . . . , xp) = (2π)
B(y−11 , . . . , y
gj(xj , yj)dϕj (10)
with yj = re
iϕj , −π < ϕ ≤ π, r−1B < r < 1, gj from (8).
Proof. The integral in (10) is evaluated by
(2π)−p
b(m1, . . . ,mp)
gjnj (xj)y
 dϕ1 . . . dϕp
b(n1, . . . , np)
j=1 gjnj (xj)
and this series has the L.t. from (5).
Some further remarks: With
Gj(xj , yj) =
gj(ξ, yj)dξ (11)
instead of the gj in (10), the corresponding representation arises for
F (x1, . . . , xp) =
. . .
f(ξ1, . . . , ξp)dξ1 . . . dξp. (12)
If the series in (8) are absolutely convergent for all y ∈ C then additionally
hj(t) = 0 (13)
is supposed to hold. Then the rhs of (9) is the L.t. of gj(x, y) for any fixed y and all
sufficiently large t.
In some cases the functions gj0 and their L.t. ĝj0 are known from univariate marginal
distributions apart from some scale factors. If the functions uj = hj(t) are explicitly
invertible then
B(u1, . . . , up) =
h−11 (u1), . . . , h
p (up)
j=1 ĝj0
h−1j (uj)
) (14)
can sometimes be found easily from the given f̂ .
3. Three representations for the Γp(α,Σ,∆)–Laplace transform
and the related generators
With any v > 0 we define
zj = (1 + v
−1tj)
−1, tj ≥ 0, uj = 1− zj = v−1tjzj , ωj = zj − uj ,
Z = diag(z1, . . . , zp), U = diag(u1, . . . , up), Ω = diag(ω1, . . . , ωp).
The scale factor v is introduced to obtain ‖B‖ < 1 for the matrices B defined in (20)
below and to effect the convergence of some series expansions. For a more general scaling
see remarks following theorem 2 in section 4.
From the relations
v−1T = UZ−1, Ip = Z + U, Ω = Z − U, (16)
it follows for the matrices I +ΣT in the L.t. (1):
I +ΣT = I + vΣUZ−1 = (Z + vΣU)Z−1 (17)
Z + vΣU =
I + (vΣ− I)U, (18a)
vΣ(I + (v−1Σ−1 − I)Z), (18b)
(I + vΣ)(I + (2(I + vΣ)−1 − I)Ω), (18c)
and therefore
|I +ΣT |−α = cα|Z|α|I +BY |−α (19)
Y = U, B = vΣ− I, c = 1, (20a)
Y = Z, B = (vΣ)−1 − I, c = |I +B|, (20b)
Y = Ω, B = 2(I + vΣ)−1 − I, c = |I +B|. (20c)
It should be noticed that ‖B‖ < 1 in (20c) for every v > 0 and Σ > 0.
Now, using (16), by a straightforward calculation the L.t. in (1) can be represented
f̂(t1, . . . , tp;α,Σ,∆) =
|Z|α|I +BU |−αetr(−(I +B)U(I +BU)−1∆), (21a)
|I +B|αetr(−∆)|Z|α|I +BZ|−αetr(Z(I +BZ)−1(I +B)∆), (21b)
|I +B|αetr(− 1
∆(I −B)) (21c)
·|Z|α|I +BΩ|−αetr(1
Ω(I +BΩ)−1(I +B)∆(I −B)),
with the corresponding matrices B from (20) and Z,U,Ω from (15).
For the former series expansions the following relations were used:
Laplace transform f̂(t): f(x): F (x) =
f(ξ)dξ:
zαun vg
α+n(vx) G
α+n(vx) (22a)
zα+n vgα+n(vx) Gα+n(vx) (22b)
zαωn vhα,n(vx) Hα,n(vx) (22c)
where z = (1 + v−1t)−1, gα+n(x) = e
−xxα−1+n/Γ(α+ n),
α+n(x) =
gα+n(x) =
α− 1 + n
L(α−1)n (x)gα(x)
with the generalized Laguerre polynomials L
(α−1)
n and
hα,n(x) = (−1)n
α− 1 + n
L(α−1)n (2x)gα(x).
The last identity is verified by L.t.
The following bounds are derived from (22.14.13) in Abramowitz and Stegun (1965):
∣∣∣g(n)α+n(x)
∣∣∣ ≤
ex/2gα(x), α ≥ 1
2nα−1ex/2gα(x), 0 < α < 1
|hα,n(x)| ≤
xα−1/Γ(α), α ≥ 1
2nxα−1/Γ(α+ 1), 0 < α < 1
, (24)
matching with the conditions in (7).
The following generators (generating functions) with the Γ(α+n)–cdf Gα+n(x) are
required for the formulas in theorem 2:
Fα(x, y) =
n=0 G
α+n(x)y
n = 1
, |y| < 1, (25a)
n=0 Gα+n(x)y
n = Gα(x, y), y ∈ C, (25b)
n=0 Hα,n(x)y
n = 1
x, 2y
, |y| < 1 (25c)
The identities (a) and (c) are verified by the L.t. of fα(x, y) =
Fα(x, y). A short
calculation shows
Gα(x, y) =
Gα(x) − y1−αe(y−1)x Gα(xy)
, y 6= 1, α > 0
Gα−1(x) − y1−αe(y−1)x Gα−1(xy)
, α ≥ 1, G0 := 1
xgα(x) + (1 + x− α) Gα(x), y = 1
gα(x, y) =
Gα(x, y) =
gα(x) + y
1−αe(y−1)x Gα(xy), α > 0
y1−αe(y−1)x Gα−1(xy), α ≥ 1
The functions Fα(x, y) are especially simple for α ∈ N since Gα(z) = 1−e−z
j=0 z
j/j!,
α ∈ N.
Besides,
Gk+1/2(z) = erf(z
1/2)− e−z
zj−1/2
Γ(j + 1/2)
, k ∈ N0.
The following simple lemma is used for the proof of theorem 2.
Lemma 1. If B is a symmetrical p× p–matrix with ‖B‖ < 1 and
Y = diag(y1, . . . , yp) then the power series expansion
|I +BY |−α =
b(n1, . . . , np)
is absolutely convergent for max |yj | < rB = ‖B‖−1.
This follows from
(n) |b(n1, . . . , np)| = O(ϑn) with any ϑ > ‖B‖, which has been
already shown in (2.1.16) . . . (2.1.18) in Royen (1991) (with the notation −C instead of
4. The integral representations
In theorem 2 below the functions F (x1, . . . , xp;α1, . . . , αn,Σ1, . . . ,Σn,∆1, . . . ,∆n)
from (2) are represented by three different integrals over Cp = (−π, π]p. Together with
the generators Fα from (25), α =
k=1 αk, the following matrices are used with a scale
factor v to enforce ‖Bk‖ < 1:
Bk = vΣk − I, Dk = ∆k(I +Bk), Fα from (25a), (27a)
Bk = (vΣk)
−1 − I, Dk = (I +Bk)∆k, Fα from (25b), (27b)
Bk = 2(I + vΣk)
−1 − I, Dk = 12 (I +Bk)∆k(I −Bk), Fα from (25c). (27c)
Furthermore, we define λmax = max ‖Σk‖ , λ−1min = max ‖Σ
k ‖, yj = reiϕj ,
−π < ϕj ≤ π, Y = diag(y1, . . . , yp),
K = K(y1, . . . , yp) =
etr(±(Y +Bk)−1Dk)|I +BkY −1|−αk ,
where the negative sign occurs only with Bk, Dk from (27a), and
Fαdϕ =
j=1 Fα(vxj , yj)dϕj .
Theorem 2. With the above notations the functions F from (2) are respresentable
by each of the following three integrals:
(2π)−p
KFαdϕ, (28)
Fα from (25a), Bk, Dk from (27a), ‖Bk‖ < 1 if v < 2λ−1max, max ‖Bk‖ < r < 1,
etr(−∆k)|I +Bk|αk
(2π)−p
KFαdϕ, (29)
Fα from (25b), Bk, Dk from (27b), ‖Bk‖ < 1 if v > 12λ
min, max ‖Bk‖ < r,
∆k(I −Bk)
|I +Bk|αk
(2π)−p
KFαdϕ, (30)
Fα from (25c), Bk, Dk from (27c), v > 0, max ‖Bk‖ < r < 1.
Proof. Because of lemma 1 the assumptions of theorem 1 are satisfied with
ĝj0(t) = z
j = (1 + v
−1tj)
−α and hj(tj) corresponding to zj or uj = v
−1tjzj = 1 − zj
or ωj = zj − uj respectively. The functions ĝj0(t)(hj(t))n are the L.t. of the functions
in the second column of (22) from which type (a) and (c) have the bounds in (23), (24),
satisfying the condition (7) for theorem 1. The series
n=0 Gα+n(x)y
n = Gα(x, y) in
(25b) is absolutely convergent for every y ∈ C. Thus, all r > max ‖Bk‖ are admissible
in (29). In (30) we have max ‖Bk‖ < 1 for every v > 0. Hence, theorem 1 together with
the respresentations of the L.t. in (21) implies (28), (29) and (30).
The univariate case of (29) provides
F (x;α1, . . . , αn, σ
1 , . . . , σ
1 , . . . , δ
n) = (31)(
σ−2αkk e
Gα(vx, e
δ2k/(1 + vσ
iϕ − 1))
1 + (v−1σ−2k − 1)e−iϕ
with 2v > maxσ−2k , r = 1, Gα from (26). With p = 1 similar formulas arise from (28)
or (30).
The cdf of a quadratic form 1
x′Ax with T ′AT = diag(λ1, . . . , λn) ≥ 0 of rank q
and a N (µ, σ2In)–random vector x is a special case of (31) with αk = 12 , σ
k = λkσ
2 and
non–centrality parameters δ2k =
µ∗2k /σ
2, k = 1, . . . , q, µ∗ = T ′µ.
Some further remarks: In (29) also ‖Bk‖ > 1 is allowed since every r = ‖Y ‖ >
max ‖Bk‖ is admissible, which entails max ‖BkY −1‖ < 1.
With ϑ = λmax/λmin it follows with special values of v:
max ‖Bk‖ ≤
in (28) with v = 2(λmin + λmax)
max ‖Bk‖ ≤
in (29) with v =
(λ−1min + λ
max),
max ‖Bk‖ ≤
ϑ− 1√
in (30) with v = (λminλmax)
−1/2.
More generally, the scale factor v = w2 can be replaced by a scale matrix
W 2 = diag(w21 , . . . , w
p) > 0. Then with Tw = W
−1TW−1, Σw = WΣW , ∆w = W∆W
the L.t. (1) equals
|I +ΣwTw|−α etr(−ΣwTw(I +ΣwTw)−1∆w). (32)
Consequently, besides the substitutions vΣk → WΣkW , ∆k → W∆kW−1, the matrices
I + Bk in theorem 2 must be replaced by WΣkW , (WΣkW )
−1 and 2(I + WΣkW )
respectively, and the generators Fα(vxj , yj) by Fα(w
jxj , yj).
In particular for a single Γp(α,Σ,∆)–distribution this more general scaling can be
used to minimize ‖B‖ or for a ”natural scaling” i.e. to standardize I+B to a correlation
matrix. However, ‖B‖ < 1 must be taken into account in (28), whereas this condition
is satisfied in (30) for every scaling. It was shown in Royen (1991) that natural scaling
can always be accomplished also in I +B = 2(I +WΣW )−1 by a unique W 2.
5. Representations of the Γp(α,Σ,∆) distribution function by (p− 1)–variate
integrals
For a single Γp(α,Σ,∆)–cdf it is always possible to perform the integration over a
single variable ϕj within the integrals from theorem 2.
We use the following functions
Gα(x, y) = e−y
Gα+n(x)
α+n(x)
(−y)n
x, y ∈ C, Gα+n, G(n)α+n from (22), and
G∗α(x, y) = eyGα(x, y).
For positive half integers α = 1/2 + k these functions can also be computed by the
erf–function and a sum of k terms which are essentially given by the modified Bessel
functions Ij−1/2(2(xy)
1/2), j = 1, . . . , k, (see e.g. Royen (1995) or (2006)).
Now let be W 2 = diag(w21 , . . . , w
p) a general scale matrix,
Y = diag(y1, . . . , yp), yj = re
iϕj , −π < ϕj ≤ π,
Bpp bp
b′p bpp
WΣW − I, (34a)
(WΣW )−1 − I, (34b)
2(I +WΣW )−1 − I, (34c)
Dpp dp
dp dpp
W∆ΣW, (35a)
W−1Σ−1∆W−1, (35b)
2(I +WΣW )−1W∆W−1(I − (I +WΣW )−1), (35c)
y0 = y0(y1, . . . , yp−1) = b
p(Ypp +Bpp)
−1bp − bpp (36)
q = q(y1, . . . , yp−1) = (b
p(Ypp +Bpp)
−1,−1)D
(Ypp +Bpp)
Kα = Kα(y1, . . . , yp−1) = etr(±(Ypp +Bpp)−1Dpp)|I +BppY −1pp |−α,
where the negative sign is only taken for Bpp from (34a).
Theorem 3. With the above notations the Γp(α,Σ,∆)–cdf F (x1, . . . , xp;α,Σ,∆)
is given by each of the following three integrals:
(2π)p−1
w2pxp
1− y0
1− y0
1− yj
w2jxj ,
yj − 1
dϕj , (38)
B from (34a), D from (35a), ‖B‖ < r < 1,
etr(−W∆W−1)
|WΣW |α ·
(2π)p−1
(1− y0)−α G∗α
(1− y0)w2pxp,
1− y0
jxj , yj)dϕj ,
B from (34b), D from (35b), ‖B‖ < r,
2αpetr(− 1
W∆W−1(I −B))
|I +WΣW |α ·
(2π)p−1
1− y0
(1− y0)−αGα
1− y0
1 + y0
w2pxp,
1− y20
1 + yj
w2jxj ,
yj + 1
dϕj ,
B from (34c), D from (35c), ‖B‖ < r < 1.
For the proof of theorem 3 the following two lemmas are required.
Lemma 2. With Y = diag(y1, . . . , yp), yj = re
iϕj , ‖B‖ < r, B,D, y0, q from (34),
(35), (36), (37) the following decomposition is obtained
etr((Y +B)−1D)|Y +B|−α
= etr
(Ypp +Bpp)
−1Dpp
|Ypp +Bpp|−α exp
yp − y0
(yp − y0)−α.
Proof. From frequently used formulas for p× p–matrices, (see e.g. complements
and problems 2.4, 2.7 in chapter 1b of Rao (1973)) it follows for
A = Y +B =
App bp
b′p yp + bpp
|A| = |App|(yp + bpp − b′pA−1pp bp) = |Ypp +Bpp|(yp − y0),
A−1 =
A−1pp +
yp−y0
A−1pp bpb
pp − 1yp−y0A
pp bp
yp−y0
yp−y0
trace(A−1D)
= trace
A−1pp Dpp +
yp − y0
A−1pp bpb
pp Dpp −A−1pp bpdp
yp − y0
(dpp − b′pA−1pp dp)
= trace(A−1pp Dpp) +
yp − y0
, which implies (41).
Lemma 3. Let be q any number, Sr = {y ∈ C
∣∣|y| = r}, y0 any number with
|y0| < r, then with Fα from (25), Gα,G∗α from (33), and the negative sign in ±q only for
(42a)
y − y0
Fα(x, y)(y − y0)−αyα−1dy
Fα from (25a), r < 1, (42a)
(1− y0)−α G∗α
(1 − y0)x, q1−y0
, Fα from (25b), (42b)
(1− y0)−α Gα
x, 2q
, Fα from (25c), r < 1, (42c)
Proof. It is sufficient to verify (42) for the corresponding derivatives fα =
At first, (42a) is shown:
With Fα from (25a) and the binomial series for (1− y0/y)−(α+n) we obtain
fα(x, y)(y − y0)−(α+n)yα−1dy
α+m(x)y
α+ n+ k − 1
y−n−1dy.
With z = (1 + t)−1, u = tz, the last integral has the L.t.
(uy)m
α+ n+ k − 1
y−n−1dy
m=n+k
α+ n+ k − 1
yk0 = z
αun(1− uy0)−(α+n).
Multiplication by (−q)n/n! and summation over n leads to the L.t.
(1− uy0)α
1− uy0
(1 + (1− y0)t)α
(1− y0)t
1 + (1− y0)t
and this is the L.t. of ∂
To verify (42b) we obtain with Fα from (25b):
fα(x, y)(y − y0)−(α+n)yα−1dy
gα+m(x)y
Γ(α+ n+ k)
Γ(α+ n)k!
y−n−1dy
m=n+k
gα+m(x)
Γ(α + n+ k)
Γ(α+ n)k!
yk0 = gα+n(x)e
= (1− y0)−(α+n)(1− y0)gα+n((1 − y0)x).
Multiplication by qn/n! and summation provides (1− y0)−α ∂∂x G
(1− y0)x, q1−y0
(42c) can be shown by L.t. in a similar way as (42a).
Proof of theorem 3. Without loss of generality yp is selected from the variables
yj = re
iϕj in Y = diag(y1, . . . , yp) with any fixed r > ‖B‖. If yp is replaced by a variable
y with any |y| then the equation
|Y +B| = |Ypp +Bpp|(yp + bpp − b′p(Ypp +Bpp)−1bp) = 0
has always a unique solution
y = y0 = b
p(Ypp +Bpp)
−1bp − bpp
with |y0| < r since ‖Bpp‖ ≤ ‖B‖.
Hence, with lemma 2 and lemma 3, theorem 3 is obtained by integration over ϕp in
the integrals of theorem 2 with n = 1.
References
Abramowitz, M. and Stegun, I.A. (1968). Handbook of Mathematical Functions, Dover, New York.
Bapat, R.B. (1989). Infinite divisibility of multivariate gamma distributions and M–matrices, Sankhyā,
Series A 51, 73–78.
Blacher, R. (2003). Multivariate quadratic forms of random vectors, Journal of Multivariate Analysis
87, 2–23.
Fang, B.Q. (2005). Noncentral quadratic forms of the skew elliptical variables, Journal of Multivariate
Analysis 95, 410–430.
Griffiths, R.C. (1984). Characterization of infinitely divisible multivariate gamma distributions, Journal
of Multivariate Analysis 15, 13–20.
Khatri, C.G., Krishnaiah, P.R. and Sen, P.K. (1977). A note on the joint distribution of correlated
quadratic forms, Journal of Statistical Planning and Inference 1, 299–307.
Krishnamoorthy, A.S. and Parthasarathy, M. (1951). A multivariate gamma type distribution, Annals
of Mathematical Statistics 22, 549–557 (correction: ibid. (1960), 31, p. 229).
Mathai, A.M. and Provost, S.B. (1992). Quadratic forms in random variables: Theory and applications,
Marcel Dekker, New York.
Rao, C.R. (1973). Linear Statistical Inference and its Applications, 2nd edition, Wiley, New York.
Royen, T. (1991). Expansions for the multivariate chi–square distribution, Journal of Multivariate
Analysis 38, 213–232.
Royen, T. (1992). On representation and computation of multivariate gamma distributions, in: Data
Analysis and Statistical Inference - Festschrift in Honour of Friedhelm Eicker, 201–216, Verlag
Josef Eul, Bergisch Gladbach, Köln.
Royen, T. (1995). On some central and non–central multivariate chi–square distributions, Statistica
Sinica 5, 373–397.
Royen, T. (1997). Multivariate gamma distributions (Update), Encyclopedia of Statistical Sciences,
Update Volume 1, 419–425, Wiley, New York.
Royen, T. (2006). Integral representations and approximations for multivariate gamma distributions,
Annals of the Institute of Statistical Mathematics, DOI 10.1007/s10463-006-0057-5.
	Introduction 
	The method 
	Three representations for the p(,,)–Laplace transform and the related generators
	The integral representations
	Representations of the p(,,) distribution function by (p-1)–variate integrals
