# A branch cut approach to the probability density and distribution functions
of a linear combination of central and non-central Chi-square random variables

Alfred Kume<EMAIL_ADDRESS>Tomonari Sei<EMAIL_ADDRESS>Andrew
T. A. Wood<EMAIL_ADDRESS>
###### Abstract

The paper considers the distribution of a general linear combination of
central and non-central chi-square random variables by exploring the branch
cut regions that appear in the standard Laplace inversion process. Due to the
original interest from the directional statistics, the focus of this paper is
on the density function of such distributions and not on their cumulative
distribution function. In fact, our results confirm that the latter is a
special case of the former. Our approach provides new insight by generating
alternative characterizations of the probability density function in terms of
a finite number of feasible univariate integrals. In particular, the central
cases seem to allow an interesting representation in terms of the branch cuts,
while general degrees of freedom and non-centrality can be easily adopted
using recursive differentiation. Numerical results confirm that the proposed
approach works well while more transparency and therefore easier control in
the accuracy is ensured.

###### keywords:

Linear combination of Chi-squares, Bingham distributions, Fisher-Bingham
distributions, directional statistics, holonomic functions.

organization=SMSAS , addressline=University of Kent, city=Canterbury,
country=U.K.

organization=The University of Tokyo , city= Tokyo, country=Japan

organization=Research School of Finance, Actuarial Studies & Statistics ,
addressline=Australian National University, country=Australia

##  1 Introduction

Let us denote by $X$ a quadratic form of normally distributed components:

$X=\sum_{i=1}^{p}Y_{i}^{\top}Y_{i}=\sum_{i=1}^{p}\lambda_{i}\chi^{2}_{n_{i}}(\delta_{i})\quad\lambda_{i}=\frac{1}{2\theta_{i}}\quad\delta_{i}=\frac{\gamma_{i}^{2}}{2\theta_{i}}$
(1)

where $Y_{i}\sim
N_{n_{i}}(\frac{\gamma_{i}}{2\theta_{i}\sqrt{n_{i}}}\textbf{1},\frac{1}{2\theta_{i}}\textbf{I}_{n_{i}})$
are $p$ independent multivariate normal rv’s of dimension $n_{i}$, each with
covariance matrix the multiple of identity
$\frac{1}{2\theta_{i}}\textbf{I}_{n_{i}}$ and mean some vector of norm
$\frac{\gamma_{i}}{2\theta_{i}}$.

Evaluating the probability density function (pdf) and the cumulative
distribution function (cdf) of such quadratic forms has important applications
in statistical theory. Many statistical tests such as the goodness of fit for
non-parametric regression (e.g. Kuonen (2005)), pseudo-likelihood ratio tests
(e.g. Liang and Self (1996)), directional statistics models using the Fisher-
Bingham distributions (see Mardia and Jupp (2000); Kume and Wood (2005)) or
the shape distributions based on the complex Bingham family (see Kent
(1994a)), are heavily relying on data modelled by these types of random
variables. Since in their most general cases there is no close form expression
for the corresponding distributions the numerical evaluation is implemented in
many situations. Representation (1) is implicitly assuming that $\theta_{i}$
as well as $\lambda_{i}$ are positive. The more general cases when some
$\theta_{i}$’s are negative can be seen as an extension of (1) since by
collecting the terms according to their signs, it can be considered as a
difference of two such positive combinations. We will consider these cases
separately.

Standard inverting arguments of the corresponding pdf’s via the moment
generating function or Laplace transform confirm that the density function of
$X$ is

$\displaystyle f_{\bm{\theta},\bm{\gamma}^{2},\bm{n}}(s)$ $\displaystyle=$
$\displaystyle\prod_{i=1}^{p}\theta_{i}^{\frac{n_{i}}{2}}\exp(-\sum_{i=1}^{p}\frac{\gamma^{2}_{i}}{4\theta_{i}})\frac{1}{2\pi\textbf{i}}\int\limits_{\textbf{i}\mathbb{R}+t_{0}}\frac{\text{exp}(\sum_{i=1}^{p}\frac{\gamma_{i}^{2}}{4(\theta_{i}+t)})}{\prod_{i=1}^{p}(\theta_{i}+t)^{\frac{n_{i}}{2}}}e^{st}dt,\quad
s>0$ (2)

where $t_{0}$ can be any real number such that the vertical line
${\textbf{i}\mathbb{R}+t_{0}}$ in the complex plane is on the right of the
poles $-\theta_{i}$ i.e. $-\min(\bm{\theta})<t_{0}$. Note that throughout the
paper we refer to $\bm{\theta},\bm{\gamma}^{2},\bm{n}$ as the vectorised
entries of $\theta_{i}$, $\gamma_{i}$ and $n_{i}$ respectively where
$\theta_{i}$ are forced to be distinct; and we apply the usual convention for
complex integration so that for the closed contours we integrate along the
anticlockwise direction and for the unbounded lines such as
${\textbf{i}\mathbb{R}+t_{0}}$ we integrate from $t_{0}-{\textbf{i}\infty}$ to
$t_{0}+{\textbf{i}\infty}$. The current methods addressing this inversion
problem could be categorised into three groups:

Standard numerical inversion: Attempts to numerically evaluate these functions
go back as far as in Imhof (1961) where a general numerical integrating
procedure along a vertical line similar to that in (2) is shown. The
alternative approach introduced in Davies (1973) and Davies (1980) is
addressing the evaluation of the cumulative distribution function calculation
with numerical performance comparable to that of Imhof (1961). A moment
matching procedure to some gamma distribution is reported in Liu et al.
(2009). However as indicated by Duchesne and De Micheaux (2010) this approach
is not always guaranteeing good performance. All of the above mentioned
methods are generally focused on the cdf’s and not on the pdf
$f_{\bm{\theta},\bm{\gamma}^{2},\bm{n}}$.

Saddle point approximation (SPA): This method relies on the extensions of the
Laplace approximation for integrals (see e.g. Daniels (1954); Barndorff-
Nielsen (1990)). The manuscript of Butler (2007) is an excellent review of
this popular method. Its implementation in our context is straightforward with
little computational cost. This approach is shown in Kuonen (1999) to perform
well with a surprising accuracy in non-parametric regression examples
reported. The same conclusion is reported by Kume and Wood (2005) in the
context of the density function evaluation. The authors implemented this for
directional statistics inference involving the evaluation of
$f_{\bm{\theta},\bm{\gamma}^{2},\bm{n}}$ at point $s=1$.

Holonomic gradient methods: The evaluation of
$f_{\bm{\theta},\bm{\gamma}^{2},\bm{n}}$ is in fact characterized as a
solution of a particular ODE equation (c.f. Nakayama et al. (2011)) where in
principle only some accurate starting point and a good ODE solver is required.
The derivation of such ODE derives from the fact that
$f_{\bm{\theta},\bm{\gamma}^{2},\bm{n}}$ is shown to be an holonomic function
for which the Pfaffian matrix equations lead to the relevant ODE. The general
adjustments for allowing multiplicities in $\bm{\theta}$ so that the
corresponding ODE do not become degenerate, are studied in Sei and Kume
(2013); Kume and Sei (2018) for central and non central cases respectively. In
fact Koyama and Takemura (2016) focus on the problem of using HGM for the cdf.

This paper belongs to the first group of these approaches. We focus on
exploring the specific structure of the Laplace transform so that a more
flexible approach both numerically and algebraically is offered. More
specifically, the main contribution aspects of this paper are:

  * 1.

Theorem 1 expresses the integrating contour in a much simpler form as seen
before in the literature, leading to numerically simpler inversion procedures.
This is achieved by exploring the simple branch cut structure due to the
square-root functions. Additionally, we establish a close connection between
the cdf of $X$ as in (1) and the pdf of some other positive linear combination
of non-central Chi-squares $X^{\prime}$ such that

$P(X\leq
x)=f_{\tilde{X}}(x)\frac{e^{x\theta_{0}}}{\theta_{0}}\prod_{i=1}^{p}\left(\frac{\theta_{i}}{\theta_{i}+\theta_{0}}\right)^{\frac{n_{i}}{2}}\text{exp}(-\sum_{i=1}^{p}\frac{\gamma^{2}_{i}\theta_{0}}{4\theta_{i}(\theta_{i}+\theta_{0})})\quad\forall\theta_{0}>0$
(3)

where $f_{\tilde{X}}$ is the pdf of
$\tilde{X}=\frac{1}{2\theta_{0}}\chi^{2}_{2}(0)+\sum_{i=1}^{p}\frac{1}{2(\theta_{i}+\theta_{0})}\chi^{2}_{n_{i}}(\frac{\gamma_{i}^{2}}{2(\theta_{i}+\theta_{0})})$
for any arbitrary $\theta_{0}$. An immediate consequence of this fact is that
any general method for evaluating the pdf of any positive linear combination
of non-central Chi-square random variables suffices for the distribution
functions too. For example the saddlepoint approximation of pdf’s as in Kume
and Wood (2005) can also be used for the tail probabilities as an alternative
to that of Kuonen (1999); Barndorff-Nielsen (1990). Additionally, the ODE
equations of Koyama and Takemura (2016) for cdf. could also be seen as a
special case of those developed in Sei and Kume (2013); Kume and Sei (2018)
for the pdf. To our knowledge property (3) has not been reported in the
literature before.

  * 2.

Theorem 2 makes the corresponding simplifications for the practically
important central cases, $\bm{\gamma}=0$, where the univariate contours are
now reduced to a finite linear combination of some bounded elementary
univariate integrals which are numerically manageable and therefore the
normalising constants of the Bingham distributions can be easily derived.

  * 3.

Theorem 3 extends the above mentioned results to the difference of two
positive linear combinations including the corresponding result as in (3). In
particular, for the case of chi-square distributions of degrees of freedom 2
which corresponds to the rescaled exponential distributions, a closed form
expression is immediately obtained using residues of simple poles.

  * 4.

The relative simplicity and efficiency of our inversion is confirmed through
numerical examples. In particular, we show that we can not just simply offer
individual entries to the tail probabilities like shown before but the whole
function is easily evaluated due to feasibility of the numerical integration.
We have successfully implemented here the standard numerical integration
routines within the R package where, if required, the multiple precision
integration is available calling the GMP library of the C programming
environment.

The Laplace inversion (2) provides additional identities for these
distributions. More specifically, using the flexibility on the choice of such
$t_{0}$ and applying a change of variable $u=s(t+c)$, for any real value of
$c$ implies the rescale and shift property

$f_{\bm{\theta},\bm{\gamma}^{2},\bm{n}}(s)=\frac{e^{-sc}}{s}f_{s\bm{\theta}-sc,s\bm{\gamma}^{2}}(1)\prod_{i=1}^{p}\frac{\theta_{i}^{\frac{n_{i}}{2}}}{(\theta_{i}-c)^{\frac{n_{i}}{2}}}\text{exp}(\sum_{i=1}^{p}\frac{\gamma^{2}_{i}}{4(\theta_{i}-c)}-\frac{\gamma^{2}_{i}}{4\theta_{i}}).$
(4)

In fact, $f_{s\bm{\theta}-sc,s\bm{\gamma}^{2}}(1)$ is also related to that of
evaluating the normalising constant of the Fisher-Bingham distributions on
$\sum_{i=1}^{p}n_{i}-1$ dimensional spheres (c.f. Kume and Wood (2005)):

$\mathcal{C}(\bm{\theta},\bm{\gamma},\bm{n})=\int\limits_{\mathcal{S}^{\sum_{i=1}^{p}n_{i}-1}}\text{exp}(\sum_{i=1}^{p}-\theta_{i}x_{i}^{2}+\gamma_{i}x_{i})d_{\mathcal{S}^{p-1}}(\bm{x})=2\pi^{\frac{\sum_{i=1}^{p}n_{i}}{2}}\frac{\text{exp}(\sum_{i=1}^{p}\frac{\gamma^{2}_{i}}{4\theta_{i}})}{\prod_{i=1}^{p}\theta_{i}^{\frac{n_{i}}{2}}}f_{\bm{\theta},\bm{\gamma}^{2},\bm{n}}(1)$
(5)

where $d_{\mathcal{S}^{\sum_{i=1}^{p}n_{i}-1}}(\bm{x})$ is the uniform measure
and if all $\gamma_{i}=0$,
$\mathcal{C}(\bm{\theta},\bm{0},\bm{n})=\mathcal{C}(\bm{\theta},\bm{n})$
reduces to that of the Bingham distribution. The relationship

$f_{\bm{\theta},\bm{\gamma}^{2},\bm{n}}(s)=\frac{\prod_{i=1}^{p}\theta_{i}^{\frac{n_{i}}{2}}\text{exp}(\sum_{i=1}^{p}-\frac{\gamma^{2}_{i}}{4\theta_{i}})}{2\pi^{\frac{\sum_{i=1}^{p}n_{i}}{2}}}s^{\frac{\sum_{i=1}^{p}n_{i}}{2}-1}\mathcal{C}(s\bm{\theta},\sqrt{s}\bm{\gamma},\bm{n})$
(6)

implies that the methods developed for evaluating the normalizing constants
for directional distributions apply immediately to those of the density
function of some linear combination of chi-squares. In fact the shift property
(4) implies that $\theta_{i}$ in (5) can be allowed to be shifted at negative
values, as it is the case for alternative parametrisations of
$\mathcal{C}(\bm{\theta},\bm{\gamma},\bm{n})$ where $\theta_{i}$ can be either
positive or negative. An additional property which derives from (2) is the
connection between multiplicities (or degrees of freedom) $n_{i}$ and
differentiation:

$f_{\bm{\theta},\bm{\gamma}^{2},\bm{n}}(s)\prod_{i=1}^{p}\theta_{i}^{-\frac{n_{i}}{2}}e^{-\sum_{i=1}^{p}\frac{\gamma^{2}_{i}}{4\theta_{i}}}=\mathcal{D}\mathcal{D^{\prime}}\left(f_{\bm{\theta}^{0},\bm{\gamma}^{2},\bm{n}^{0}}(s)\prod_{i=1}^{p}\theta_{i}^{-\frac{n^{0}_{i}}{2}}\text{exp}(-\sum_{i=1}^{p}\frac{\gamma^{2}_{i}}{4\theta_{i}})\right)$
(7)

where

$\mathcal{D^{\prime}}=\prod_{i=1}^{p_{1}}\frac{(-1)^{r_{i}}}{n^{0}_{i}(n^{0}_{i}+1)\cdots(n^{0}_{i}+r_{i}-1)}\frac{\partial^{r_{i}}}{\partial^{r_{i}}(\theta^{\prime}_{i})}\quad\text{}\quad\mathcal{D}=\prod_{i=1}^{p-p_{1}}\frac{\partial^{r_{i}}}{\partial^{r_{i}}(\gamma_{i}^{2}/4)}$

and $f_{\bm{\theta}^{0},\bm{\gamma}^{2},\bm{n}^{0}}(s)$ representing the pdf
with each $\theta_{i}$ having multiplicity $n^{0}_{i}\in\\{1,2\\}$,
$n_{i}=2r_{i}+n^{0}_{i}$ and the first $p_{1}$ terms are collected to
represent the central cases i.e. $\gamma_{i}=0$.

These differentiation rules generalize similar results of Kume and Wood (2007)
for the normalizing constants of Bingham distributions and those of Imhof
(1961) for the pdf.

While there have been ongoing attempts in the literature (see the recent work
in Chen and Tanaka (2021) and the references therein), to address this problem
numerically, they are mostly focused only in the case of positive
$\theta_{i}$. In fact, the properties reported here imply alternative
computation methods which are easily implemented terms of some elementary
integrating terms and also hold even if in (2) some $\theta_{i}\leq 0$ .
Therefore if extreme accuracy is required, standard multiple precision
integration can be easily adopted to the resulting elementary functions.

The paper is organized as follows. In Section 2 we explore the branch cuts of
the integrand function of (2) and then adjust the contour accordingly. This
gives rise to the connection of the pdf expressions with those for the cdf
stated in (3). In the following section we report the special adaptations for
the relevant quantities for practical implementation. We then extend in
Section 4 the results to that of a difference between two quadratic forms of
type (1). We report some numerical illustrations in Section 5 and conclude
with some general remarks.

##  2 Exploring the branch cut structure

We focus initially on the complex valued function

$h(t)=\frac{1}{\sqrt{(\theta_{1}+t)(\theta_{2}+t)\cdots(\theta_{p}+t)}}=\prod_{i=1}^{p}r_{i}^{-1}e^{\text{-{i}}\alpha}\quad\alpha=\sum_{i=1}^{p}\alpha_{i}$
(8)

which corresponds to the central case, $\bm{\gamma}=\bm{0}$. For each term
above we use the parametrization

$\sqrt{\theta_{i}+t}=r_{i}e^{\text{{i}}\alpha_{i}}\quad\text{where}\quad\alpha_{i}\in(+\pi,-\pi).$

It is important to note that these terms are individually not analytic for
each positive $r_{i}$ and $\alpha_{i}=\pm\pi$ and the corresponding branch cut
is the line of real numbers on the left of $-\theta_{i}$. The product of
square roots however leads to pairwise branch cuts as

$C_{r}=[-\theta_{2r},-\theta_{2r-1}]$

which are segments determined by pairs of $-\theta_{i}$. If $p$ is odd then
the last value $-\theta_{p}$ will be paired with $-\infty$ as

$C^{\infty}_{[\frac{p+1}{2}]}=[-\theta_{p},-\infty].$

One can easily see that in general the product

$\frac{1}{\sqrt{(\theta_{1}+t)(\theta_{2}+t)\cdots(\theta_{p}+t)}}$

has $k=[\frac{p+1}{2}]$ branch cuts

$C_{1},C_{2}\cdots C_{k}$

such that if $p=2k-1$, $C_{k+1}$ is unbounded.

Please note that the choice of these branch cuts is somewhat arbitrary as each
square root term can allow branches on either side of the respective pole
$-\theta_{i}$ and therefore each pair of products
$\sqrt{(\theta_{i}+t)(\theta_{j}+t)}$ can generate either an individual branch
cut as the segment $[\theta_{i},\theta_{j}]$ or its complement in the real
line. More specifically, using a different parametrization, one could
alternatively generate the unbounded branch cut extending from $-\theta_{1}$
to $+\infty$ instead.

We are using here the ordered pairs so that the interpretation is meaningful
such that it allows for the odd number $p$ to generate the leftmost branch
extending to $-\infty$. Since outside these branch cuts the non-central
components are analytical functions except the respective poles $-\theta_{i}$
then the branch cut structure is not affected also for

$g(t)=\prod_{i=1}^{p}\frac{\text{exp}(\sum_{i=1}^{p}\frac{\gamma_{i}^{2}}{4(\theta_{i}+t)})}{(\theta_{i}+t)^{k_{i}}}\frac{1}{\prod_{i=1}^{p}\sqrt{\theta_{i}+t}}=\prod_{i=1}^{p}\frac{\text{exp}(\sum_{i=1}^{p}\frac{\gamma_{i}^{2}}{4(\theta_{i}+t)})}{(\theta_{i}+t)^{k_{i}}}h(t)$

which corresponds also to the cases of odd degrees of freedom, i.e.
$n_{i}=2k_{i}+1$. If we have even multiplicities for some particular
$\theta_{i}$, say $n_{i}=2k_{i}$ then there is no presence in the $h(t)$ part
of the integrand and the contribution of the corresponding term
$\frac{1}{(\theta_{i}+t)^{k_{i}}}$ is analytical. Therefore the branch cut
structure is not affected except for the appropriate degeneracy at the
respective pole at $-\theta_{i}$. In fact, an even multiplicity can appear in
two types:

  1. Type 1

appears if the branch cut $C_{r}$ reduces to a single point i.e.
$\theta_{2r}=\theta_{2r-1}$ (see the pole $-\theta_{1}=-\theta_{2}$ in Figure
(2), for $p=3$ and $r=1$). This occurs in general if the number of distinct
$\theta_{i}$ less than $\theta_{2r-1}$ is even.

  2. Type 2

appears if two adjacent branch cuts merge so that the joining extreme appears
with multiplicity 2 (see the pole $-\theta_{3}=-\theta_{2}$ in Figure (2)).

Figure 1: Contour limit as $\theta_{2}\to\theta_{1}$ for $p=3$.

Figure 2: Contour limit as $\theta_{2}\to\theta_{3}$ for $p=3$.

###### Theorem 1

Assume that $p$ distinct and increasingly ordered $\theta_{i}$ have
multiplicities $n_{i}$ generating $r$ disjoint segments of branch cuts $C_{k}$
and $l$ particular poles of even multiplicities of type 1. Then the
corresponding pdf and cdf expressions of the random variable
$X=\sum_{i=1}^{p}\frac{1}{2\theta_{i}}\chi^{2}_{n_{i}}(\frac{\gamma_{i}^{2}}{2\theta_{i}})$
are

$f_{\bm{\theta},\bm{\gamma}^{2},\bm{n}}(s)=\frac{f_{\bm{\theta}s,\bm{\gamma}^{2}s,\bm{n}}(1)}{s}=\kappa\left(\frac{-1}{2\pi\textbf{i}}\sum_{i=1}^{r}\int\limits_{\Gamma^{i}}g(t)e^{st}dt+\frac{-1}{2\pi\textbf{i}}\sum_{j=1}^{l}\int\limits_{\Gamma^{j}_{\varepsilon}}g(t)e^{st}dt\right)$

where

$g(t)=\frac{\exp(\sum_{i=1}^{p}\frac{\gamma_{i}^{2}}{4(\theta_{i}+t)})}{\prod_{i=1}^{p}(\theta_{i}+t)^{\frac{n_{i}}{2}}}\quad\kappa=\prod_{i=1}^{p}\theta_{i}^{\frac{n_{i}}{2}}\exp(-\sum_{i=1}^{p}\frac{\gamma^{2}_{i}}{4\theta_{i}})$

and $\Gamma^{i}$ are the non overlapping contours enveloping the branch cut
regions $C_{i}$ determined by $\theta_{i}$ while $\Gamma^{j}_{\varepsilon}$
are the non overlapping contours around the points with even multiplicities
$\theta_{j}$ of type 1. Additionally,

$\displaystyle P(X\leq x)$ $\displaystyle=$
$\displaystyle\frac{e^{x\theta_{0}}}{\theta_{0}}f_{\tilde{\bm{\theta}},\tilde{\bm{\gamma}}^{2}}(x)\prod_{i=1}^{p}\frac{\sqrt{\theta_{i}}}{\sqrt{\theta_{i}+\theta_{0}}}\exp(-\sum_{i=1}^{p}\frac{\gamma^{2}_{i}\theta_{0}}{4\theta_{i}(\theta_{i}+\theta_{0})})$
$\displaystyle=$
$\displaystyle\frac{e^{x\theta_{0}}}{x\theta_{0}}f_{x\tilde{\bm{\theta}},x\tilde{\bm{\gamma}}^{2}}(1)\prod_{i=1}^{p}\frac{\sqrt{\theta_{i}}}{\sqrt{\theta_{i}+\theta_{0}}}\exp(-\sum_{i=1}^{p}\frac{\gamma^{2}_{i}\theta_{0}}{4\theta_{i}(\theta_{i}+\theta_{0})})$

where $f_{\tilde{\bm{\theta}},\tilde{\bm{\gamma}}^{2}}$ is the density
function of the linear combination
$\tilde{X}=\frac{1}{2\theta_{0}}\chi^{2}_{2}(0)+\sum_{i=1}^{p}\frac{1}{2(\theta_{i}+\theta_{0})}\chi^{2}_{n_{i}}(\frac{\gamma_{i}^{2}}{2(\theta_{i}+\theta_{0})}),\forall\theta_{0}>0$.

The proof is presented in Appendix.

Theorem 1 suggests that it is sufficient to integrate along the non-degenerate
contours like “dog-bone” type as $\Gamma^{1}$ or “keyhole” type as
$\Gamma^{2}$ around branch cuts $C_{1}$ and $C_{2}$, see Figure (4).
Integration along such curves, if appropriately parametrized, gives rise to
univariate integrals. For example, the relevant term here
$\frac{-1}{2\pi\textbf{i}}\int\limits_{\Gamma^{i}}g(t)e^{st}dt$ is only the
imaginary part of the integrand so that one could use
$\frac{-1}{2\pi}\int\limits_{\Gamma^{i}}\textbf{Im}(g(t)e^{st})dt$ instead. In
our practical implementation reported in Section 5, we use the Romberg
integration method on spherical contours which seems to work well within the R
package which can be easily utilised to use the multiple precision routines.
If more control or efficiency is required other lower level programming tools
can in principle be used.

In order to obtain the contours $\Gamma^{i}$ and $\Gamma^{j}_{\varepsilon}$ of
Theorem above, for a given set of $p$ distinct $\theta_{i}$ with
multiplicities $n_{i}$, one could perform the following two steps:

  1. 1.

Consider initially only those $\theta_{i}$ which have odd multiplicities and
pair them to determine the corresponding branch cuts $C_{i}$ as in
$\Gamma^{i}$’s of Theorem.

  2. 2.

Consider next only those $\theta_{i}$ with even multiplicities(type 1) which
are outside $\Gamma^{i}$ regions and calculate their residues either
numerically by integrating along contours $\Gamma^{j}_{\varepsilon}$ or
analytically if $\gamma_{i}=0$.

In fact the integrating contours reduce to line segments for the practically
important case of all terms being non-central chi squares with single degrees
of freedom, namely, the case of $p$ distinct $\theta_{i}$, each having
multiplicity $n_{i}=1$ and $\gamma_{i}=0$.

###### Theorem 2

For distinct $\theta_{i}>0$ and $\gamma_{i}=0$, the corresponding pdf and cdf
expressions for the random variable
$X=\sum_{i=1}^{p}\frac{1}{2\theta_{i}}\chi^{2}_{1}(0)$ are:

$\displaystyle f_{\bm{\theta},\bm{0},\bm{1}}(s)$ $\displaystyle=$
$\displaystyle\frac{\prod_{i=1}^{p}\sqrt{\theta_{i}}}{\pi}\left(\int_{\theta_{1}}^{\theta_{2}}\frac{e^{-st}}{\sqrt{-\prod_{i=1}^{p}(\theta_{i}-t)}}dt-\int_{\theta_{3}}^{\theta_{4}}\frac{e^{-st}}{\sqrt{-\prod_{i=1}^{p}(\theta_{i}-t)}}dt+\int_{\theta_{5}}^{\theta_{6}}\frac{e^{-st}}{\sqrt{-\prod_{i=1}^{p}(\theta_{i}-t)}}dt\cdots\right)$
(9) $\displaystyle=$
$\displaystyle\frac{\prod_{i=1}^{p}\sqrt{\theta_{i}}}{\pi}\sum_{r=1}^{[\frac{p+1}{2}]}(-1)^{r+1}\int_{\theta_{2r-1}}^{\theta_{2r}}\frac{e^{-st}}{\sqrt{-\prod_{i=1}^{p}(\theta_{i}-t)}}dt$

where

$\theta_{2[\frac{p+1}{2}]}=\left\\{\begin{array}[]{cc}\theta_{p}&\text{p
even}\\\ \infty&\text{p odd}\end{array}\right.$

and

$P(X\leq
x)=f_{\tilde{\bm{\theta}},\bm{0},\bm{1}}(x)\frac{e^{x\theta_{0}}\prod\limits_{i=1}^{p}\sqrt{\theta_{i}}}{\theta_{0}\prod_{i=1}^{p}\sqrt{\theta_{i}+\theta_{0}}}=f_{\tilde{\bm{\theta}}x,\bm{0},\bm{1}}(1)\frac{e^{x\theta_{0}}\prod\limits_{i=1}^{p}\sqrt{\theta_{i}}}{x\theta_{0}\prod_{i=1}^{p}\sqrt{\theta_{i}+\theta_{0}}}\quad\forall\theta_{0}>0$

where $f_{\tilde{\bm{\theta}},\bm{0},\bm{1}}$ is the pdf of
$\tilde{X}=\frac{1}{2\theta_{0}}\chi^{2}_{2}(0)+\sum_{i=1}^{p}\frac{1}{2(\theta_{i}+\theta_{0})}\chi^{2}_{n_{i}}(0)$.

The proof is presented in Appendix.

Figure 3: Contours for non central cases

Figure 4: Contour limit as $\nu,\varepsilon\to 0$

In fact these theorems indicate that the regions where the relevant quantities
need to be integrated are only those around the branch cuts appearing from
$\theta$’s with odd multiplicities (or chi-squared with odd degrees of
freedom). In the central cases as in Theorem 2, these regions can be also
reduced to simple lines. Contour integration is also valid along any lines
containing all these branch cut regions and in fact any such choice will lead
to a different integrating parametrization.

###  2.1 Computation remarks for the central cases $\gamma_{i}=0$

The univariate integrating terms of Theorem 2 are easily seen to be finite
while can have alternative representation:

$\displaystyle\int\limits_{\theta_{2r-1}}^{\theta_{2r}}\frac{e^{-st}}{\sqrt{-\prod_{i=1}^{p}(\theta_{i}-t)}}dt$
$\displaystyle=$ $\displaystyle
e^{-s\theta_{2r-1}}\int_{0}^{1}\frac{e^{-su\alpha_{r}}du}{\sqrt{u(1-u)\prod\limits_{2r-1\neq
i\neq 2r}^{p}(\theta_{i}-\theta_{2r-1}-u\alpha_{r})}}$ $\displaystyle=$
$\displaystyle
2e^{-s\theta_{2r-1}}\int_{0}^{\pi/2}\frac{\text{exp}(-s\alpha_{r}\sin^{2}u)du}{\sqrt{\prod\limits_{2r-1\neq
i\neq 2r}^{p}(\theta_{i}-\theta_{2r-1}\cos^{2}u-\theta_{2r}\sin^{2}u)}}$

where $\alpha_{r}=\theta_{2r}-\theta_{2r-1}$ and
$t=\theta_{2r-1}\cos^{2}u+\theta_{2r}\sin^{2}u$. In the type one multiplicity
case i.e. $\alpha_{r}=0$, one can progress by taking the limit on (2.1) as
$\alpha_{r}\to 0$:

$\int\limits_{\theta_{2r-1}}^{\theta_{2r}}\frac{e^{-st}}{\sqrt{-\prod_{i=1}^{p}(\theta_{i}-t)}}dt=\xrightarrow[\alpha_{r}\to
0]{}\frac{\pi e^{-s\theta_{2r-1}}}{\prod\limits_{2r-1\neq i\neq
2r}^{p}\sqrt{\theta_{i}-\theta_{2r}}}$ (11)

and therefore there is no need to numerically integrate this term. This is not
surprising as in this case the pole is of order 1 and therefore the residue
theorem implies the result above.

Similarly, when $p$ is odd, the last term corresponds to the branch cut
$(-\infty,-\theta_{p})$

$\int\limits^{\infty}_{\theta_{p}}\frac{e^{-st}}{\sqrt{-\prod_{i=1}^{p}(\theta_{i}-t)}}dt=\int\limits^{\infty}_{0}\frac{\text{exp}(-s\theta_{p}-sv)}{\sqrt{\prod_{i=1}^{p-1}(\theta_{i}-\theta_{p}-v)}\sqrt{v}}dt=2e^{-s\theta_{p}}\int\limits_{0}^{+\infty}\frac{e^{-su^{2}}du}{\sqrt{\prod_{i=1}^{p-1}(\theta_{p}-\theta_{i}+u^{2})}}$
(12)

where and $t=\theta_{p}+v$, $v=u^{2}$. These terms together with those of
(2.1), can also be efficiently integrated numerically as these functions are
either decaying exponentially to zero or having a “U” shape with only one
critical point. See Figure 5 for the log-scale behaviour of these elementary
functions.

Figure 5: Function $g(t)e^{-st}$ in log-scale for the central case with
$\bm{\theta}=c(1,2,3,4,5)$ and $s=1$.

###  2.2 Adding degrees of freedom (multiplicities of $\theta_{1}$) higher
than 2

A recursive application of the derivatives to the normalising constants as in
Kume and Wood (2007) or as generalised in (7) could in principle generate the
required expressions. The differentiation needs to be applied to each of the
elementary terms of Theorem 2 and a simple differentiation has the effect
adding the multiplicity by 2. For example, the differentiation as below for
each of the terms of Theorem 2 leads to the cases of multiplicity 3:

$\frac{\partial}{\partial\theta_{j}}\int\limits_{C_{r}}\frac{e^{t}}{\sqrt{-\prod_{i=1}^{p}(\theta_{i}+t)}}dt=\left\\{\begin{array}[]{lr}-\frac{1}{2}\int\limits_{-\theta_{2r}}^{-\theta_{2r-1}}\frac{1}{(\theta_{j}+t)\sqrt{-\prod_{i=1}^{p}(\theta_{i}+t)}}e^{st}dt&\theta_{2r}\neq\theta_{j}\neq\theta_{2r-1}\\\
\frac{1}{\theta_{2r}-\theta_{2r-1}}\int\limits_{-\theta_{2r}}^{-\theta_{2r-1}}\frac{\theta_{2r-1}+t}{\sqrt{-\prod_{i=1}^{p}(\theta_{i}+t)}}e^{st}\left(-1+\frac{1}{2}\sum\limits_{2r-1\neq
i\neq 2r}\frac{1}{\theta_{i}+t}\right)dt&j=\theta_{2r-1}\\\
\frac{1}{\theta_{2r}-\theta_{2r-1}}\int\limits_{-\theta_{2r}}^{-\theta_{2r-1}}\frac{\theta_{2r}+t}{\sqrt{-\prod_{i=1}^{p}(\theta_{i}+t)}}e^{st}\left(-1+\frac{1}{2}\sum\limits_{2r-1\neq
i\neq 2r}\frac{1}{\theta_{i}+t}\right)dt&j=\theta_{2r}\end{array}\right.$

##  3 Cases of practical importance

###  3.1 Bingham distribution

For the Bingham distributions with ordered and distinct $\theta_{i}$ we have
that

$\displaystyle\mathcal{C}(\bm{\theta})$ $\displaystyle=$
$\displaystyle\int\limits_{\mathcal{S}^{p-1}}\text{exp}(\sum_{i=1}^{p}-\theta_{i}x_{i}^{2})d_{\mathcal{S}^{p-1}}(\bm{x})=2\pi^{\frac{\sum_{i=1}^{p}n_{i}}{2}-1}\sum_{r=1}^{[\frac{p+1}{2}]}(-1)^{r+1}\int_{\theta_{2r-1}}^{\theta_{2r}}\frac{e^{-t}}{\sqrt{-\prod_{i=1}^{p}(\theta_{i}-t)}}dt.$
(13)

If $p=3$, the Bingham on the ordinary sphere has normalising constant

$\frac{\mathcal{C}(\bm{\theta})}{4\sqrt{\pi}}=e^{-\theta_{1}}\int_{0}^{\pi/2}\frac{\text{exp}(-(\theta_{2}-\theta_{1})\sin^{2}u)du}{\sqrt{\theta_{3}-\theta_{1}\cos^{2}u-\theta_{2}\sin^{2}u}}du-e^{-\theta_{3}}\int\limits_{0}^{+\infty}\frac{e^{-u^{2}}du}{\sqrt{\prod_{i=1}^{2}(\theta_{3}-\theta_{i}+u^{2})}}$

For $p=4$, there are only two finite branch cuts, and from (2.1) we can
generate this expression as

$\displaystyle\frac{\mathcal{C}(\bm{\theta})}{4\pi}$ $\displaystyle=$
$\displaystyle
e^{-\theta_{1}}\int_{0}^{\pi/2}\frac{\text{exp}\left(-(\theta_{2}-\theta_{1})\sin^{2}u\right)du}{\sqrt{\prod\limits_{i=3}^{4}(\theta_{i}-\theta_{1}\cos^{2}u-\theta_{2}\sin^{2}u)}}du-e^{-\theta_{3}}\int_{0}^{\pi/2}\frac{\text{exp}\left(-(\theta_{4}-\theta_{3})\sin^{2}u\right)du}{\sqrt{\prod\limits_{i=1}^{2}(\theta_{i}-\theta_{3}\cos^{2}u-\theta_{4}\sin^{2}u)}}du$
(14) $\displaystyle=$
$\displaystyle\frac{e^{-\theta_{1}}}{2}\int_{0}^{1}\frac{e^{-u\alpha_{1}}du}{\sqrt{u(1-u)(\theta_{3}-\theta_{1}-u\alpha_{1})(\theta_{4}-\theta_{1}-u\alpha_{1})}}$
$\displaystyle-\frac{e^{-\theta_{3}}}{2}\int_{0}^{1}\frac{e^{-u\alpha_{2}}du}{\sqrt{u(1-u)(\theta_{1}-\theta_{3}-u\alpha_{2})(\theta_{2}-\theta_{3}-u\alpha_{2})}}$

where $\alpha_{1}=\theta_{2}-\theta_{1}$ and
$\alpha_{2}=\theta_{4}-\theta_{3}$. In the light of property (7),
differentiating for various $\theta_{i}$ here leads to other representations
corresponding to additional multiplicities.

###  3.2 Fisher distribution of rotations in $\mathbb{R}^{3}$ and its use in
3-D shape inference

The expression of Bingham distributions for $p=4$, can be adopted as in Wood
(1993) so that for some diagonal matrix

$\Phi=\left(\begin{array}[]{ccc}\phi_{1}&0&0\\\ 0&\phi_{2}&0\\\
0&0&\phi_{3}\end{array}\right)\quad\text{and}\quad\bm{\theta}=-2(\phi_{1},\phi_{2},\phi_{3},\phi_{1}+\phi_{2}+\phi_{3})$

the normalizing constant of the Fisher distribution in $SO(3)$ is

$\int_{SO(3)}e^{\text{tr}(\Phi\mathbf{R})}d\mathbf{R}=\frac{1}{2}e^{\phi_{1}+\phi_{2}+\phi_{3}}\mathcal{C}(\bm{\theta})\quad$

where the integration is taking place along the Haar uniform measure
$d\mathbf{R}$ on the group of rotations $SO(3)$ and the evaluation of
$\mathcal{C}(\bm{\theta})$ is as in (3.1). This formula is similar to the one
shown in Wood (1993).

In Dryden et al. (2021), the likelihood function for a very general set of
regression shape models is based on the efficient evaluation of this
normalising constant. In addition they address the bias correction of the well
known procrustes mean shape estimators for the 3-d shape models based on the
principles of the maximum likelihood estimation. For establishing the
resulting estimators the authors suggest using the gradient of
$\log\int_{SO(3)}e^{\text{tr}(\Phi R)}dR$ as a bias correcting term. In
particular they use this result

$\nabla_{\Phi}\log\int_{SO(3)}e^{\text{tr}(\Phi\mathbf{R})}d\mathbf{R}=\left(\begin{matrix}1-\frac{2}{\mathcal{C}(\bm{\theta})}\frac{\partial\mathcal{C}(\bm{\theta})}{\partial\theta_{1}}-\frac{2}{\mathcal{C}(\bm{\theta})}\frac{\partial\mathcal{C}(\bm{\theta})}{\partial\theta_{4}}\\\
1-\frac{2}{\mathcal{C}(\bm{\theta})}\frac{\partial\mathcal{C}(\bm{\theta)}}{\partial\theta_{2}}-\frac{2}{\mathcal{C}(\bm{\theta})}\frac{\partial\mathcal{C}(\bm{\theta})}{\partial\theta_{4}}\\\
1-\frac{2}{\mathcal{C}(\bm{\theta})}\frac{\partial\mathcal{C}(\bm{\theta)}}{\partial\theta_{3}}-\frac{2}{\mathcal{C}(\bm{\theta})}\frac{\partial\mathcal{C}(\bm{\theta})}{\partial\theta_{4}}\end{matrix})\right)$

where the partial derivative expressions are easily obtained from the integral
representations as in (13). In fact Dryden et al. (2021) show that that for
practical purposes, these components are accurately and very easily calculated
using the saddlepoint approximation.

###  3.3 Kent’s formula for Complex Bingham

If even multiplicities are introduced, the limit of the expression (13) will
suffice. For example, if all $\theta_{i}$ as in Lemma 1 are of even
multiplicities then the corresponding function $g(t)$ will not have any branch
cuts but just simple poles therefore a simple application of the residue
theorem generates the wellknown result of Kent for the complex Bingham
distributions. More specifically, if $p=2k$ while $\theta_{i}$’s are
coalescing in pairs i.e. $\theta_{2r}\to\theta_{2r-1}$, and following similar
arguments to those in (11), each term of (13) becomes

$\frac{e^{-\theta_{2r}}}{\prod\limits_{2r-1\neq i\neq
2r}^{p}(\theta_{i}-\theta_{2r})^{1/2}}\int_{0}^{1}\frac{du}{\sqrt{u(1-u)}}=\frac{e^{-\theta_{2r}}}{\prod\limits_{2r-1\neq
i\neq 2r}^{k}(\theta_{i}-\theta_{2r})}\pi$

and therefore for $k$ distinct pairwise $\theta$’s.

$\mathcal{C}(\bm{\theta})=2\pi^{k}\sum_{r=1}^{k}(-1)^{r+1}\frac{e^{-\theta_{r}}}{\prod\limits_{i\neq
r}^{p}(\theta_{i}-\theta_{r})}=2\pi^{k}\sum_{r=1}^{k}\frac{e^{-\theta_{r}}}{\prod\limits_{i\neq
r}^{p}(\theta_{r}-\theta_{i})}$ (16)

which is consistent with the result of Kent for the complex Bingham
distribution. Note that in the light of (6) for $n_{i}=2$ and $\gamma_{i}=0$,
Kent’s normalising constant above can also be used for close form expressions
of the pdf and cdf of a positive linear combination of the exponential
distributing terms see e.g. Kent (1994b).

###  3.4 Kent distributions

In these distributions, we have only one $\gamma_{r}\neq 0$ and in particular
for the $p=3$ case, $\gamma_{2}\neq 0$ and therefore, the integration along
the unbounded branch cut is reduced as in Theorem 2.

An additional feature for the Kent distribution is that the middle pole
$\theta_{2}$ is equally distant from the other two so that
$\bm{\theta}=(\theta_{1},\theta_{1}+\alpha,\theta_{1}+2\alpha)$. Since the
non-central term is not on the unbounded branch $(-\infty,-\theta_{3})$, one
can reduce that integration along the real line throughout this interval and
use an additional circular contour parametrisation
$t=-\theta_{1}-\frac{\alpha}{2}+re^{\textbf{i}u}$ for the finite branch, and
so for the Kent distribution:

$\displaystyle\frac{\mathcal{C}(\bm{\theta},\gamma_{2})}{2e^{-\theta_{1}}\sqrt{\pi}}$
$\displaystyle=$
$\displaystyle\frac{-re^{-\frac{\alpha}{2}}}{2\pi}\int\limits_{0}^{2\pi}\textbf{Im}\left(\frac{\text{exp}\left(\frac{\gamma_{2}^{2}}{4(\frac{\alpha}{2}+re^{\textbf{i}u})}+re^{\textbf{i}u}\right)}{\sqrt{(-\frac{\alpha}{2}+re^{\textbf{i}u})(\frac{\alpha}{2}+re^{\textbf{i}u})(3\frac{\alpha}{2}+re^{\textbf{i}u})}}\right)du$
$\displaystyle-$ $\displaystyle
2e^{-2\alpha}\int\limits_{0}^{+\infty}\frac{\text{exp}\left(-u^{2}-\frac{\gamma_{2}^{2}}{4(\alpha+u^{2})}\right)du}{\sqrt{(\alpha+u^{2})(2\alpha+u^{2})}}$

where the radius $\frac{\alpha}{2}<r<\alpha$ is chosen so that the circle
contains only the first branch cut. Note that for implementing this
distribution in directional statistics, $\theta_{1}$ can be allowed to be $0$
and therefore the normalizing constant of this distribution is in fact defined
in terms of only two parameters $\alpha$ and $\gamma_{2}$. These parameters
correspond to $\beta$ and $\kappa$ in Kent (1982).

##  4 Difference of two positive linear combinations

If in general some of the $\lambda_{i}$ in (1) are negative then the
characterisation of the problem in terms of normal distribution components is
as the difference between two quadratic forms. The contour integration however
can be defined along the vertical lines such that the strictly negative poles
are on the left of the contours $\textbf{i}\mathbb{R}+t_{0}$ while $t_{0}<0$
as in (2). More specifically, if $Z=X-Y$ where both $X$ and $Y$ have the same
general expression as in (1) for some positive parameters $\bm{\theta}$,
$\bm{\gamma}$, $\bm{\theta}^{\prime}$ and $\bm{\gamma}^{\prime}$, with
respective densities

$f_{\bm{\theta},\bm{\gamma}^{2},\bm{n}}(s)=\kappa\frac{1}{2\pi\textbf{i}}\int\limits_{\textbf{i}\mathbb{R}+t_{0}}g(t)e^{st}dt\quad
f_{\bm{\theta}^{\prime},\bm{\gamma}^{\prime
2},\bm{n}^{\prime}}(s)=\kappa^{\prime}\frac{1}{2\pi\textbf{i}}\int\limits_{\textbf{i}\mathbb{R}+v_{0}}g^{\prime}(v)e^{sv}dv\quad
s>0$ (17)

where
$\kappa=\prod_{i=1}^{p}\theta_{i}^{\frac{n_{i}}{2}}\text{exp}(-\sum_{i=1}^{p}\frac{\gamma^{2}_{i}}{4\theta_{i}})$,
$\kappa^{\prime}=\prod_{i=1}^{p^{\prime}}\theta_{i}^{\prime\frac{n^{\prime}_{i}}{2}}\text{exp}(-\sum_{i=1}^{p}\frac{\gamma^{\prime
2}_{i}}{4\theta^{\prime}_{i}})$ and

$g(t)=\frac{\text{exp}(\sum_{i=1}^{p}\frac{\gamma_{i}^{2}}{4(\theta_{i}+t)})}{\prod_{i=1}^{p}(\theta_{i}+t)^{\frac{n_{i}}{2}}}\quad
g^{\prime}(v)=\frac{\text{exp}(\sum_{i=1}^{p}\frac{\gamma^{\prime
2}_{i}}{4(\theta^{\prime}_{i}+v)})}{\prod_{i=1}^{p}(\theta^{\prime}_{i}+v)^{\frac{n_{i}^{\prime}}{2}}}$

the following result holds:

###### Theorem 3

If
$X=\sum_{i=1}^{p}\frac{1}{2(\theta_{i})}\chi^{2}_{n_{i}}(\frac{\gamma_{i}^{2}}{2(\theta_{i})})$
and
$Y=\sum_{i=1}^{p^{\prime}}\frac{1}{2(\theta^{\prime}_{i})}\chi^{2}_{n^{\prime}_{i}}(\frac{\gamma_{i}^{\prime
2}}{2(\theta^{\prime}_{i})})$ then the pdf of $Z=X-Y$ at $z\geq 0$ is

$f_{Z}(z)=\prod_{i=1}^{p}\theta_{i}^{\frac{n_{i}}{2}}\exp(-\sum_{i=1}^{p}\frac{\gamma^{2}_{i}}{4\theta_{i}})\prod_{i=1}^{p}\theta_{i}^{\prime\frac{n^{\prime}_{i}}{2}}\exp(-\sum_{i=1}^{p}\frac{\gamma^{\prime
2}_{i}}{4\theta^{\prime}_{i}})\frac{-1}{2\pi\textbf{i}}\int\limits_{\Gamma}g(t)g^{\prime}(-t)e^{zt}dt$

where the integrating contour $\Gamma$ can be simplified as in Theorems 1 and
2 along branch cuts involving only the values of $\bm{\theta}$. Additionally

$P(X-Y>z\geq
0)=\frac{e^{\theta_{0}z}}{\theta_{0}}f_{\tilde{Z}}(z)\prod_{i=1}^{p}\frac{\sqrt{\theta_{i}}}{\sqrt{\theta_{i}+\theta_{0}}}\prod_{i=1}^{p^{\prime}}\frac{\sqrt{\theta^{\prime}_{i}}}{\sqrt{\theta^{\prime}_{i}-\theta_{0}}}\exp(\sum_{i=1}^{p}\frac{\gamma^{2}_{i}}{4\theta_{i}}-\frac{\gamma^{2}_{i}}{4(\theta_{i}+\theta_{0})}+\sum_{i=1}^{p^{\prime}}\frac{\gamma^{\prime
2}_{i}}{4\theta^{\prime}_{i}}-\frac{\gamma^{\prime
2}_{i}}{4(\theta^{\prime}_{i}-\theta_{0})})$

for any arbitrary $0<\theta_{0}<min(\bm{\theta}^{\prime})$ and

$\tilde{Z}=\frac{1}{2\theta_{0}}\chi^{2}_{2}(0)+\sum_{i=1}^{p}\frac{1}{2(\theta_{i}+\theta_{0})}\chi^{2}_{n_{i}}(\frac{\gamma_{i}^{2}}{2(\theta_{i}+\theta_{0})})-\sum_{i=1}^{p^{\prime}}\frac{1}{2(\theta^{\prime}_{i}-\theta_{0})}\chi^{2}_{n^{\prime}_{i}}(\frac{\gamma_{i}^{\prime
2}}{2(\theta^{\prime}_{i}-\theta_{0})}).$

The corresponding expressions for $z<0$ can be obtained similarly for
$f_{Y-X}(-z)$ by switching the roles of $X$ and $Y$ as in $-Z=Y-X$.

As a special case one could think of the exponential distributions which
correspond modulo some scaling to a central chi-square of degree. These cases
have practical implications due to their connections with the Erlang
distributions and Phase Type distributions in the probability theory. A recent
paper for generating the distribution for the positive entries as in Neumüller
et al. (2022) and also the proof for the extended case of real coefficients
can be found in Mathai (1983).

In our construction however, these expressions are simply the residues of
ordinary poles. Similar to the case of complex Bingham distributions which are
also related to the exponentially distributed components, the respective
branch cuts reduce to simple poles and the residues generate closed form
expressions:

###### Corollary 4

If the chi-square components in Theorem 3 are central chi-squares r.v’s with
degrees of freedom 2, namely $Z=X-Y$ represents a linear combination of
exponential random variables, then the corresponding density and distribution
functions at $z\geq 0$ are evaluated as

$f_{Z}(z)=\prod_{i=1}^{p}\theta_{i}\prod_{j=1}^{p^{\prime}}\theta_{j}^{\prime}\sum_{i=1}^{p}(-1)^{i}\frac{e^{-\theta_{i}z}}{\prod_{j\neq
i}^{p}(\theta_{j}-\theta_{i})\prod_{k=1}^{p^{\prime}}(\theta^{\prime}_{k}-\theta_{i})}$

and

$P(Z>z\geq
0)=\prod_{i=1}^{p}\theta_{i}\prod_{j=1}^{p^{\prime}}\theta_{j}^{\prime}\sum_{i=1}^{p}(-1)^{i}\frac{e^{-\theta_{i}z}}{\theta_{i}\prod_{j\neq
i}^{p}(\theta_{j}-\theta_{i})\prod_{k=1}^{p^{\prime}}(\theta^{\prime}_{k}-\theta_{i})}.$

For the cases when $z<0$ the roles of $\bm{\theta}$ and $\bm{\theta^{\prime}}$
are reversed.

##  5 Numerical evidence

We illustrate our method by replicating the numerical figures reported in
Imhof (1961). We focus initially on the first two central cases reported
there, where the coefficients lambda of $\eqref{eqn:gen:rep}$ are defined as
$\mathbf{\lambda}_{i}=\frac{1}{2\bm{\theta}_{i}}$ which in our parametrization
they correspond to:

$\bm{\theta}_{1}=(0.8333333,1.6666667,5.0000000)\quad\bm{\theta}_{2}=(0.8333333,0.8333333,1.6666667,1.6666667,5,5)$

Note in particular that for $\bm{\theta}_{2}$ above, the three distinct values
have multiplicity $2$ i.e. simple poles and therefore as we have noted in the
corresponding inversion is immediately available in closed form for both pdf
and cdf using (16) and (6). In Imhof (1961) only the cdf values for specific
entries are shown while we can see that the whole function can be easily
obtained using our simple contour integrating terms. The comparison of both
exact and SPA methods for density, pdf. and cdf. together with the
corresponding relative error are shown in Figures (7) and (7). It can be
easily seen that the density SPA used for the cdf calculations as stated in
our Theorem 1, performs very well especially in the extreme values of the
variable $s$.

Figure 6: Comparisons for $\bm{\theta}_{1}$

Figure 7: Comparisons for $\bm{\theta}_{2}$

As illustrative evidence for the difference of two positive linear
combinations, we can mention cases reported in Imhof (1961) where in
particular, the three example entries on our Table 1 here relate to the cases
$\frac{Q_{3}}{3}-\frac{2Q_{4}}{3}$ and $\frac{Q_{5}}{2}-\frac{Q_{6}}{2}$ and
$\frac{1}{6}Q_{3}+\frac{2}{6}Q_{6}-\frac{1}{6}Q_{5}-\frac{2}{6}Q_{4}$ there.

The approximation method reported in Imhof (1961) is illustrated against some
exact figures which are also reported in our Table 1. These exact values which
the Imhof’s method is comparing to are in full agreement with our method of
bounded contour integration as seen for various values $s$. Our method can
easily produce the pdf and cdf of these three cases as shown on Figure 8.

Table 1: Some numerical checks with those of Imhof (1961)

Model |  Our parametrisation $\bm{\theta}^{\mathbf{n}}_{\mathbf{\gamma}}$ | s | Imhof (1961) | Our Method
---|---|---|---|---
$\frac{1}{3}Q_{3}-\frac{2}{3}Q_{4}$ | $0.4^{2}_{0}+0.2^{4}_{0}+(\frac{2}{30})^{6}_{0}-0.35^{1}_{6}-0.15^{1}_{2}$ | -2 | 0.9102 | 0.9102254
|  | 0 | 0.4061 | 0.4061061
|  | 2.5 | 0.097 | 0.0097598
$\frac{Q_{5}}{2}-\frac{Q_{6}}{2}$ | $3.5^{6}_{6}+0.15^{2}_{2}-0.35^{1}_{6}-0.15^{1}_{2}$ | -2 | 0.9218 | 0.921792
|  | 2 | 0.4779 | 0.4778933
|  | 7 | 0.0396 | 0.0396319
$\frac{1}{6}(Q_{3}-Q_{5})+\frac{2}{6}(Q_{6}-Q_{4})$ | $0.1_{2}^{7}+0.05^{4}_{0}+(\frac{0.1}{6})^{2}_{0}+(\frac{1.4}{6})^{1}_{6}$ | -3 | 0.9861 | 0.9861469
| $-0.2_{0}^{2}-0.1^{4}_{0}-(\frac{0.2}{6})^{6}_{0}-(\frac{0.7}{6})^{6}_{6}-0.05^{2}_{2}$ | 0 | 0.5170 | 0.5170232
|  | 4 | 0.0152 | 0.0152041

Figure 8:  A plot of pdf and cdf for$\frac{1}{3}Q_{3}-\frac{2}{3}Q_{4}$ in
black; $\frac{Q_{5}}{2}-\frac{Q_{6}}{2}$ in blue and
$\frac{1}{6}(Q_{3}-Q_{5})+\frac{2}{6}(Q_{6}-Q_{4})$ in red.

##  6 Concluding remarks

In this paper we explore alternative expressions for the integral
representation of this important family of distributions. It is clear that the
branch cut approach covered here has simplicity in expression leading to an
explicit connection between their pdf and df. Therefore one can apply the
standard saddle point approximation method in evaluating tail probabilities.
The resulting expressions in terms or branch cuts are also easy to implement
using the standard numerical routines available in many computer packages. By
identifying these cuts enables the user to vary the integrating contours and
not just follow standard inversion methods or its improvements like Talbot
methods. The saddle point approximation is shown to work well for a range of
values, with a relative improvement in the tails of the distributions but such
approximations are generally not performing well in the case of negative
entries.

## References

  * Ablowitz and Fokas (2003) Ablowitz, M., Fokas, A., 2003. Introduction and applications of complex variables. Cambridge University Press.
  * Barndorff-Nielsen (1990) Barndorff-Nielsen, O. E., 1990. Approximate interval probabilities. Journal of the Royal Statistical Society. Series B 52, 485–496.
  * Butler (2007) Butler, R. W., 2007. Saddlepoint approximations with applications. Vol. 22 of Camb. Ser. Stat. Probab. Math. Cambridge: Cambridge University Press.
  * Chen and Tanaka (2021) Chen, Y., Tanaka, K., 2021. Maximum likelihood estimation of the fisher–bingham distribution via efficient calculation of its normalizing constant. Statistics and Computing 31.
  * Daniels (1954) Daniels, H. E., 1954. Saddlepoint approximations in statistics. Annals of Mathematical Statistics 25, 631–650.
  * Davies (1973) Davies, R. B., 1973. Numerical inversion of a characteristic function. Biometrika, 415–417.
  * Davies (1980) Davies, R. B., 1980. The distribution of a linear combination of x2 random variables. Applied Statistics 29 (3), 323–333.
  * Dryden et al. (2021) Dryden, I. L., Kume, A., Paine, P. J., Wood, A. T. A., 2021. Regression modeling for size-and-shape data based on a Gaussian model for landmarks. Journal of the American Statistical Association 116 (534), 1011–1022.
  * Duchesne and De Micheaux (2010) Duchesne, P., De Micheaux, P. L., 2010. Computing the distribution of quadratic forms: Further comparisons between the liu–tang–zhang approximation and exact methods. Computational Statistics & Data Analysis 54 (4), 858–862.
  * Imhof (1961) Imhof, J., 1961. Computing the distribution of quadratic forms in normal variables. Biometrika.
  * Kent (1982) Kent, J. T., 1982. The Fisher–Bingham distribution on the sphere. Journal of the Royal Statistical Society: Series B 44, 71–80.
  * Kent (1994a) Kent, J. T., 1994a. The complex bingham distribution and shape analysis. Journal of the Royal Statistical Society. Series B (Methodological), 285–299.
URL http://www.jstor.org/stable/10.2307/2345900

  * Kent (1994b) Kent, J. T., 1994b. The complex bingham distribution and shape analysis. J. Roy. Stat. Soc., Series B 56, 285–299.
  * Koyama and Takemura (2016) Koyama, T., Takemura, A., 2016. Holonomic gradient method for distribution function of a weighted sum of noncentral chi-square random variables. Computational Statistics 31 (4), 1645–1659.
  * Kume and Sei (2018) Kume, A., Sei, T., Jul. 2018. On the exact maximum likelihood inference of fisher–bingham distributions using an adjusted holonomic gradient method. Statistics and Computing 28 (4), 835.
URL http://dx.doi.org/10.1007/s11222-017-9765-3

  * Kume and Wood (2005) Kume, A., Wood, A. T. A., 2005. Saddlepoint approximations for the Bingham and Fisher-Bingham normalising constants. Biometrika 92, 465–476.
  * Kume and Wood (2007) Kume, A., Wood, A. T. A., 2007. On the derivatives of the normalising constant of the Bingham distribution. Statistics and Probability Letters 77, 832–837.
  * Kuonen (1999) Kuonen, D., 1999. Miscellanea. saddlepoint approximations for distributions of quadratic forms in normal variables. Biometrika 86 (4), 929–935.
  * Kuonen (2005) Kuonen, D., 2005. Studentized bootstrap confidence intervals based on m-estimates. Journal of Applied Statistics 32 (5), 443–460.
  * Liang and Self (1996) Liang, K.-Y., Self, S. G., 1996. On the asymptotic behaviour of the pseudolikelihood ratio test statistic. Journal of the Royal Statistical Society. Series B 58, 785–796.
  * Liu et al. (2009) Liu, H., Tang, Y., Zhang, H. H., 2009. A new chi-square approximation to the distribution of non-negative definite quadratic forms in non-central normal variables. Computational Statistics & Data Analysis 53 (4), 853–856.
  * Mardia and Jupp (2000) Mardia, K. V., Jupp, P. E., 2000. Directional Statistics. Wiley Series in Probability and Statistics. John Wiley & Sons Ltd., Chichester.
  * Mathai (1983) Mathai, A. M., 1983. On linear combinations of independent exponential variables. Communications in Statistics - Theory and Methods 12, 625–632.
  * Nakayama et al. (2011) Nakayama, H., Nishiyama, K., Noro, M., Ohara, K., Sei, T., Takayama, N., Takemura, A., 2011. Holonomic gradient descent and its application to the Fisher-Bingham integral. Advances in Applied Mathematics 47, 639–658.
  * Neumüller et al. (2022) Neumüller, C., Robert, J., Heuberger, A., 9 2022. Convolving hyper-erlang with hyper-exponential distributions using linear algebra. In: Proceedings of the 2022 5th International Conference on Mathematics and Statistics. ICoMS ’22. Association for Computing Machinery, New York, NY, USA, p. 1–6.
URL https://doi.org/10.1145/3545839.3545840

  * Sei and Kume (2013) Sei, T., Kume, A., 2013. Calculating the normalising constant of the bingham distribution on the sphere using the holonomic gradient method. Statistics and Computing, 1–12.
URL http://link.springer.com/article/10.1007/s11222-013-9434-0

  * Wood (1993) Wood, A. T., 1993. Estimation of the concentration parameters of the fisher matrix distribution on 50 3) and the bingham distribution on sq, q 2. Australian Journal of Statistics 35 (1), 69–79.
URL
http://onlinelibrary.wiley.com/doi/10.1111/j.1467-842X.1993.tb01313.x/abstract

##  Appendix A Proof of Theorem 1

In order to give now the idea of the proof, it is sufficient to consider as
above the case $p=3$ while the generalization to larger $p$ easily follows. In
order to see the derivation of the branch cuts first, let us assume for
simplicity that $p=3$ and for a point $t$ in the complex plane just above the
real line segment $[-\theta_{2},-\theta_{1}]$ we have that
$arg(\sqrt{(\theta_{1}+t)}=r_{1}e^{\text{{i}}\alpha_{1}})=\alpha_{1}$ where
$\alpha_{1}\approx\pi/2$ but $\alpha_{2}=arg(\sqrt{(\theta_{2}+t)}\approx 0$
and similarly $\alpha_{3}\approx 0$ and therefore
$\alpha=\alpha_{1}+\alpha_{2}+\alpha_{3}\approx\pi/2$. However for a point
just below $[-\theta_{2},-\theta_{1}]$ the same reasoning leads to
$\alpha=\alpha_{1}+\alpha_{2}+\alpha_{3}\approx-\pi/2$. So
$\alpha\approx\pm\pi/2$ while taking different signs on either sides of
$[-\theta_{2},-\theta_{1}]$ and hence the function $g(t)$ is not analytic in
$[-\theta_{2},-\theta_{1}]$ implying a branch cut. However, if we are to
follow similar reasoning as above for points $t$ on either side of
$[-\theta_{3},-\theta_{2}]$ we then have
$\alpha_{1}\approx\alpha_{2}\approx\pi/2$ and $\alpha_{3}\approx 0$, i.e.
$\alpha\approx\pi$ for the points just above and $\alpha\approx-\pi$ for the
points just below this segment. Since the sign choice for
$\alpha\approx\pm\pi$ does not mater the function $g(t)$ remains analytic in
$(-\theta_{3},-\theta_{2})$. While for the points around
$[-\infty,-\theta_{3}]$,
$\alpha_{1}\approx\alpha_{2}\approx\alpha_{3}\approx\pm\pi/2$, i.e.
$\alpha\approx\pm\frac{3}{2}\pi$ which suggests a branch cut.

Figure 9: Contour integration

Figure 10: Contour limits as $R\to\infty$ and $\nu\to 0$

We will focus now on the limiting behaviour of contours in Figure (10) as
$R\to+\infty$ and $\nu\to 0$ while for the moment $\varepsilon$ is taken to be
very small so that the regions containing the poles $\theta_{i}$ resemble the
“dogbone” or “keyhole” types as in Figure (10). The closed contour $\Gamma$ in
Figure (10) is comprised of the vertical line $(-R,+R)$ and the elementary
parts
$\Gamma_{R},\Gamma_{-R},\Gamma_{\varepsilon},\Gamma_{\pm},\hat{\Gamma}_{\pm}$.
Applying the integration along the whole closed curve $\Gamma$, one can see
that the function is analytic while it does not contain any poles and
therefore the contour integration along this closed contour is zero. Now from
theorem 4.3.2 in Ablowitz and Fokas (2003) since the function $g(t)e^{st}$ is
exponentially bounded i.e. $\lim_{t\to\infty}g(t)e^{st}e^{\kappa}=0$ for some
$\kappa$, the curve integrals along the part-circle curves
$\Gamma_{R},\Gamma_{-R}$ disappear:

$\lim_{R\to\infty}\int\limits_{\Gamma_{R}\cup\,\Gamma_{-R}}g(t)e^{st}dt=0$

Since the integration on the regions $\hat{\Gamma}_{+}$ $\hat{\Gamma}_{-}$ is
along the opposing directions it then follows that the integral components
where the function values coincide cancel out:

$\lim_{\nu\to
0}\int\limits_{\hat{\Gamma}_{+}\cup\,\hat{\Gamma}_{-}}g(t)e^{st}dt=0$

Therefore the integration along the vertical line is opposite the integration
along the contours $\Gamma^{i}$. We can also easily see that without loss of
generality if some of these contours circle some single point representing
multiplicity of type 1.

For deriving the second result of Theorem 1, one can notice that the survival
function $P(X>x)$ is evaluated as

$P(X>x)=\int_{x}^{\infty}f_{\bm{\theta},\bm{\gamma}^{2}}(s)ds=\kappa\left(\frac{1}{2\pi\textbf{i}}\sum_{i=1}^{r}\int\limits_{\Gamma^{i}}\frac{g(t)}{t}e^{tx}dt+\frac{1}{2\pi\textbf{i}}\sum_{j=1}^{l}\int\limits_{\Gamma^{j}_{\varepsilon}}\frac{g(t)}{t}e^{tx}dt\right)$

which is obtained by replacing the expression for
$f_{\bm{\theta},\bm{\gamma}^{2}}(s)$ as in the first equation of Theorem 1
while $\frac{1}{t}$ appears by changing the order of integration such that
$\int_{x}^{\infty}g(t)e^{ts}ds=-\frac{g(t)}{t}e^{xt}$ for complex $t$ such
that $\text{Re}(t)<0.$

Note the fact that the integrand $\frac{g(t)}{t}$ has a simple pole at $0$
with residue
$\kappa^{-1}=\frac{e^{\sum_{i=1}^{p}\frac{\gamma_{i}^{2}}{4\theta_{i}}}}{\prod_{i=1}^{p}\theta_{i}^{\frac{n_{i}}{2}}}$.
The original integrating contour was the vertical line
$\textbf{i}\mathbb{R}+t_{0}$ at some $-min(\bm{\theta})<t_{0}<0$ which was
then reduced to the elementary contours $\Gamma_{i}$ as in Theorem 1.

Let us chose another integrating contour by changing $t_{0}$ to some
$t^{\prime}_{0}>0$. Now $\textbf{i}\mathbb{R}+t^{\prime}_{0}$ will be on the
right of the pole $0$ and the residue of $\frac{g(t)}{t}$ at $t=0$ is the
difference between these two integrals:

$\displaystyle\frac{1}{2\pi\textbf{i}}\int\limits_{\textbf{i}\mathbb{R}+t^{\prime}_{0}}\frac{g(t)}{t}e^{st}dt$
$\displaystyle-$
$\displaystyle\frac{1}{2\pi\textbf{i}}\int\limits_{\textbf{i}\mathbb{R}+t_{0}}\frac{g(t)}{t}e^{st}dt=\kappa^{-1}$
$\displaystyle\frac{\kappa}{2\pi\textbf{i}}\int\limits_{\textbf{i}\mathbb{R}+t^{\prime}_{0}}\frac{g(t)}{t}e^{st}dt$
$\displaystyle=$ $\displaystyle
1+\frac{\kappa}{2\pi\textbf{i}}\int\limits_{\textbf{i}\mathbb{R}+t_{0}}\frac{g(t)}{t}e^{st}dt.$
(18)

The appearance of the additional pole at $0$ can be seen as taking some
particular $\theta_{i}$ to be $0$ and with multiplicity $n_{i}=2$. In fact,
one can easily apply a shift parametrisation such that the pole $0$ is shifted
to some $-\theta_{0}$:

$\frac{1}{2\pi\textbf{i}}\int\limits_{\textbf{i}\mathbb{R}+t^{\prime}_{0}}\frac{g(t)}{t}e^{st}dt=e^{s\theta_{0}}\frac{1}{2\pi\textbf{i}}\int\limits_{\textbf{i}\mathbb{R}+t^{\prime}_{0}-\theta_{0}}\frac{g(t+\theta_{0})}{t+\theta_{0}}e^{st}dt=f_{\tilde{\bm{\theta}},\bm{0}}(x)e^{x\theta_{0}}\frac{\kappa}{\kappa^{\prime}}\quad\text{with}\quad\kappa^{\prime}=\frac{\theta_{0}\prod_{i=1}^{p}\theta_{i}^{\frac{n_{i}}{2}}}{e^{\sum_{i=1}^{p}\frac{\gamma_{i}^{2}}{4(\theta_{i}-\theta_{0})}}}$

While noting that for $t_{0}<0$,

$P(X\leq
x)=1-P(X>x)=1-\frac{\kappa}{2\pi\textbf{i}}\int_{x}^{+\infty}\int\limits_{\textbf{i}\mathbb{R}+t_{0}}g(t)e^{st}dtds=1+\frac{\kappa}{2\pi\textbf{i}}\int\limits_{\textbf{i}\mathbb{R}+t_{0}}g(t)\frac{e^{xt}}{t}dt$

identity (18) implies for the distribution function

$\displaystyle P(X\leq
x)=1-\int_{x}^{+\infty}f_{\bm{\theta},\bm{\gamma}^{2}}(s)ds=\frac{e^{x\theta_{0}}}{x\theta_{0}}f_{x\tilde{\bm{\theta}},x\tilde{\bm{\gamma}}^{2}}(1)\prod_{i=1}^{p}\frac{\theta_{i}^{\frac{n_{i}}{2}}}{(\theta_{i}+\theta_{0})^{\frac{n_{i}}{2}}}e^{\sum_{i=1}^{p}\frac{\gamma^{2}_{i}}{4(\theta_{i}+\theta_{0})}-\frac{\gamma^{2}_{i}}{4(\theta_{i})}}$
(19)

where $\tilde{\bm{\theta}}=(0,0,\bm{\theta})+\theta_{0}$ and
$\tilde{\bm{\gamma}}=(0,0,\bm{\gamma})$.

The second result follows from the discussion following equation (18) whose
derivation is now easily confirmed since the integration along the unbounded
contours $\textbf{i}\mathbb{R}+t_{0}$ and
$\textbf{i}\mathbb{R}+t^{\prime}_{0}$ is reduced to the elementary contours
$\Gamma_{i}$ and $\Gamma_{\epsilon}^{j}$ as described earlier in the proof the
theorem.

##  Appendix B Proof of Theorem 2

Note now some important observations related to the contour integration around
a simple pole. Theorem 4.3.1 of Ablowitz and Fokas (2003) states that if at a
particular pole $-\theta$, $g(t)(t+\theta)\to 0$ as $t\to-\theta_{i}$
then111If $\theta$ is a pole of finite order $n+1$ then the residue theorem
confirms  $\lim_{\varepsilon\to
0}\frac{1}{2\pi\textbf{i}}\int_{\Gamma_{\varepsilon}}g(t)dt=\frac{1}{n!}\frac{d^{n}}{d^{n}t}\left(g(\theta)(t-\theta)^{n}\right)$

$\lim_{\varepsilon\to 0}\int_{\Gamma_{\varepsilon}}g(t)dt=0$

In our case, since $\theta_{i}$ have multiplicity 1 the term in $g(t)$ related
to $\theta_{i}$ is $\frac{1}{(\theta_{i}+t)^{\frac{1}{2}}}$ and therefore
$g(t)(t+\theta)\to 0$ as $t\to-\theta_{i}$ and so $\Gamma_{\varepsilon}$
contours around the poles are zero and in addition:

$\int\limits_{\Gamma_{\varepsilon}}g(t)e^{st}dt\to 0\quad\nu,\varepsilon\to 0$

where the integration is applied at each $\Gamma_{\varepsilon}$ around
$-\theta_{i}$ and $\nu$ here represents the distance of the horizontal parts
of $\Gamma^{i}$ from the real axis (see Figure 4) . This implies that the
contours $\Gamma^{i}$ as in Figure 10 could be simplified to interval parts of
real axis while poles having residue zero. This limit holds since at each pole
$-\theta_{i}$, we have multiplicity 1 and therefore
$g(t)e^{st}(t+\theta_{i})\to 0$ as $t\to-\theta_{i}$. As a result the Bromwich
contour integral along the vertical line is the opposite of the value of the
contour integration along both horizontal contours which are on either sides
of the real axis of negative values. Now as $\nu\to 0$ both sides of
$\Gamma^{i}$’s can be taken as close as possible to the real line until they
both align with each other.

However, while the branch cut regions $C_{i}$, (bolded line regions in the
Figure 4) are the limits of $\Gamma^{i}$ as $\nu\to 0$, the argument values
$\alpha$ of the function $\frac{g(t)}{2\pi\mathbf{i}}e^{t}$ are such that
$\alpha\to\pm(2r-1)\frac{\pi}{2}$. And, more importantly as $\nu\to 0$, the
argument $\alpha$ approaches with opposing signs while the regions
$\Gamma_{\pm}$ go to $C_{r}$ . For example, $C_{1}$ in the figure is the limit
of $\Gamma^{1}_{+}$, as $\nu\to 0$, but $\alpha$ approaches the value
$\frac{\pi}{2}$ while in the region $\Gamma^{1}_{-}$, $\alpha$ approaches
$-\frac{\pi}{2}$. This feature repeats itself in every branch cut region
$C_{r}$ such that at regions $\Gamma_{+}$ $\Gamma_{-}$, $\alpha$ becomes
either $r\frac{\pi}{2}$ or $-r\frac{\pi}{2}$ respectively. While in such
regions the limiting behavior of $g(t)e^{st}$ has a sign change this is then
reversed by the opposing integration directions. This implies that the
integration along the branch cut is doubled while in the other parts or the
real axis is zero. Therefore after converting the elementary integration in
real coordinates and by introducing $-1$ in the square root the result
follows.

##  Appendix C Proof of Theorem 3

The density function of $Z$ for $z\geq 0$ is

$\displaystyle f_{Z}(z)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}f_{\bm{\theta},\bm{\gamma}^{2}}(y+z)f_{\bm{\theta}^{\prime},\bm{\gamma}^{\prime
2}}(y)dy$ $\displaystyle=$
$\displaystyle\kappa\kappa^{\prime}\int_{0}^{\infty}\int\limits_{\textbf{i}\mathbb{R}+t_{0}}\int\limits_{\textbf{i}\mathbb{R}+v_{0}}e^{(y+z)t}g(t)e^{yv}g^{\prime}(v)dvdtdy$
$\displaystyle=$
$\displaystyle\kappa\kappa^{\prime}\int\limits_{\textbf{i}\mathbb{R}+t_{0}}\int\limits_{\textbf{i}\mathbb{R}+v_{0}}e^{zt}g(t)dtg^{\prime}(v)dv\int_{0}^{\infty}e^{y(t+v)}dy$
$\displaystyle=$
$\displaystyle\frac{\kappa}{2\pi\textbf{i}}\frac{\kappa^{\prime}}{{2\pi\textbf{i}}}\int\limits_{\textbf{i}\mathbb{R}+t_{0}}\int\limits_{\textbf{i}\mathbb{R}+v_{0}}e^{zt}g(t)dtg^{\prime}(v)dv\frac{-1}{t+v}$
$\displaystyle=$
$\displaystyle\kappa\kappa^{\prime}\frac{1}{2\pi\textbf{i}}\int\limits_{\textbf{i}\mathbb{R}+t_{0}}e^{zt}g(t)g^{\prime}(-t)dt$

in the identities above we used the fact that
$\int_{0}^{\infty}e^{y(t+v)}dy=-\frac{1}{t+v}$ for complex numbers $t$ and $v$
which are both on the left of the line $\textbf{i}\mathbb{R}$ and from the
residue theorem
$\frac{1}{2\pi\textbf{i}}\int\limits_{\textbf{i}\mathbb{R}+v_{0}}\frac{1}{t+v}g^{\prime}(v)dv=-g^{\prime}(-t)$
for any $t$ on the right of $\textbf{i}\mathbb{R}+v_{0}$ (or the poles of
$g(t)$). As in the previous cases of Theorems 1 and 2, the integrating contour
could then be reduced to the poles or branch cuts on the left of the
integrating line which involve only those poles related to $\bm{\theta}$.
Additionally, by integrating on $(z,+\infty)$ one can obtain the right tail
probability so that

$P(Z>z>0)=\kappa\kappa^{\prime}\frac{-1}{2\pi\textbf{i}}\int\limits_{\textbf{i}\mathbb{R}+t_{0}}\frac{e^{zt}}{t}g(t)g^{\prime}(-t)dt$

which has the same effect as in Theorem 1 case where the additional pole at
$t=0$ on the left of the integrating contour leads to

$\kappa\kappa^{\prime}\frac{1}{2\pi\textbf{i}}\int\limits_{\textbf{i}\mathbb{R}+t_{0}}\frac{e^{zt}}{t}g(t)g^{\prime}(-t)dt+\kappa\kappa^{\prime}\frac{1}{2\pi\textbf{i}}\int\limits_{\textbf{i}\mathbb{R}+t^{\prime}_{0}}\frac{e^{zt}}{t}g(t)g^{\prime}(-t)dt=1$

for some $0<t^{\prime}_{0}<min(\bm{\theta^{\prime}})$. As a result for $z>0$

$P(Z<z)=\frac{e^{\theta_{0}z}}{\theta_{0}}f_{\tilde{Z}}(z)\prod_{i=1}^{p}\frac{\sqrt{\theta_{i}}}{\sqrt{\theta_{i}+\theta_{0}}}\prod_{i=1}^{p^{\prime}}\frac{\sqrt{\theta^{\prime}_{i}}}{\sqrt{\theta^{\prime}_{i}-\theta_{0}}}e^{\sum_{i=1}^{p}\frac{\gamma^{2}_{i}}{4\theta_{i}}-\frac{\gamma^{2}_{i}}{4(\theta_{i}+\theta_{0})}}e^{\sum_{i=1}^{p^{\prime}}\frac{\gamma^{\prime
2}_{i}}{4\theta^{\prime}_{i}}-\frac{\gamma^{\prime
2}_{i}}{4(\theta^{\prime}_{i}-\theta_{0})}}$

where

$\tilde{Z}=\frac{1}{2\theta_{0}}\chi^{2}_{2}(0)+\sum_{i=1}^{p}\frac{1}{2(\theta_{i}+\theta_{0})}\chi^{2}_{n_{i}}(\frac{\gamma_{i}^{2}}{2(\theta_{i}+\theta_{0})})-\sum_{i=1}^{p^{\prime}}\frac{1}{2(\theta^{\prime}_{i}-\theta_{0})}\chi^{2}_{n^{\prime}_{i}}(\frac{\gamma_{i}^{\prime
2}}{2(\theta^{\prime}_{i}-\theta_{0})})$

with $0<\theta_{0}<min(\bm{\theta}^{\prime})$. Noting that for $-Z=Y-X$, a
similar argument as above can be applied by switching the roles of $g$ and
$g^{\prime}$ so that for $z<0$ we obtain

$f_{z}(z)=\kappa\kappa^{\prime}\int\limits_{\textbf{i}\mathbb{R}+v_{0}}e^{-zt}g(-t)g^{\prime}(t)dt\quad
z<0$