Density and Distribution Evaluation for Convolution of Independent Gamma Variables

1 Introduction

Fast and precise evaluation of the density and distribution function of convolution of independent gamma variables is important in many applications such as storage capacity measurement (e.g., Mathai, 1982), reliability analysis (e.g., Kadri et al., 2015), and point processes (e.g., Sim, 1992). Let $X_{1}, \dots, X_{n}$ be $n$

mutually independent random variables that have gamma distributions with shape parameters

α_{i} > 0

and scale parameters

β_{i} > 0

j = 1, \dots, n

. Then random variable

Y = \sum_{i = 1}^{n} X_{i}

is the convolution of independent gamma variables. Without loss of generality, the scale parameters

β_{i}

’s can be assumed to of distinctive values; otherwise, variables with the same scales can be summed before the convolution. Evaluation of the density and distribution of

Y

is the focus of this paper.

Exact evaluations, which do not have simple closed forms, are challenging. Mathai (1982) was the first to give an expression of the density in terms of multiple infinite series for the general case of arbitrary shape and scale parameters; the multiple series is a confluent hypergeometric function of $n - 1$ variables. When $n = 2$ , the confluent hypergeometric function is univariate with efficient implementations in the GNU Scientific Library (GSL) Galassi et al. (2009), which was adopted by Di Salvo (2006). Mathai (1982) also gave relatively easier-to-compute expressions in the special cases where all shapes are integers or are identical. Moschopoulos (1985) simplified the complex expression into a single gamma-series representation along with a formula for the truncation error to evaluate the precision of numerical computation. Akkouchi (2005) indicated that the density of $Y$ can be expressed through an integral of the generalized beta function, but did not give explicit ways to numerically evaluate the integral. More recently, Vellaisamy and Upadhye (2009) proposed a random parameter representation of the gamma-series of Moschopoulos (1985), where the weights define the probability mass function of a discrete distribution on non-negative integers. Implementations of the gamma-series methods are simple but the computation is too much CPU-time consuming when the variability of the scale parameters is large and the shape parameters are small. Built on the representation of Vellaisamy and Upadhye (2009), Barnabani (2017) proposed a fast approximation which approximates the weights of the gamma-series by a discrete distribution.

The contribution of this article is two-fold. First, we give a computationally review of the methods of Mathai (1982) and Moschopoulos (1985) and provide their implementations in our R package coga (Hu et al., 2017). Their speeds are compared in cases of $n = 2$ and $n = 3$ . In the case of $n = 3$ , the accuracy of the fast approximation of Barnabani (2017)

is also assessed using our implementation. Second, in an application to renewal processes with holding times following a mixture of exponential distributions, we derive a new formula for the probability mass function of the number of renewals by a given amount of time, which provides very fast exact evaluations in a numerical study.

2 Exact Evaluations

Let us introduce some notations first. Let $G (y; α, β)$ and $g (y; α, β)$ be, respectively, the distribution and density function of a gamma variable with shape $α$ , and scale $β$ . For the special case where $β = 1$ , we use $G (y; α)$ and $g (y; α)$ , respectively. Let $F (x; (α_{1}, β_{1}), (α_{2}, β_{2}))$ and $f (x; (α_{1}, β_{1}), (α_{2}, β_{2}))$ , respectively, be the distribution and density function of $Y$ in the case of $n = 2$ . Finally, let $(x)_{m} = x (x + 1) \dots (x + m - 1)$ be the Pochhammer polynomial (Abramowitz and Stegun, 1972, Equation 6.1.22).

2.1 Mathai’s Method

Mathai (1982) expresses the density of $Y$ via a multiple infinite series

f (x) = {[n \prod j = 1 β_{j}^{α_{j}} Γ (γ)]}^{- 1} x^{γ - 1} e^{- x / β_{1}} ϕ (α_{2}, \dots, α_{n}; γ; (1 / β_{1} - 1 / β_{2}) x, \dots, (1 / β_{1} - 1 / β_{n}) x),

(1)

where $β_{1} = {min}_{j} (β_{j})$ , $γ = \sum_{j = 1}^{n} α_{j}$ , and $ϕ$ is a confluent hypergeometric function of $n - 1$ variables defined by a multiple series

		$ϕ (α_{2}, \dots, α_{n}; γ; (1 / β_{1} - 1 / β_{2}) x, \dots, (1 / β_{1} - 1 / β_{n}) x)$
	$=$	$\infty \sum r_{2} = 0 \dots \infty \sum r_{n} = 0 {(α_{2})_{r_{2}} \dots (α_{n})_{r_{n}} {[(1 / β_{1} - 1 / β_{2}) x]}_{1}^{r_{2}} \dots {[(1 / β_{1} - 1 / β_{n}) x]}_{1}^{r_{n}} / [r_{2}! \dots r_{n}! (γ)_{r}]},$

a special function which has been studied in the literature (Mathai and Saxena, 1978). With the gamma function kernels, $x^{γ + r_{i} - 1} e^{- x / β_{1}}$ , The distribution function can be expressed in terms of incomplete gamma functions by term-by-term integration of Equation (1).

For the special case of $n = 2$ , the density is expressed in terms of the Kummer confluent hypergeometric function, $_{1} F_{1}$ as (Abramowitz and Stegun, 1972, Formula 13.1.2),

	$f (x; (α_{1}, β_{1}), (α_{2}, β_{2}))$	$= \frac{x^{γ - 1} e^{- x / β_{1}}}{β_{1}^{α_{1}} β_{2}^{α_{2}} Γ (γ)}_{1} F_{1} (α_{2}; γ; (1 / β_{1} - 1 / β_{2}) x)$		(2)
		$= {(\frac{β_{1}}{β_{2}})}^{α_{2}} g (x; γ, β_{1})_{1} F_{1} (α_{2}; γ; (1 / β_{1} - 1 / β_{2}) x) .$		(2)

The benefit of Equation (2) is that the GSL (Galassi et al., 2009) has an implementation of $_{1} F_{1}$ . Note that the condition $β_{1} < β_{2}$ is not needed in (2). Indeed, if $β_{1} > β_{2}$ , then

	$f (x; (α_{1}, β_{1}), (α_{2}, β_{2}))$	$= f (x; (α_{2}, β_{2}), (α_{1}, β_{1}))$
		$= \frac{x^{γ - 1} e^{- x / β_{2}}}{β_{1}^{α_{1}} β_{2}^{α_{2}} Γ (γ)}_{1} F_{1} (α_{1}; γ; (1 / β_{2} - 1 / β_{1}) x)$
		$= \frac{x^{γ - 1} e^{- x / β_{2}}}{β_{1}^{α_{1}} β_{2}^{α_{2}} Γ (γ)} e^{- (1 / β_{2} - 1 / β_{1}) x}_{1} F_{1} (α_{2}; γ; (1 / β_{1} - 1 / β_{2}) x)$
		$= \frac{x^{γ - 1} e^{- x / β_{1}}}{β_{1}^{α_{1}} β_{2}^{α_{2}} Γ (γ)}_{1} F_{1} (α_{2}; γ; (1 / β_{1} - 1 / β_{2}) x),$

where the third equation follows from $_{1} F_{1} (a; b; z) \equiv e^{z}_{1} F_{1} (b - a; b; - z)$ (Abramowitz and Stegun, 1972, Equation 13.1.27).

The distribution function when $n = 2$ can be explicitly expressed as, for $y > 0$ ,

	$F (y; (α_{1}, β_{1}), (α_{2}, β_{2}))$	(3)
$=$	$\frac{1}{β_{1}^{α_{1}} β_{2}^{α_{2}}} \infty \sum k = 0 \frac{(\frac{α_{2} + k - 1}{k})}{Γ (γ + k)} (1 / β_{1} - 1 / β_{2})^{k} \int_{0}^{y} x^{k + γ - 1} e^{- x / β_{1}} d x$
$=$	$\frac{1}{β_{1}^{α_{1}} β_{2}^{α_{2}}} \infty \sum k = 0 \frac{(\frac{α_{2} + k - 1}{k})}{Γ (γ + k)} (1 / β_{1} - 1 / β_{2})^{k} β_{1}^{k + γ} G (y / β_{1}; k + γ) Γ (k + γ)$
$=$	${(\frac{β_{1}}{β_{2}})}^{α_{2}} \infty \sum k = 0 (\frac{α_{2} + k - 1}{k}) (1 - β_{1} / β_{2})^{k} G (y / β_{1}; k + γ) .$

2.2 Moschopoulos’ Method

Moschopoulos (1985) expresses the density of $Y$ by a single gamma series with coefficients that can be calculated recursively:

	$f (x) =$
	$=$	$C \infty \sum k = 0 δ_{k} g (x; ρ + k, β_{1}), x > 0,$

where $β_{1} = {min}_{i} (β_{i})$ , $C = \prod_{i = 1}^{n} (β_{1} / β_{i})^{α_{i}}$ , $ρ = \sum_{i = 1}^{n} α_{i} > 0$ , and $δ_{k}$ is given by the recursive relations

δ_{k + 1} = \frac{1}{k + 1} k + 1 \sum i = 1 i γ_{i} δ_{k + 1 - i}, k = 0, 1, 2, \dots,

with $δ_{0} = 1$ and $γ_{k} = \sum_{i = 1}^{n} α_{i} (1 - β_{1} / β_{i})^{k} / k$ for $k = 1, 2, \dots$ . This expression facilitates distribution function evaluation as

F (y) = C \infty \sum k = 0 δ_{k} G (y; ρ + k, β_{1}), y > 0.

The weights $C δ_{k}$

’s can be viewed as the probability masses of a discrete random variable on non-negative integers

(Vellaisamy and Upadhye, 2009). When

n = 2

, this discrete distribution is negative binomial. For

n > 2

, Barnabani (2017)

proposed to approximate the discrete distribution by a three-parameter generalized negative binomial distribution defined by

Jain and Consul (1971)

through moment matching.

2.3 Timing Comparison

We implemented the methods of Mathai (1982) and Moschopoulos (1985) in an open source R package coga (Hu et al., 2017). The computation is done in C++ code and interfaced to R (R Core Team, 2017) in the coga package. In addition, the fast approximation of Barnabani (2017) is also available in the package for $n > 2$ .

Parameters			Density		Distribution Function
$α$	$β_{1}$	$β_{2}$	Moschopoulos	Mathai	Moschopoulos	Mathai
0.2	0.4	0.3	23,030	86	25,011	2,648
	4	0.3	103,987	176	110,804	9,018
	4	3	24,566	88	25,696	2,849
2	0.4	0.3	29,163	94	31,030	3,086
	4	0.3	165,901	96	173,916	13,590
	4	3	30,378	90	33,489	3,397
20	0.4	0.3	53,231	103	57,468	6,471
	4	0.3	538,807	175	566,857	37,527
	4	3	53,259	108	58,479	6,604

Table 1: Timing comparison (in microseconds) of Mathai’s and Moschopoulos’ methods when

n = 2

in evaluating the density and distribution function of convolutions of independent gamma variables.

We first compare the speed of the two methods in the case of $n = 2$ . The shape parameters were set to be $α_{1} = α_{2} \in {0.2, 2, 20}$ . The scale parameters were set to be $(β_{1}, β_{2}) \in {(0.4, 0.3), (4, 0.3), (4, 3)}$ . For each configuration, we used a large number ( $100, 000$ ) of simulated observations from the distribution to determine the bulk range of the observations. Then we evaluate the density and the distribution of the convolution over 100 equally spaced grid points in the bulk range. The evaluations were repeated 100 times.

Table 1 summarizes the median time to evaluate density and distribution at the 100-point grid from 100 replicates obtained on an Intel 2.50GHz computer. Density evaluation using Mathai’s method (implemented with the $_{1} F_{1}$ function from the GSL) performs much faster than Moschopoulos’ method in all settings; in some settings, it is up to 3,000 times faster Distribution evaluation takes much longer than density evaluation using Mathai’s method, but is still up to 16 times faster than Moschopoulos’ method. Moschopoulos’ method takes longer when the scale parameters are very different, $(β_{1}, β_{2}) = (4, 0.3)$ , or the shape parameters bigger. Mathai’s method is much less sensitive to the parameter settings.

Parameters				Density			Distribution Function
$α$	$β_{1}$	$β_{2}$	$β_{3}$	Mosch.	Mathai	Approx.	Mosch.	Mathai	Approx.
0.2	0.4	0.3	0.2	37	1,223	6	39	1,428	8
	4	0.3	0.2	167	5,796	18	181	8,067	28
	4	3	0.2	186	10,208	19	197	14,000	30
	4	3	2	38	1,245	6	40	1,474	8
2	0.4	0.3	0.2	48	1,044	8	51	1,328	11
	4	0.3	0.2	242	5,333	21	253	8,975	35
	4	3	0.2	313	12,597	25	331	20,646	43
	4	3	2	49	1,082	8	53	1,415	11
20	0.4	0.3	0.2	109	3,250	14	119	4,721	22
	4	0.3	0.2	780	16,083	29	950	40,704	78
	4	3	0.2	596	21,553	18	1,418	133,699	101
	4	3	2	110	3,329	14	123	4,953	23

Table 2: Timing comparison (in milliseconds) of Mathai’s, Moschopoulos’ exact methods and Barnabani’s approximation method when

n = 3

Following the design and steps in the case of $n = 2$ , we conducted a numerical analysis for $n = 3$ . The shape parameters were set to be $α_{1} = α_{2} = α_{3} \in {0.2, 2, 20}$ . The scale parameters were set to be $(β_{1}, β_{2}, β_{3}) \in {(0.4, 0.3, 0.2), (4, 0.3, 0.2), (4, 3, 0.2), (4, 3, 2)}$ . The two methods gave numerically indistinguishable results in both density and distribution evaluations, verifying each other. Table 2 summarize the median timing results from 100 replicates. When $n > 2$ , the multivariable confluent hypergeometric function has no efficient implementations yet to the best of our knowledge. Therefore, Mathai’s method needs to evaluate $n - 1$ nested infinite series, which makes it quite complicated (Jasiulewicz and Kordecki, 2003; Sen and Balakrishnan, 1999). In this case, Moschopoulos’ method is preferred for its single-series. The approximation of Barnabani (2017) was also included in the comparison, which is up to $30$ times faster in density evaluation and $14$ times faster in distribution evaluation than Moschopoulos’ method.

Figure 1: Differences between the approximation method and the exact methods in evaluating the density and distribution function of convolution of three independent gamma variables. The parameter in setup 1 are $α_{1} = α_{2} = α_{3} = 0.2$ , $β_{1} = 4$ , $β_{2} = 0.3$ , and $β_{3} = 0.2$ . The parameter in setup 2 are $α_{1} = α_{2} = α_{3} = 2$ , $β_{1} = 0.4$ , $β_{2} = 0.3$ , and $β_{3} = 0.2$ .

The exact implementations make it possible to assess the accuracy of the approximation approach of Barnabani (2017). The differences between the approximation and exact evaluation in the two worst cases out of a collection of settings we experimented are shown in Figure 1. Apparently, the approximation is very accurate.

3 Application to a Renewal Process

Let ${U_{k}}_{k \geq 1}$ be independent and identically distributed random variables following a mixture of two exponential distributions, $E x p (β_{1})$ and $E x p (β_{2})$ , with weights $p \in (0, 1)$ and $1 - p$ , respectively. For any $t > 0$ , define

N (t) = sup {n \geq 0 : n \sum k = 1 U_{k} \leq t},

(4)

where, by convention, the summation over an empty set is $0$ . The process $N (t)$ , $t > 0$ , is a renewal process. The following result gives an expression of the distribution of $N (t)$ in terms of the Kummer confluent hypergeometric function, which allows a very efficient numerical evaluation.

Proposition 1.

For the renewal process defined in (4) and integer $n \geq 0$ ,

	$Pr (N (t) = n) = n \sum k = 0$	${p ψ (β_{1}, β_{2}, t, k, n) + (1 - p) ψ (β_{2}, β_{1}, t, n - k, n)} \times$
		$\frac{t^{n}}{β_{1}^{k} β_{2}^{n - k} Γ (n + 1)} (\frac{n}{k}) p^{k} (1 - p)^{n - k},$

where

ψ (β_{1}, β_{2}, t, k, n) = e^{- t / β_{1}}_{1} F_{1} (n - k; n + 1; t (1 / β_{1} - 1 / β_{2})) .

Proof.

First, note that

		$Pr (N (t) = n) = Pr (n \sum k = 1 U_{k} \leq t, n + 1 \sum k = 1 U_{k} > t)$
	$=$	$p Pr (n \sum k = 1 U_{k} \leq t, n \sum k = 1 U_{k} + E_{1} > t) + (1 - p) Pr (n \sum k = 1 U_{k} \leq t, n \sum k = 1 U_{k} + E_{2} > t),$

where $E_{1}$ and $E_{2}$ are independent of ${U_{k}}_{k \geq 1}$ with $E x p (β_{1})$ and $E x p (β_{2})$ distributions, respectively. Next, using mixture randomization we get

	$Pr (n \sum k = 1 U_{k} \leq t, n \sum k = 1 U_{k} + E_{1} > t)$	$= Pr (n \sum k = 1 U_{k} \leq t) - Pr (n \sum k = 1 U_{k} + E_{1} \leq t)$
		$= n \sum k = 0 H (t; (k, β_{1}), (n - k, β_{2})) (\frac{n}{k}) p^{k} (1 - p)^{n - k},$

where $H (x; α_{1}, β_{1}, α_{2}, β_{2}) = F (x; α_{1}, β_{1}, α_{2}, β_{2}) - F (x; α_{1} + 1, β_{1}, α_{2}, β_{2})$ .

Similarly,

	$Pr (n \sum k = 1 U_{k} \leq t, n \sum k = 1 U_{k} + E_{2} > t)$	$= Pr (n \sum k = 1 U_{k} \leq t) - Pr (n \sum k = 1 U_{k} + E_{2} \leq t)$
		$= n \sum k = 0 H (t; (n - k, β_{2}), (k, β_{1})) (\frac{n}{k}) p^{k} (1 - p)^{n - k} .$

That is

		$Pr (N (t) = n)$		(5)
	$=$	$n \sum k = 0 [p H (t; (k, β_{1}), (n - k, β_{2})) + (1 - p) H (t; (n - k, β_{2}), (k, β_{1}))] (\frac{n}{k}) p^{k} (1 - p)^{n - k}$		(5)

Now it is sufficient to show that for any positive integers $α_{1}$ , $α_{2}$ and $y, β_{1}, β_{2} > 0$ ,

H (y; (α_{1}, β_{1}), (α_{2}, β_{2})) = \frac{y^{α_{1} + α_{2}} e^{- y / β_{1}}}{β_{1}^{α_{1}} β_{2}^{α_{2}} Γ (α_{1} + α_{2} + 1)}_{1} F_{1} (α_{2}; α_{1} + α_{2} + 1; y (1 / β_{1} - 1 / β_{2})) .

(6)

Firstly, from equation (3) we have

	$H (y; (α_{1}, β_{1}), (α_{2}, β_{2})) =$	${(\frac{β_{1}}{β_{2}})}^{α_{2}} \infty \sum k = 0 (\frac{α_{2} + k - 1}{k}) {(1 - β_{1} / β_{2})}_{1}^{k} G (y / β_{1}; k + α_{1} + α_{2})$
		$- {(\frac{β_{1}}{β_{2}})}^{α_{2}} \infty \sum k = 0 (\frac{α_{2} + k - 1}{k}) {(1 - β_{1} / β_{2})}_{1}^{k} G (y / β_{1}; k + α_{1} + α_{2} + 1)$
	$=$	${(\frac{β_{1}}{β_{2}})}^{α_{2}} \infty \sum k = 0 (\frac{α_{2} + k - 1}{k}) (1 - β_{1} / β_{2})^{k} \frac{(y / β_{1})^{k + α_{1} + α_{2}} e^{- y / β_{1}}}{Γ (k + α_{1} + α_{2} + 1)}$
	$=$	${(\frac{β_{1}}{β_{2}})}^{α_{2}} {(\frac{y}{β_{1}})}^{α_{1} + α_{2}} e^{- y / β_{1}} \infty \sum k = 0 \frac{(\frac{α_{2} + k - 1}{k}) {[y (1 / β_{1} - 1 / β_{2})]}_{1}^{k}}{Γ (k + α_{1} + α_{2} + 1)},$

where the second equation follows from

G (y; α) - G (y; α + 1) = \frac{y^{α} e^{- y}}{Γ (α + 1)} .

Because of the identities

(\frac{α_{2} + k + 1}{k}) = \frac{(α_{2})_{k}}{k!}

and

Γ (k + α_{1} + α_{2} + 1) = (α_{1} + α_{2} + 1)_{k} Γ (α_{1} + α_{2} + 1),

we finally obtain

	$H (y; (α_{1}, β_{1}), (α_{2}, β_{2})) =$	$\frac{y^{α_{1} + α_{2}} e^{- y / β_{1}}}{β_{1}^{α_{1}} β_{2}^{α_{2}} Γ (α_{1} + α_{2} + 1)} \infty \sum k = 0 \frac{(α_{2})_{k} [y (1 / β_{1} - 1 / β_{2})]^{k}}{(α_{1} + α_{2} + 1)_{k} k!}$
	$=$	$\frac{y^{α_{1} + α_{2}} e^{- y / β_{1}}}{β_{1}^{α_{1}} β_{2}^{α_{2}} Γ (α_{1} + α_{2} + 1)}_{1} F_{1} (α_{2}; α_{1} + α_{2} + 1; y (1 / β_{1} - 1 / β_{2})) .$

∎

Corollary 1.

Let ${U_{k}}_{k \geq 1}$ be independent and identically distributed random variables distributed as a mixture of $S$ exponential distributions, $E x p (β_{s})$ with weight $p_{s}$ , $s = 1, \dots, S$ , and $\sum_{s = 1}^{S} p_{s} = 1$ . For $n \geq 0$ ,

		$Pr (N (t) = n)$
	$=$	$n \sum k_{1} = 0 n - k_{1} \sum k_{2} = 0 \dots n - k_{1} - \dots - k_{S - 2} \sum k_{S - 1} = 0 [(S \sum s = 1 p_{s} H_{s} (t; (k_{1}, β_{1}), \dots, (k_{S}, β_{S}))) \frac{n!}{k_{1}! k_{2}! \dots k_{S}!} p_{1}^{k_{1}} \dots p_{S}^{k_{S}}],$

where

	$H_{s} (t; (k_{1}, β_{1}), \dots, (k_{S}, β_{S})) =$	$F_{S} (t; (k_{1}, β_{1}), \dots, (k_{s}, β_{s}), \dots, (k_{S}, β_{S}))$
		$- F_{S} (t; (k_{1}, β_{1}), \dots, (k_{s} + 1, β_{s}), \dots, (k_{S}, β_{S})),$

$s = 1, \dots, S$ , $k_{S} = n - k_{1} - \dots - k_{S - 1}$ , and $F_{S}$ is the distribution function of the convolution of $S$ independent gamma variables with shape parameter $(k_{1}, \dots, k_{S})$ and scale parameter $(β_{1}, \dots, β_{S})$ .

Proof.

Note that

Pr (N (t) = n) = Pr (n \sum k = 1 U_{k} \leq t, n + 1 \sum k = 1 U_{k} > t) = S \sum s = 1 p_{s} Pr (n \sum k = 1 U_{k} \leq t, n \sum k = 1 U_{k} + E_{s} > t),

where $E_{s}$ are independent of ${U_{k}}_{k \geq 1}$ with an $E x p (β_{s})$ distribution and

		$Pr (n \sum k = 1 U_{k} \leq t, n \sum k = 1 U_{k} + E_{s} > t) = Pr (n \sum k = 1 U_{k} \leq t) - Pr (n \sum k = 1 U_{k} + E_{s} \leq t)$
	$=$	$n \sum k_{1} = 0 n - k_{1} \sum k_{2} = 0 \dots n - k_{1} - \dots - k_{S - 2} \sum k_{S - 1} = 0 H_{s} (t, (k_{1}, β_{1}), \dots, (k_{S}, β_{S})) \frac{n!}{k_{1}! k_{2}! \dots k_{S}!} p_{1}^{k_{1}} \dots p_{S}^{k_{S}} .$

By substitution, we can complete proof. ∎

Parameters			Time Informations
$β_{1}$	$β_{2}$	$n$	Mathai	Moschopoulos	Proposed
0.4	0.3	27	4,993	36,314	33
		32	5,788	43,129	34
		40	6,981	54,637	41
4	0.3	10	2,675	23,613	15
		18	4,667	42,529	23
		30	7,525	65,408	58
4	3	2	308	951	9
		3	380	1,497	15
		5	556	3,226	16

Table 3: Timing comparison (in microseconds) using different methods in

H

to evaluate

Pr (N (10) = n)

in the renewal process application when

S = 2

To demonstrate the computational efficiency of the result of Proposition 1, we performed a numerical study for $S = 2$ . The scale parameters of the exponential distributions were set to be $(β_{1}, β_{2}) \in {(0.4, 0.3), (4, 0.3), (4, 3)}$ . For each configuration, we evaluated $Pr (N (t) = n}$ for $t = 10$ and three $n$ values corresponding to the 20th, 50th, and 90th percentiles of a random sample of size 1,000 from the distribution of $N (10)$ . The $H$ function was evaluated with Mathai’s method, Moschopoulos’ method and the proposed result in Proposition 1. The numerical results are identical and the median timing results from 100 replicates are summarized in Table 3. The method in Proposition 1, which bypasses evaluating two distribution function of the convolution of two exponential variables, shows a huge advantage, speeding up Mathai’s method by a factor of about 60–200.

$β_{1}$	$β_{2}$	$β_{3}$	$n$	Mathai	Mosch.	Approx.	Exact Value	Relative Error
Parameters				Time Informations
0.4	0.3	0.2	36	95,440	3,033	669	4.2456e-02	3.3887e-03
			42	132,284	4,097	895	5.7594e-02	-2.5825e-03
			51	189,053	5,901	1,285	2.2793e-02	-4.7991e-04
4	0.3	0.2	10	8,729	376	63	2.8303e-02	3.6682e-03
			19	33,312	1,249	218	3.3972e-02	-2.9601e-05
			35	113,068	3,906	674	1.4896e-02	-1.0625e-02
4	3	0.2	5	2,791	94	15	5.8889e-02	3.3020e-04
			10	12,937	382	62	6.2835e-02	-2.0693e-03
			19	49,028	1,229	203	2.1189e-02	1.5217e-03
4	3	2	2	31	6	1	1.2854e-01	1.2242e-03
			4	233	26	3	1.8740e-01	-1.9957e-03
			7	829	62	11	7.2131e-02	1.9813e-03

Table 4: Timing comparison (in milliseconds) using different methods in

H

to evaluate

Pr (N (10) = n)

in the renewal process application when

S = 3

For $S > 2$ , we performed a similar study with Mathai’s and Moschopoulos’ method for evaluating $H_{3}$ using the result in Corollary 1. For comparison in accuracy and speed, the approximation method of Barnabani (2017) in evaluating $H_{3}$ was also included. The scale parameters of the exponential distributions were set to be

(β_{1}, β_{2}, β_{3}) \in {(0.4, 0.3, 0.2), (4, 0.3, 0.2), (4, 3, 0.2), (4, 3, 2)} .

Again, for each configuration, we evaluated $Pr (N (10) = n)$ for three $n$ values corresponding to the 20th, 50th, and 90th percentiles of a random sample of size 1000 from the distribution of $N (10)$ . The median timing results from 100 replicates and the relative error of the approximation method are summarized in Table 4. Similar to the comparison reported in Section 2, Moschopoulos’ method is over 10 times faster than Mathai’s method in all the cases. A small sumation of $β$ ’s requires longer computaion time. The approximation method is over 6 times faster than Moschopoulo’s method, with relative error less than 1 percent.

4 Conclusions

We reviewed two exact methods and one approximation method for evaluating the density and distribution of convolutions of independent gamma variables. From our study, the method of Mathai (1982) is the fastest in the case of $n = 2$ because of efficient GSL implementation of the univariate Kummer confluent hypergeometric function; when $n > 3$ , the method of Moschopoulos (1985) is faster than Mathai’s method. The fast approximation method Barnabani (2017) is quite accurate, which provides a useful tool for applications with $n > 3$ . Implementations of the reviewed methods are available in R package coga (Hu et al., 2017), which is built on C++ code for fast speed.

The result in Proposition 1 for the distribution of the event count in the renew process application with holding time following a mixture of exponential distribution is an extremely fast evaluation. One side-benefit of the proposition is formula (6) for the difference of distribution functions of convolution of two independent Erlang distributions with shape parameters that differ by $1$ . It can be evaluated accurately and efficiently using the GSL implementation of the confluent hypergeometric function $_{1} F_{1}$

. This difference plays an important role in computation of distribution of the occupation times for a certain continuous time Markov chain

(Pozdnyakov et al., 2017), the defective density in generalized integrated telegraph processes (Zacks, 2004), and the first-exit time in a compound process (Perry et al., 1999).

References

Abramowitz and Stegun (1972) Abramowitz, M. and I. A. Stegun (1972). Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Volume 9. Dover, New York.
Akkouchi (2005) Akkouchi, M. (2005). On the convolution of gamma distributions. Soochow Journal of Mathematics 31(2), 205–211.
Barnabani (2017) Barnabani, M. (2017). An approximation to the convolution of gamma distributions. Communications in Statistics - Simulation and Computation 46(1), 331–343.
Di Salvo (2006) Di Salvo, F. (2006). The exact distribution of the weighted convolution of two gamma distributions.
Galassi et al. (2009) Galassi, M., J. Davies, J. Theiler, B. Gough, and G. Jungman (2009). GNU Scientific Library: Reference Manual (3 ed.). Network Theory Ltd.
Hu et al. (2017) Hu, C., V. Pozdnyakov, and J. Yan (2017). coga: Convolution of Gamma Distributions. R package version 0.2.1.
Jain and Consul (1971) Jain, G. C. and P. C. Consul (1971). A generalized negative binomial distribution. SIAM Journal on Applied Mathematics 21(4), 501–513.
Jasiulewicz and Kordecki (2003) Jasiulewicz, H. and W. Kordecki (2003). Convolutions of Erlang and of pascal distributions with applications to reliability. Demonstratio Mathematica. Warsaw Technical University Institute of Mathematics 36(1), 231–238.
Kadri et al. (2015) Kadri, T., K. Smaili, and S. Kadry (2015). Markov modeling for reliability analysis using hypoexponential distribution. In S. Kadry and A. El Hami (Eds.), Numerical Methods for Reliability and Safety Assessment: Multiscale and Multiphysics Systems, pp. 599–620. Cham: Springer.
Mathai (1982) Mathai, A. M. (1982). Storage capacity of a dam with gamma type inputs. Annals of the Institute of Statistical Mathematics 34(1), 591–597.
Mathai and Saxena (1978) Mathai, A. M. and R. K. Saxena (1978). The H function with Applications in Statistics and Other Disciplines. New Delhi: Wiley.
Moschopoulos (1985) Moschopoulos, P. G. (1985). The distribution of the sum of independent gamma random variables. Annals of the Institute of Statistical Mathematics 37(1), 541–544.
Perry et al. (1999) Perry, D., W. Stadje, and S. Zacks (1999). First-exit times for increasing compound processes. Communications in Statistics. Stochastic Models 15(5), 977–992.
Pozdnyakov et al. (2017) Pozdnyakov, V., C. Hu, T. Meyer, and J. Yan (2017). On discretely observed Brownian motion governed by a continuous time Markov chain. Technical report, University of Connecticut, Department of Statistics.
R Core Team (2017) R Core Team (2017). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing.
Sen and Balakrishnan (1999) Sen, A. and N. Balakrishnan (1999). Convolution of geometrics and a reliability problem. Statistics & Probability Letters 43(4), 421–426.
Sim (1992) Sim, C. (1992). Point processes with correlated gamma interarrival times. Statistics & Probability Letters 15(2), 135–141.
Vellaisamy and Upadhye (2009) Vellaisamy, P. and N. Upadhye (2009). On the sums of compound negative binomial and gamma random variables. Journal of Applied Probability 46(1), 272–283.
Zacks (2004) Zacks, S. (2004). Generalized integrated telegraph processes and the distribution of related stopping times. Journal of Applied Probability 41(2), 497–507.

Density and Distribution Evaluation for Convolution of Independent Gamma Variables

Authors

Entropies of sums of independent gamma random variables

On the discrete analog of gamma-Lomax distribution: properties and applications

Deep-Learned Event Variables for Collider Phenomenology

A series acceleration algorithm for the gamma-Pareto (type I) convolution and related functions of interest for pharmacokinetics

A general statistical model for waiting times until collapse of a system

Gamma Spaces and Information

A Cluster Model for Growth of Random Trees

1 Introduction

2 Exact Evaluations

2.1 Mathai’s Method

2.2 Moschopoulos’ Method

2.3 Timing Comparison

3 Application to a Renewal Process

Proposition 1.

Proof.

Corollary 1.

Proof.

4 Conclusions

References

Density and Distribution Evaluation for Convolution of Independent Gamma Variables

Authors

Related Research

Entropies of sums of independent gamma random variables

On the discrete analog of gamma-Lomax distribution: properties and applications

Deep-Learned Event Variables for Collider Phenomenology

A series acceleration algorithm for the gamma-Pareto (type I) convolution and related functions of interest for pharmacokinetics

A general statistical model for waiting times until collapse of a system

Gamma Spaces and Information

A Cluster Model for Growth of Random Trees

This week in AI

1 Introduction

2 Exact Evaluations

2.1 Mathai’s Method

2.2 Moschopoulos’ Method

2.3 Timing Comparison

3 Application to a Renewal Process

Proposition 1.

Proof.

Corollary 1.

Proof.

4 Conclusions

References