Density and Distribution Evaluation for Convolution of Independent Gamma Variables

by   Chaoran Hu, et al.

Several numerical evaluations of the density and distribution of convolution of independent gamma variables are compared in their accuracy and speed. In application to renewal processes, an efficient formula is derived for the probability mass function of the event count.



page 1

page 2

page 3

page 4


Entropies of sums of independent gamma random variables

We establish several Schur-convexity type results under fixed variance f...

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

A two parameter discrete gamma-Lomax distribution is derived as a discre...

Deep-Learned Event Variables for Collider Phenomenology

The choice of optimal event variables is crucial for achieving the maxim...

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

The gamma-Pareto type I convolution (GPC type I) distribution, which has...

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

The distribution of waiting times until the occurrence of a critical eve...

Gamma Spaces and Information

We investigate the role of Segal's Gamma-spaces in the context of classi...

A Cluster Model for Growth of Random Trees

We first consider the growth of trees by probabilistic attachment of new...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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


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),


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 ,


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
Table 1: Timing comparison (in microseconds) of Mathai’s and Moschopoulos’ methods when 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 . 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
Table 2: Timing comparison (in milliseconds) of Mathai’s, Moschopoulos’ exact methods and Barnabani’s approximation method when .

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.

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 , , , and . The parameter in setup 2 are , , , and .

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


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.

For the renewal process defined in (4) and integer ,



First, note that

where and are independent of with and distributions, respectively. Next, using mixture randomization we get

where .


That is


Now it is sufficient to show that for any positive integers , and ,


Firstly, from equation (3) we have

where the second equation follows from

Because of the identities


we finally obtain

Corollary 1.

Let be independent and identically distributed random variables distributed as a mixture of exponential distributions, with weight , , and . For ,


, , and is the distribution function of the convolution of independent gamma variables with shape parameter and scale parameter .


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
Table 3: Timing comparison (in microseconds) using different methods in to evaluate in the renewal process application when .

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
Table 4: Timing comparison (in milliseconds) using different methods in to evaluate in the renewal process application when .

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).


  • 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.