Integral representations for convolutions of non-central multivariate gamma distributions

TL;DR

This paper presents new integral representations for the cumulative distribution functions of convolutions of non-central multivariate gamma distributions, facilitating more efficient numerical computation and analysis.

Contribution

It introduces three types of integral formulas for these distributions, including joint distributions of quadratic forms, improving computational efficiency over traditional methods.

Findings

Derived integral representations for convolutions of non-central multivariate gamma distributions.

Obtained joint distribution of diagonal elements of quadratic forms involving normal vectors.

Integral formulas are numerically more favorable than Fourier or Laplace inversion methods.

Abstract

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.

Equations

$$\hat{f}(t_1,\dots ,t_p;\alpha ,\Sigma ,\Delta )$$

$$T$$

$$F(x_1,\dots ,x_p;{\alpha }_1,\dots ,{\alpha }_n,{\Sigma }_1,\dots ,{\Sigma }_n,{\Delta }_1,\dots ,{\Delta }_n)$$

$$\prod _{k=1}^n|I_p+{\Sigma }_{k}T{|}^{-{\alpha }_k}\mathrm{etr}(-{\Sigma }_{k}T(I+{\Sigma }_{k}T{)}^{-1}{\Delta }_k),$$

$$(2\pi {)}^{-p}{\int }_{{\mathcal{C}}_p}A(y_1,\dots ,y_p)B(y_1^{-1},\dots ,y_p^{-1})d{\phi }_1\dots d{\phi }_p$$

$$\hat{f}(t_1,\dots ,t_p)$$

$$\sum b(n_1,\dots ,n_p)\prod _{j=1}^p{z}_j^{n_j}$$

$$\displaystyle|g_{jn}(x)|\leq n^{c}k(x)\mbox{ with a constant $c$ and}$$

$$\displaystyle\int^{\infty}_{0}k(x)e^{-tx}dx<\infty\mbox{ for all $t>0$ .}$$

$$g_j(x,y)=\sum _{n=0}^{\infty }{g}_{jn}(x)y^n,j=1,\dots ,p,$$

$${\hat{g}}_j(t,y)=\frac{{\hat{g}}_{j_0}(t)}{1-y{h}_j(t)},t\ge 0.$$

$$f(x_1,\dots ,x_p)=(2\pi {)}^{-p}{\int }_{{\mathcal{C}}_p}B(y_1^{-1},\dots ,y_p^{-1})\prod _{j=1}^p{g}_j(x_j,y_j)d{\phi }_j$$

$$\displaystyle\begin{array}[]{l}(2\pi)^{-p}\displaystyle\int_{\mathcal{C}_{p}}\left(\sum b(m_{1},\ldots,m_{p})\prod^{p}_{j=1}y^{-m_{j}}_{j}\right)\left(\sum\prod^{p}_{j=1}g_{jn_{j}}(x_{j})y_{j}^{n_{j}}\right)d\varphi_{1}\ldots d\varphi_{p}\\ =\sum b(n_{1},\ldots,n_{p})\prod^{p}_{j=1}g_{jn_{j}}(x_{j})\end{array}$$

$$G_j(x_j,y_j)={\int }_0^{x_j}{g}_j(\xi ,y_j)d\xi$$

$$F(x_1,\dots ,x_p)={\int }_0^{x_1}\dots {\int }_0^{x_p}f({\xi }_1,\dots ,{\xi }_p)d{\xi }_1\dots d{\xi }_p.$$

$$\underset{t\to \infty }{\lim }{h}_j(t)=0$$

$$B(u_1,\dots ,u_p)=\frac{\hat{f}\left(h_1^{-1}(u_1),\dots ,h_p^{-1}(u_p)\right)}{{\prod }_{j=1}^p{\hat{g}}_{j_0}\left(h_j^{-1}(u_j)\right)}$$

$$z_j=(1+v^{-1}{t}_j{)}^{-1},t_j\ge 0,u_j=1-z_j=v^{-1}{t}_j{z}_j,{\omega }_j=z_j-u_j,$$

$$Z=\mathrm{diag}(z_1,\dots ,z_p),U=\mathrm{diag}(u_1,\dots ,u_p),\Omega =\mathrm{diag}({\omega }_1,\dots ,{\omega }_p).$$

$$v^{-1}T=U{Z}^{-1},I_p=Z+U,\Omega =Z-U,$$

$$I+\Sigma T=I+v\Sigma U{Z}^{-1}=(Z+v\Sigma U)Z^{-1}$$

$$\displaystyle\hskip 56.9055ptZ+v\Sigma U=\left\{\begin{array}[]{l@{\hspace{2.2cm}}r}I+(v\Sigma-I)U,\hfil\hskip 62.59596pt&(20a)\\ v\Sigma(I+(v^{-1}\Sigma^{-1}-I)Z),\hfil\hskip 62.59596pt&(20b)\\ \frac{1}{2}(I+v\Sigma)(I+(2(I+v\Sigma)^{-1}-I)\Omega),\hfil\hskip 62.59596pt&(20c)\end{array}\right.$$

$$|I+\Sigma T{|}^{-\alpha }=c^{\alpha }|Z{|}^{\alpha }|I+BY{|}^{-\alpha }$$

$$\displaystyle\hskip 71.13188pt\begin{array}[]{lll@{\hspace{2.6cm}}r}Y=U,&B=v\Sigma-I,&c=1,\hfil\hskip 73.97733pt&(23a)\\ Y=Z,&B=(v\Sigma)^{-1}-I,&c=|I+B|,\hfil\hskip 73.97733pt&(23b)\\ Y=\Omega,&B=2(I+v\Sigma)^{-1}-I,&c=|I+B|.\hfil\hskip 73.97733pt&(23c)\end{array}$$

$$\hat{f}(t_1,\dots ,t_p;\alpha ,\Sigma ,\Delta )=$$

$$\displaystyle\left\{\begin{array}[]{r@{\hspace{1.7cm}}r}|Z|^{\alpha}|I+BU|^{-\alpha}\mathrm{etr}(-(I+B)U(I+BU)^{-1}\Delta),\hskip 48.36958pt&(25a)\\ |I+B|^{\alpha}\mathrm{etr}(-\Delta)|Z|^{\alpha}|I+BZ|^{-\alpha}\mathrm{etr}(Z(I+BZ)^{-1}(I+B)\Delta),\hskip 48.36958pt&(25b)\\ |I+B|^{\alpha}\mathrm{etr}(-\frac{1}{2}\Delta(I-B))\hskip 48.36958pt&(25c)\\ \cdot|Z|^{\alpha}|I+B\Omega|^{-\alpha}\mathrm{etr}(\frac{1}{2}\Omega(I+B\Omega)^{-1}(I+B)\Delta(I-B)),\hskip 48.36958pt&\end{array}\right.$$

$$\displaystyle\hskip 28.45274pt\begin{array}[]{l@{\hspace{3.5cm}}l@{\hspace{2cm}}l@{\hspace{2.2cm}}r}z^{\alpha}u^{n}\hfil\hskip 99.58464pt&vg_{\alpha+n}^{(n)}(vx)\hfil\hskip 56.9055pt&G_{\alpha+n}^{(n)}(vx)\hfil\hskip 62.59596pt&(27a)\\ z^{\alpha+n}\hfil\hskip 99.58464pt&vg_{\alpha+n}(vx)\hfil\hskip 56.9055pt&G_{\alpha+n}(vx)\hfil\hskip 62.59596pt&(27b)\\ z^{\alpha}\omega^{n}\hfil\hskip 99.58464pt&vh_{\alpha,n}(vx)\hfil\hskip 56.9055pt&H_{\alpha,n}(vx)\hfil\hskip 62.59596pt&(27c)\end{array}$$

$$g_{\alpha +n}^{(n)}(x)=\frac{d^n}{d{x}^n}{g}_{\alpha +n}(x)={\left(\genfrac{}{}{0ex}{}{\alpha -1+n}{n}\right)}^{-1}{L}_n^{(\alpha -1)}(x)g_{\alpha }(x)$$

$$h_{\alpha ,n}(x)=(-1{)}^n{\left(\genfrac{}{}{0ex}{}{\alpha -1+n}{n}\right)}^{-1}{L}_n^{(\alpha -1)}(2x)g_{\alpha }(x).$$

$$\displaystyle\left|g_{\alpha+n}^{(n)}(x)\right|\leq\left\{\begin{array}[]{l@{\hspace{.5cm}}l}e^{x/2}g_{\alpha}(x),\hfil\hskip 14.22636pt&\alpha\geq 1\\ 2n\alpha^{-1}e^{x/2}g_{\alpha}(x),\hfil\hskip 14.22636pt&0<\alpha<1\end{array}\right\}$$

$$\displaystyle\left|h_{\alpha,n}(x)\right|\leq\left\{\begin{array}[]{l@{\hspace{.5cm}}l}x^{\alpha-1}/\Gamma(\alpha),\hfil\hskip 14.22636pt&\alpha\geq 1\\ 2nx^{\alpha-1}/\Gamma(\alpha+1),\hfil\hskip 14.22636pt&0<\alpha<1\end{array}\right\},$$

$$\displaystyle\hskip 28.45274ptF_{\alpha}(x,y)=\left\{\begin{array}[]{l@{\hspace{.5cm}}l@{\hspace{2cm}}r}\sum^{\infty}_{n=0}G^{(n)}_{\alpha+n}(x)y^{n}=\frac{1}{1-y}\ G_{\alpha}\left(x,\frac{y}{y-1}\right),\hfil\hskip 14.22636pt&|y|<1,\hfil\hskip 56.9055pt&(35a)\\ \sum^{\infty}_{n=0}G_{\alpha+n}(x)y^{n}=G_{\alpha}(x,y),\hfil\hskip 14.22636pt&y\in\mathbb{C},\hfil\hskip 56.9055pt&(35b)\\ \sum^{\infty}_{n=0}H_{\alpha,n}(x)y^{n}=\frac{1}{1+y}\ G_{\alpha}\left(x,\frac{2y}{y+1}\right),\hfil\hskip 14.22636pt&|y|<1\hfil\hskip 56.9055pt&(35c)\end{array}\right.$$

Full text

1. Introduction

The following notations are used: $\sum_{(n)}$ stands for $\sum_{n_{1}+\ldots+n_{p}=n}$ with $n_{1},\ldots,n_{p}\in\mathbb{N}_{0}$ and $\sum$ without any indices means $\sum^{\infty}_{n=0}\sum_{(n)}$. The notation $D\geq 0$ is also used for non–symmetrical matrices $D_{p\times p}$ with only non–negative eigenvalues. The spectral norm of a $p\times p$–matrix $B$ is denoted by $\|B\|$, $I$ or $I_{p}$ is always an identity matrix and $\mathcal{C}_{p}$ is the $p$–cube $(-\pi,\pi]^{p}$.

The Laplace transform (L.t.) of a $p$–variate non–central $\Gamma_{p}(\alpha,\Sigma,\Delta)$–density with $\alpha>0$, $\Sigma>0$ and a non–centrality matrix $\Delta\geq 0$ was originally obtained from the L.t. of a non-central $W_{p}(2\alpha,\Sigma,\Delta)$–Wishart distribution (with an additional scale factor 2) and is given by

[TABLE]

This function $\widehat{f}$ is generally the L.t. of the density of a real measure on $(0,\infty)^{p}$ which is not always a probability measure. The term ”$\Gamma_{p}(\alpha,\Sigma,\Delta)$–distribution” is used here in this general sense. The exact set of values $\alpha$, leading to a probability density (pdf) $f(x_{1},\ldots,x_{p};\alpha,\Sigma,\Delta)$, depends on $\Sigma$ and presumedly on $\Delta$. To obtain a pdf, all positive integers $2\alpha$ (degrees of freedom) are admissible and all $2\alpha>p-1$. Moreover, in the central case all non–integer values $2\alpha>p-2\geq 0$ are allowed. For $p-2<2\alpha<p-1$ see Royen (1997). Furthermore all $\alpha>0$ are admissible if $|I+\Sigma T|^{-1}$ is infinitely divisible. Two characterizations of infinite divisibility of a $\Gamma_{p}(\alpha,\Sigma)$–distribution are found in Griffiths (1984) and Bapat (1989). Further conditions for admissible non–integer $2\alpha<p-2$ are given in Royen (1997), (2006).

Three integral representations by integration over $\mathcal{C}_{p}$ are provided by theorem 2 in section 4 for the functions

[TABLE]

where $f$ has the L.t.

[TABLE]

Thus, $F$ is not always the cumulative distribution function (cdf) of a probability measure.

In particular let $X_{p\times n}$ be a $\mathcal{N}_{p\times n}(M_{p\times n},\Sigma_{p\times p}\otimes I_{n})$–random matrix and $A_{n\times n}\geq 0$ of rank $q$ with $T^{\prime}AT=\Lambda=\mathrm{diag}(\lambda_{1},\ldots,\lambda_{n})$, $\lambda_{1}\geq\ldots\geq\lambda_{n}$. Then the joint distribution of the diagonal elements of the generalized quadratic form $\frac{1}{2}XAX^{\prime}$ equals the distribution of the diagonal of $\frac{1}{2}Y\Lambda Y^{\prime}$ with a $\mathcal{N}_{p\times n}(MT,\Sigma_{p\times p}\otimes I_{n})$–distributed $Y=XT$. This is the distribution of a sum of $q$ independent $\Gamma_{p}(\frac{1}{2},\lambda_{k}\Sigma,\Delta_{k}=\frac{1}{2}\mu_{k}^{*}\mu_{k}^{*^{\prime}}\Sigma^{-1})$–random vectors, where $\mu_{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 $\alpha_{k}=\frac{1}{2}$, $\Sigma_{k}=\lambda_{k}\Sigma$, $k=1,\ldots,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 $\chi^{2}_{p}(2\alpha,\Sigma)$–cdf were derived from three different representations of the $\chi^{2}_{p}(2\alpha,\Sigma)$–L.t. which are extended to the general $\Gamma_{p}(\alpha,\Sigma,\Delta)$–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 $\mathcal{C}_{p}$ is explained in more detail in section 2 and summarized in theorem 1. The idea is as follows:

If $A(z_{1},\ldots,z_{p})$ and $B(z_{1},\ldots,z_{p})$ are analytical functions whose power series have the coefficients $a(m_{1},\ldots,m_{p})$ and $b(n_{1},\ldots,n_{p})$ and which are absolutely convergent for $\max|z_{j}|<r_{A}$ and $\max|z_{j}|<r_{B}$ respectively, where $r^{-1}_{B}<r_{A}$, then

[TABLE]

holds with $y_{j}=re^{i\varphi_{j}}$, $-\pi<\varphi_{j}\leq\pi$, $j=1,\ldots,p$ and $r^{-1}_{B}<r<r_{A}$.

The integrals in (1.) 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 (1.). 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 $\mathcal{C}_{p}$ and not over $\mathbb{R}^{p}$ as by the Fourier or Laplace inversion formula. A single $\Gamma_{p}(\alpha,\Sigma,\Delta)$–cdf is represented by a $(p-1)$–variate integral over $\mathcal{C}_{p-1}$ in section 5.

A totally different ${m+1\choose 2}$–variate integral representation of the $\Gamma_{p}(\alpha,\Sigma)$–cdf has been given recently by Royen (2006), which is based on $m$–factorial decompositions $\sum^{-1}_{p\times p}=D-BB^{\prime}$, where $D$ is a real or complex diagonal matrix minimizing the rank $m$ of $\Sigma^{-1}-D$. Approximations to a $\Gamma_{p}(\alpha,\Sigma)$–cdf are obtained by $m$–factorial approximations to $\Sigma$ 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 transforms, but the version below is sufficient for the purpose of the underlying paper.

Let $\widehat{f}(t_{1},\ldots,t_{p})$, $t_{1},\ldots,t_{p}\geq 0$, be a given L.t. of an unknown function $f(x_{1},\ldots,x_{p})$ with $f=0$ for $\min x_{j}<0$. It is assumed that there are univariate L.t. $\widehat{g}_{j_{0}}(t)$ of some probability densities $g_{j_{0}}(x)$ on $(0,\infty)$ and further functions $h_{j}(t)$ with $|h_{j}(t)|\leq 1$, uniformly for $t\geq 0$, which enable a representation

[TABLE]

with an analytical function $B(z_{1},\ldots,z_{p})$ whose power series expansion

[TABLE]

is absolutely convergent for $|z_{1}|,\ldots,|z_{p}|<r_{B}$ with a certain value $r_{B}>1$.

Furthermore, the products $\widehat{g}_{j_{0}}(t)(h_{j}(t))^{n}$ are supposed to be the L.t. of continuous functions $g_{jn}(x)$, $x>0$, which satisfy the conditions

[TABLE]

Hence, the generating functions (generators)

[TABLE]

are defined for all $x>0$ and $|y|<1$, and they have the L.t.

[TABLE]

Theorem 1

Under the assumptions from (6) and (2.) $\widehat{f}$ in (5) is the L.t. of

[TABLE]

with $y_{j}=re^{i\varphi_{j}}$,  $-\pi<\varphi\leq\pi$,  $r^{-1}_{B}<r<1$,  $g_{j}$ from (8).

• Proof. 

The integral in (10) is evaluated by

[TABLE]

and this series has the L.t. from (5).

Some further remarks: With

[TABLE]

instead of the $g_{j}$ in (10), the corresponding representation arises for

[TABLE]

If the series in (8) are absolutely convergent for all $y\in\mathbb{C}$ then additionally

[TABLE]

is supposed to hold. Then the rhs of (9) is the L.t. of $g_{j}(x,y)$ for any fixed $y$ and all sufficiently large $t$.

In some cases the functions $g_{j_{0}}$ and their L.t. $\widehat{g}_{j_{0}}$ are known from univariate marginal distributions apart from some scale factors. If the functions $u_{j}=h_{j}(t)$ are explicitly invertible then

[TABLE]

can sometimes be found easily from the given $\widehat{f}$.

3. Three representations for the $\Gamma_{p}(\alpha,\Sigma,\Delta)$–Laplace transform

and the related generators

With any $v>0$ we define

[TABLE]

The scale factor $v$ is introduced to obtain $\|B\|<1$ for the matrices $B$ defined in (3.) 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

[TABLE]

it follows for the matrices $I+\Sigma T$ in the L.t. (1):

[TABLE]

and

[TABLE]

and therefore

[TABLE]

with

[TABLE]

It should be noticed that $\|B\|<1$ in (3.c) for every $v>0$ and $\Sigma>0$.

Now, using (17), by a straightforward calculation the L.t. in (1) can be represented by

[TABLE]

with the corresponding matrices $B$ from (3.) and $Z,U,\Omega$ from (3.).

For the former series expansions the following relations were used:

Laplace transform $\widehat{f}(t)$:             $f(x)$:               $F(x)\,=\,\int_{0}^{x}f(\xi)d\xi$:

[TABLE]

where $z=(1+v^{-1}t)^{-1}$, $g_{\alpha+n}(x)=e^{-x}x^{\alpha-1+n}/\Gamma(\alpha+n)$,

[TABLE]

with the generalized Laguerre polynomials $L^{(\alpha-1)}_{n}$ and

[TABLE]

The last identity is verified by L.t.

The following bounds are derived from (22.14.13) in Abramowitz and Stegun (1965):

[TABLE]

[TABLE]

matching with the conditions in (2.).

The following generators (generating functions) with the $\Gamma(\alpha+n)$–cdf $G_{\alpha+n}(x)$ are required for the formulas in theorem 2:

[TABLE]

The identities (a) and (c) are verified by the L.t. of $f_{\alpha}(x,y)=\frac{\partial}{\partial x}F_{\alpha}(x,y)$. A short calculation shows

[TABLE]

and

[TABLE]

The functions $F_{\alpha}(x,y)$ are especially simple for $\alpha\in\mathbb{N}$ since $G_{\alpha}(z)=1-e^{-z}\sum^{\alpha-1}_{j=0}{z^{j}}/{j!}$, $\alpha\in\mathbb{N}$.

Besides,

[TABLE]

The following simple lemma is used for the proof of theorem 2.

Lemma 1

*If $B$ is a symmetrical $p\times p$–matrix with $\|B\|<1$ and

$Y=\mathop{\rm diag}\nolimits(y_{1},\ldots,y_{p})$ then the power series expansion*

[TABLE]

is absolutely convergent for $\max|y_{j}|<r_{B}=\|B\|^{-1}$.

This follows from $\sum_{(n)}|b(n_{1},\ldots,n_{p})|=O(\vartheta^{n})$ with any $\vartheta>\|B\|$, which has been already shown in (2.1.16) $\ldots$ (2.1.18) in Royen (1991) (with the notation $-C$ instead of $B$).

4. The integral representations

In theorem 2 below the functions $F(x_{1},\ldots,x_{p};\alpha_{1},\ldots,\alpha_{n},\Sigma_{1},\ldots,\Sigma_{n},\Delta_{1},\ldots,\Delta_{n})$ from (1.) are represented by three different integrals over $\mathcal{C}_{p}=(-\pi,\pi]^{p}$. Together with the generators $F_{\alpha}$ from (3.), $\alpha=\sum^{n}_{k=1}\alpha_{k}$, the following matrices are used with a scale factor $v$ to enforce $\|B_{k}\|<1$:

[TABLE]

Furthermore, we define $\lambda_{\max}=\max\|\Sigma_{k}\|$ ,  $\lambda_{\min}^{-1}=\max\|\Sigma_{k}^{-1}\|$,  $y_{j}=re^{i\varphi_{j}}$,

$-\pi<\varphi_{j}\leq\pi$,  $Y=\mathop{\rm diag}\nolimits(y_{1},\ldots,y_{p})$,

[TABLE]

where the negative sign occurs only with $B_{k},D_{k}$ from (4.a), and

$\mathcal{F}_{\alpha}d\varphi=\prod^{p}_{j=1}F_{\alpha}(vx_{j},y_{j})d\varphi_{j}$.

Theorem 2

With the above notations the functions $F$ from (1.) are respresentable by each of the following three integrals:

[TABLE]

$F_{\alpha}$* from (3.a), $B_{k},D_{k}$ from (4.a), $\|B_{k}\|<1$ if $v<2\lambda^{-1}_{\max}$, $\max\|B_{k}\|<r<1$,*

[TABLE]

$F_{\alpha}$* from (3.b), $B_{k},D_{k}$ from (4.b), $\|B_{k}\|<1$ if $v>\frac{1}{2}\lambda^{-1}_{\min}$, $\max\|B_{k}\|<r$,*

[TABLE]

$F_{\alpha}$* from (3.c), $B_{k},D_{k}$ from (4.c), $v>0$, $\max\|B_{k}\|<r<1$.*

• Proof. 

Because of lemma 1 the assumptions of theorem 1 are satisfied with $\widehat{g}_{j_{0}}(t)=z^{\alpha}_{j}=(1+v^{-1}t_{j})^{-\alpha}$ and $h_{j}(t_{j})$ corresponding to $z_{j}$ or $u_{j}=v^{-1}t_{j}z_{j}=1-z_{j}$ or $\omega_{j}=z_{j}-u_{j}$ respectively. The functions $\widehat{g}_{j_{0}}(t)(h_{j}(t))^{n}$ are the L.t. of the functions in the second column of (3.) from which type (a) and (c) have the bounds in (30), (33), satisfying the condition (2.) for theorem 1. The series $\sum^{\infty}_{n=0}G_{\alpha+n}(x)y^{n}=G_{\alpha}(x,y)$ in (3.b) is absolutely convergent for every $y\in\mathbb{C}$. Thus, all $r>\max\|B_{k}\|$ are admissible in (44). In (45) we have $\max\|B_{k}\|<1$ for every $v>0$. Hence, theorem 1 together with the respresentations of the L.t. in (3.) implies (43), (44) and (45).

The univariate case of (44) provides

[TABLE]

with $2v>\max\sigma^{-2}_{k}$, $r=1$, $G_{\alpha}$ from (39). With $p=1$ similar formulas arise from (43) or (45).

The cdf of a quadratic form $\frac{1}{2}x^{\prime}Ax$ with $T^{\prime}AT=\mathop{\rm diag}\nolimits(\lambda_{1},\ldots,\lambda_{n})\geq 0$ of rank $q$ and a $\mathcal{N}(\mu,\sigma^{2}I_{n})$–random vector $x$ is a special case of (4.) with $\alpha_{k}=\frac{1}{2}$, $\sigma^{2}_{k}=\lambda_{k}\sigma^{2}$ and non–centrality parameters $\delta^{2}_{k}=\frac{1}{2}\mu^{*2}_{k}/\sigma^{2}$, $k=1,\ldots,q$, $\mu^{*}=T^{\prime}\mu$.

Some further remarks: In (44) also $\|B_{k}\|>1$ is allowed since every $r=\|Y\|>\max\|B_{k}\|$ is admissible, which entails $\max\|B_{k}Y^{-1}\|<1$.

With $\vartheta=\lambda_{\max}/\lambda_{\min}$ it follows with special values of $v$:

[TABLE]

but

[TABLE]

More generally, the scale factor $v=w^{2}$ can be replaced by a scale matrix $W^{2}=\mathop{\rm diag}\nolimits(w^{2}_{1},\ldots,w^{2}_{p})>0$. Then with $T_{w}=W^{-1}TW^{-1}$, $\Sigma_{w}=W\Sigma W$, $\Delta_{w}=W\Delta W^{-1}$ the L.t. (1) equals

[TABLE]

Consequently, besides the substitutions $v\Sigma_{k}\to W\Sigma_{k}W$, $\Delta_{k}\to W\Delta_{k}W^{-1}$, the matrices $I+B_{k}$ in theorem 2 must be replaced by $W\Sigma_{k}W$, $(W\Sigma_{k}W)^{-1}$ and $2(I+W\Sigma_{k}W)^{-1}$ respectively, and the generators $F_{\alpha}(vx_{j},y_{j})$ by $F_{\alpha}(w_{j}^{2}x_{j},y_{j})$.

In particular for a single $\Gamma_{p}(\alpha,\Sigma,\Delta)$–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 (43), whereas this condition is satisfied in (45) for every scaling. It was shown in Royen (1991) that natural scaling can always be accomplished also in $I+B=2(I+W\Sigma W)^{-1}$ by a unique $W^{2}$.

5. Representations of the $\Gamma_{p}(\alpha,\Sigma,\Delta)$ distribution function

by $(p-1)$–variate integrals

For a single $\Gamma_{p}(\alpha,\Sigma,\Delta)$–cdf it is always possible to perform the integration over a single variable $\varphi_{j}$ within the integrals from theorem 2.

We use the following functions

[TABLE]

For positive half integers $\alpha=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 $I_{j-1/2}(2(xy)^{1/2})$, $j=1,\ldots,k$, (see e.g. Royen (1995) or (2006)).

Now let be $W^{2}=\mathop{\rm diag}\nolimits(w^{2}_{1},\ldots,w^{2}_{p})$ a general scale matrix,

$Y=\mathop{\rm diag}\nolimits(y_{1},\ldots,y_{p})$, $y_{j}=re^{i\varphi_{j}}$, $-\pi<\varphi_{j}\leq\pi$,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where the negative sign is only taken for $B_{pp}$ from (5.a).

Theorem 3

With the above notations the $\Gamma_{p}(\alpha,\Sigma,\Delta)$–cdf $F(x_{1},\ldots,x_{p};\alpha,\Sigma,\Delta)$ is given by each of the following three integrals:

[TABLE]

$B$* from (5.a), $D$ from (5.a), $\|B\|<r<1$,*

[TABLE]

$B$* from (5.b), $D$ from (5.b), $\|B\|<r$,*

[TABLE]

$B$* from (5.c), $D$ from (5.c), $\|B\|<r<1$.*

For the proof of theorem 3 the following two lemmas are required.

Lemma 2

With $Y=\mathop{\rm diag}\nolimits(y_{1},\ldots,y_{p})$, $y_{j}=re^{i\varphi_{j}}$, $\|B\|<r$, $B,D,y_{0},q$ from (5.), (5.), (56), (57) the following decomposition is obtained

[TABLE]

• Proof. 

From frequently used formulas for $p\times p$–matrices, (see e.g. complements and problems 2.4, 2.7 in chapter 1b of Rao (1973)) it follows for

$A=Y+B=\scriptstyle{\begin{pmatrix}A_{pp}&b_{p}\\ b^{\prime}_{p}&\qquad y_{p}+b_{pp}\end{pmatrix}}$:

[TABLE]

and

[TABLE]

Lemma 3

Let be $q$ any number, S_{r}=\{y\in\mathbb{C}\big{|}|y|=r\}, $y_{0}$ any number with $|y_{0}|<r$, then with $F_{\alpha}$ from (3.), $\mathcal{G}_{\alpha},\mathcal{G}^{*}_{\alpha}$ from (51), and the negative sign in $\pm q$ only for (3a)

[TABLE]

• Proof. 

It is sufficient to verify (3) for the corresponding derivatives $f_{\alpha}=\frac{\partial}{\partial x}F_{\alpha}$. At first, (3a) is shown:

With $F_{\alpha}$ from (3.a) and the binomial series for $(1-y_{0}/y)^{-(\alpha+n)}$ we obtain

[TABLE]

With $z=(1+t)^{-1}$, $u=tz$, the last integral has the L.t.

[TABLE]

Multiplication by $(-q)^{n}/n!$ and summation over $n$ leads to the L.t.

[TABLE]

and this is the L.t. of $\frac{\partial}{\partial x}\ \mathcal{G}_{\alpha}\left(\frac{x}{1-y_{0}},\frac{q}{1-y_{0}}\right)$.

To verify (3b) we obtain with $F_{\alpha}$ from (3.b):

[TABLE]

Multiplication by $q^{n}/n!$ and summation provides $(1-y_{0})^{-\alpha}\frac{\partial}{\partial x}\ \mathcal{G}^{*}\left((1-y_{0})x,\frac{q}{1-y_{0}}\right)$.

(3c) can be shown by L.t. in a similar way as (3a).

Proof of theorem 3.

Without loss of generality $y_{p}$ is selected from the variables $y_{j}=re^{i\varphi_{j}}$ in $Y=diag(y_{1},\ldots,y_{p})$ with any fixed $r>\|B\|$. If $y_{p}$ is replaced by a variable $y$ with any $|y|$ then the equation

[TABLE]

has always a unique solution

[TABLE]

with $|y_{0}|<r$ since $\|B_{pp}\|\leq\|B\|$.

Hence, with lemma 2 and lemma 3, theorem 3 is obtained by integration over $\varphi_{p}$ in the integrals of theorem 2 with $n=1$.