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 be
mutually independent random variables that have gamma distributions with shape parameters
and scale parameters , . Then random variable is the convolution of independent gamma variables. Without loss of generality, the scale parameters ’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 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 variables. When , 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 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 and . In the case of , 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 and be, respectively, the distribution and density function of a gamma variable with shape , and scale . For the special case where , we use and , respectively. Let and , respectively, be the distribution and density function of in the case of . Finally, let be the Pochhammer polynomial (Abramowitz and Stegun, 1972, Equation 6.1.22).
2.1 Mathai’s Method
Mathai (1982) expresses the density of via a multiple infinite series
(1) |
where , , and is a confluent hypergeometric function of variables defined by a multiple series
a special function which has been studied in the literature (Mathai and Saxena, 1978). With the gamma function kernels, , 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 , the density is expressed in terms of the Kummer confluent hypergeometric function, as (Abramowitz and Stegun, 1972, Formula 13.1.2),
(2) | ||||
The benefit of Equation (2) is that the GSL (Galassi et al., 2009) has an implementation of . Note that the condition is not needed in (2). Indeed, if , then
where the third equation follows from (Abramowitz and Stegun, 1972, Equation 13.1.27).
The distribution function when can be explicitly expressed as, for ,
(3) | ||||
2.2 Moschopoulos’ Method
Moschopoulos (1985) expresses the density of by a single gamma series with coefficients that can be calculated recursively:
where , , , and is given by the recursive relations
with and for . This expression facilitates distribution function evaluation as
The weights
’s can be viewed as the probability masses of a discrete random variable on non-negative integers
(Vellaisamy and Upadhye, 2009). When , this discrete distribution is negative binomial. For , 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 .
Parameters | Density | Distribution Function | ||||
---|---|---|---|---|---|---|
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 |
We first compare the speed of the two methods in the case of . The shape parameters were set to be . The scale parameters were set to be . For each configuration, we used a large number () 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 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, , or the shape parameters bigger. Mathai’s method is much less sensitive to the parameter settings.
Parameters | Density | Distribution Function | |||||||
---|---|---|---|---|---|---|---|---|---|
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 |
Following the design and steps in the case of , we conducted a numerical analysis for . The shape parameters were set to be . The scale parameters were set to be . 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 , the multivariable confluent hypergeometric function has no efficient implementations yet to the best of our knowledge. Therefore, Mathai’s method needs to evaluate 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 times faster in density evaluation and times faster in distribution evaluation than Moschopoulos’ method.

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 be independent and identically distributed random variables following a mixture of two exponential distributions, and , with weights and , respectively. For any , define
(4) |
where, by convention, the summation over an empty set is . The process , , is a renewal process. The following result gives an expression of the distribution of in terms of the Kummer confluent hypergeometric function, which allows a very efficient numerical evaluation.
Proposition 1.
Proof.
First, note that
where and are independent of with and distributions, respectively. Next, using mixture randomization we get
where .
Similarly,
That is
(5) | ||||
Now it is sufficient to show that for any positive integers , and ,
(6) |
Firstly, from equation (3) we have
where the second equation follows from
Because of the identities
and
we finally obtain
∎
Corollary 1.
Let be independent and identically distributed random variables distributed as a mixture of exponential distributions, with weight , , and . For ,
where
, , and is the distribution function of the convolution of independent gamma variables with shape parameter and scale parameter .
Proof.
Note that
where are independent of with an distribution and
By substitution, we can complete proof. ∎
Parameters | Time Informations | ||||
---|---|---|---|---|---|
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 |
To demonstrate the computational efficiency of the result of Proposition 1, we performed a numerical study for . The scale parameters of the exponential distributions were set to be . For each configuration, we evaluated for and three values corresponding to the 20th, 50th, and 90th percentiles of a random sample of size 1,000 from the distribution of . The 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.
Parameters | Time Informations | |||||||
---|---|---|---|---|---|---|---|---|
Mathai | Mosch. | Approx. | Exact Value | Relative Error | ||||
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 |
For , we performed a similar study with Mathai’s and Moschopoulos’ method for evaluating using the result in Corollary 1. For comparison in accuracy and speed, the approximation method of Barnabani (2017) in evaluating was also included. The scale parameters of the exponential distributions were set to be
Again, for each configuration, we evaluated for three values corresponding to the 20th, 50th, and 90th percentiles of a random sample of size 1000 from the distribution of . 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 because of efficient GSL implementation of the univariate Kummer confluent hypergeometric function; when , 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 . 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 . It can be evaluated accurately and efficiently using the GSL implementation of the confluent hypergeometric function
. 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.