6
$\begingroup$

The gamma distribution with parameters $m > 0$ and $\lambda > 0$ (denoted $\Gamma(m, \lambda)$) has density function $$f(x) = \frac{\lambda e^{-\lambda x} (\lambda x)^{m - 1}}{\Gamma(m)}, x > 0$$ Given two independent gamma random variables $X = \Gamma(m,\lambda)$ and $Y = \Gamma(n, \mu)$ with integer numbers $m \neq n, \lambda \neq \mu$, what is the density function of their sum $X + Y = \Gamma(m,\lambda) + \Gamma(n,\mu)$?

Notice that both $X$ and $Y$ are also Erlang distribution since $m,n$ are positive integers.

My attempt

First, I searched for well-known results about gamma distribution and I got:

(1) If $\lambda = \mu$, the sum random variable is a Gamma distribution $\sim \Gamma(m+n, \lambda)$ (See Math.SE).

(2) $\Gamma(m, \lambda)$ (or $\Gamma(n, \mu)$) is the sum of $m$ (or $n$) independent exponential random variables each having rate $\lambda$ (or $\mu$). The hypoexponential distribution is related to the sum of independent exponential random variables. However, it require all the rates distinct.

(3) This site is devoted to the problem of sums of gamma random variables. In section 3.1, it claims that if $m$ and $n$ are integer numbers (which is my case), the density function can be expressed in terms of elementary functions (proved in section 3.4). The answer is likely buried under a haystack of formulas (however, I failed to find it; you are recommended to have a try).

Then, I try to calculate it:

$$f_{X+Y}(a) = \int_{0}^{a} f_{X}(a-y) f_{Y}(y) dy \\ = \int_{0}^{a} \frac{\lambda e^{-\lambda (a-y)} (\lambda (a-y))^{m-1}}{\Gamma(m)} \frac{\mu e^{-\mu y} (\mu y)^{n-1}}{\Gamma(n)} dy \\ = e^{-\lambda a} \frac{\lambda^m \mu^n}{\Gamma(m) \Gamma(n)} \int_{0}^{a} e^{(\lambda - \mu) y} (a-y)^{m-1} y^{n-1} dy$$

Here, I am stuck with the integral and gain nothing ... Therefore,

  1. How to compute the density function of $\Gamma(m,\lambda) + \Gamma(n,\mu)$ with integer numbers $m \neq n, \lambda \neq \mu$?

  2. Added: The answers assuming $m = n$ ($\lambda \neq \mu$) are also appreciated.

$\endgroup$
  • $\begingroup$ I deleted my answer as I had missed (a quite crucial) parenthesis! Will try to have a closer look later instead. $\endgroup$ – hejseb May 21 '14 at 9:53
  • $\begingroup$ @hejseb Thank you all the same. And you are recommended to refer to the material: Sums of Gamma Random Variables mentioned in the post if you want to come back. The answer is likely buried under a haystack of formulas (however, I failed to find it). $\endgroup$ – hengxin May 21 '14 at 11:37
  • 2
    $\begingroup$ You are almost there. Since $m$ and $n$ are integers, expand $(a-y)^{m-1}y^{n-1}$ via the binomial theorem into a polynomial in $y$. Then you are left with a sum of integrals of the form $\displaystyle \int_0^a y^ie^{\nu y}\, dy$ each of which can be integrated by parts. $\endgroup$ – Dilip Sarwate May 21 '14 at 13:14
  • $\begingroup$ @DilipSarwate Following your instruction, I get the integrals of the form $\int_{0}^{a} e^{(\lambda - \mu) y} y^{n+x-1} dy$. Here $x$ is related to the general term of binomial extension of $(a-y)^{m-1}$ . Using Mathematica, I get the Gamma distribution (i.e., $\Gamma(n+x, \mu - \lambda)$) back. In addition, I find it hard to combine the binomial terms together after computing the integrals. Stuck again... $\endgroup$ – hengxin May 21 '14 at 14:22
  • 1
    $\begingroup$ The Maple code $$with(Statistics); X := RandomVariable(GammaDistribution(m, lambda)): Y := RandomVariable(GammaDistribution(n, mu)): PDF(X+Y, t)\, assuming\, m::posint, n::posint, m <> n; $$ produces the answer in terms of LaguerreL. $\endgroup$ – user64494 May 22 '14 at 2:33
3
$\begingroup$

A closed form expression is provided in the following paper.

SV Amari, RB Misra, Closed-form expressions for distribution of sum of exponential random variables, IEEE Transactions on Reliability, 46 (4), 519-522.

$\endgroup$
2
$\begingroup$

Update:

We summarize the development in follows:

Step 1: We simplify the case to be $l(\Gamma(m, 1)+k\Gamma(n,1))$ by choosing appropriate $k,l$ for a scale transformation.

Step 2. We want to calculate $$ \Gamma(m,1)+k\Gamma(n,1)=\sum^{m}_{i=1}X_{i}+k\sum^{n}_{i=1}Y_{j} $$

Step 3. For $m=n$ case, we only need to calculate $X+kY$, where $X,Y\sim \Gamma(1,1)$. In the case of $k=1$, we let $$ Z=X+Y,W=\frac{X}{X+Y}, X=ZW, Y=Z-ZW $$ The Jacobian is $Z$. Therefore we have $$ f_{Z,W}(z,w)=ze^{-zw}e^{zw-z}dzdw=ze^{-z}dzdw $$ and $$ Z\sim \Gamma(2,1) $$ as desired. In the general case we have $$ Z=X+kY, W=\frac{X}{X+kY}, X=ZW, Y=\frac{1}{k}(Z-ZW) $$ The Jacobian is $\frac{1}{k}Z$. We thus have $$ f_{Z,W}(z,w)=\frac{1}{k}ze^{-zw}e^{\frac{1}{k}(zw-z)}=\frac{1}{k}ze^{-\frac{k-1}{k}zw-\frac{1}{k}z} $$ and I do not have a good way to factorize it.

A reason this technique might not work in general is the moment generating function does not change when we use the scale transformation, and for different $\beta$ the moment generating function is different. Thus the problem may be better to be attacked numerically.

$\endgroup$
  • $\begingroup$ As you suggest, the computation of the density function of $X + Y$ can be reduced to that of $X' + kY''$. However, why is $X', Y'' \sim \text{exp}(1)$? In my calculation, it is of form $\sim \Gamma(n, 1) = \frac{e^{-x} x^{n-1}}{\Gamma(n)} = \frac{e^{-x} x^{n-1}}{(n-1)!}$ since $n$ is a positive integer. $\endgroup$ – hengxin May 22 '14 at 6:41
  • $\begingroup$ I updated the computation. I think you are right, the factorization now seems rather complicated. $\endgroup$ – Bombyx mori May 22 '14 at 15:29
  • $\begingroup$ I think I found a way out. $\endgroup$ – Bombyx mori May 22 '14 at 15:50
  • $\begingroup$ Some constants may still be off, I need to fix it. $\endgroup$ – Bombyx mori May 22 '14 at 15:59
  • $\begingroup$ Thanks for your efforts. The factorization trick is rather impressive. However, it seems that you have missed the requirement that $m \neq n$. Following your instruction, I get: $$X + Y = \lambda \Gamma(m,1) + \mu \Gamma(n,1) = \lambda (\Gamma(m,1) + \frac{\mu}{\lambda} \Gamma(n,1)) = \lambda (\sum_{i=1}^{m} \text{exp}(1) + \frac{\mu}{\lambda} \sum_{i=1}^{n} \text{exp}(1))$$. Unfortunately, we cannot factor $\sum_{i=1}^{i=n}$ out and focus on $$\text{exp}(1) + \frac{\mu}{\lambda} \text{exp}(1)$$. What do you think of it? (Now it has been reduced to the weighted sum of $\text{exp}(1)$s). $\endgroup$ – hengxin May 23 '14 at 2:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.