10

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.

Community
  • 1
hengxin
  • 3,777
  • I deleted my answer as I had missed (a quite crucial) parenthesis! Will try to have a closer look later instead. – hejseb May 21 '14 at 09:53
  • @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). – hengxin May 21 '14 at 11:37
  • 3
    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. – Dilip Sarwate May 21 '14 at 13:14
  • @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... – hengxin May 21 '14 at 14:22
  • @DilipSarwate This problem actually arises from a Renewal process. I have done some calculation and come at the sum of two gamma random variables presented in this post. It is likely that I am taking a detour. I will think over my original problem. Thanks. – hengxin May 21 '14 at 14:35
  • 2
    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. – user64494 May 22 '14 at 02:33
  • @user64494 It is helpful. Thanks. If it is hard to get a closed formula, I will consider its numerical solution following your code. By the way, will the other assumption that integer numbers $\lambda \neq \mu$ lead to an alternate (maybe simpler) form? – hengxin May 22 '14 at 02:56

3 Answers3

4

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.

Suprasad Amari
  • 56
2

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.

Bombyx mori
  • 20,152
  • 1
    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. – hengxin May 22 '14 at 06:41
  • I updated the computation. I think you are right, the factorization now seems rather complicated. – Bombyx mori May 22 '14 at 15:29
  • I think I found a way out. – Bombyx mori May 22 '14 at 15:50
  • Some constants may still be off, I need to fix it. – Bombyx mori May 22 '14 at 15:59
  • 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). – hengxin May 23 '14 at 02:10
  • By my previous calculations, if $m,n$ are different then the result is messy. I do not see a clear way out. The same trick applies, but the result is not "impressive" anymore as I cannot factorize it. – Bombyx mori May 23 '14 at 02:22
  • Thanks. Your result $$exp(1) + k exp(1) \sim \Gamma(1,\frac{1}{2}+\frac{1}{2k})$$ is useful. Some other distributions than Gamma or Erlang, e.g., Hyper-exponential (weighted sum of exponentials) and Hypoexponential (general sum of exponentials) may be useful too. And, I believe, the answer (not likely in a neat formula) is hidden in site: Sums of Gamma Random Variables, although I failed to find it. – hengxin May 23 '14 at 02:45
  • A little question: In the calculation of the probability density function of $X_i + k Y_i$, how to get the constant $(\frac{1}{4} - \frac{1}{4k^2})$? I only obtain $\frac{1}{2k}$ which is $|J(X_i, Y_i)|^{-1}$. – hengxin May 24 '14 at 07:51
  • Yes, this is what I meant "some constants may still be off". By kernel method I should get $\frac{1}{4}-\frac{1}{4k^{2}}$, but direct calculation showed it to be $\frac{1}{2k}$. I was puzzled by it as well. I thought it is most likely to be a computation mistake on my side. – Bombyx mori May 24 '14 at 17:39
  • Does it mean that the intermediate conclusion that it is a Gamma random variable is debatable? – hengxin May 25 '14 at 01:22
  • The factorization should still hold, so it must be a Gamma distribution . But I am not entirely sure what was wrong with the constants. For example, if $k=1$, then the result I give obviously would not make sense. This means the whole result may be wrong. – Bombyx mori May 25 '14 at 01:38
  • For $k=1$, I am expecting $X_{i}+Y_{i}=\Gamma(1,1)+\Gamma(1,1)\sim \Gamma(2,1)$, and so is $X_{i}-Y_{i}$. But the result obviously mismatches. – Bombyx mori May 25 '14 at 01:41
  • To the specific problem, we may get stumped and frustrated. Generally, your efforts and attention are invaluable to me. – hengxin May 25 '14 at 01:50
  • I am quite embarrassed (as you can tell), for I will be the TA for probability class next semester. Allow me sometime to think about this carefully. – Bombyx mori May 25 '14 at 01:52
  • It would be OK to make mistakes on calculation (and really, who doesn't?). I am asking for help and you have helped me a lot. Hope that it will not upset you too much. – hengxin May 25 '14 at 01:59
  • Updated, but the result is still not good enough. It turns out the mistake with the constants is because of follows: $X-Y$'s range may no longer be in $[0,\infty)$. Therefore one has to integrate over $(-\infty, \infty)$ instead, and the integral would not converge. Here I used a scale transformation for $W$ instead. – Bombyx mori May 25 '14 at 18:14
  • Thanks a lot. I have posted a related problem at MO site and got some references which concern the sum of Gamma distributions: P. G. Moschopoulos (1985) The distribution of the sum of independent gamma random variables, Annals of the Institute of Statistical Mathematics, 37, 541-544 and A. M. Mathai (1982) Storage capacity of a dam with gamma type inputs, Annals of the Institute of Statistical Mathematics, 34, 591-597. They use moment generating functions too. As you suggested, I am considering its numerical solutions. Thanks again. – hengxin May 26 '14 at 01:41
  • I checked the first paper, which is in the author's webpage. It seems he used some ingenious algebraic identities to deduce the pdf from the mgf. The notation is a bit hard to read and there are no examples. But hopefully this solved the problem to some extent. – Bombyx mori May 26 '14 at 01:50
  • I have got another answer from MO site. However, it is beyond my current knowledge. Would you mind reviewing it at your convenience? Thanks. – hengxin May 27 '14 at 06:54
  • I read it briefly. The first part is clear, but I am slightly lost with the second part. I think since the probability $P(A<L<B)$ is assumed to be a constant, taking expectation of it should be itself. The first part comes from computation $P(L\le A)$ and $P(L\le B)$, then take the difference. The second part then take expectation of the two numbers with respect to $L$. Since $L\sim \exp(\lambda)$, its pdf is of the form $e^{-\lambda x}dx$(in some conventions). Therefore the expectation with respect to $L$ is a function of $\lambda$, namely $P(A\le L)$'s Laplace transform. – Bombyx mori May 27 '14 at 08:06
  • For the last step we use the fact that $E(XY)=E(X)E(Y)$, here we may let $X_{i}=\lambda (S_{i}+A_{i})$. Then we only need to calculate the term $E(e^{-p(S_i+A_i)})$, which he/she gave a hint on the final result. I do not have time to check the details, but I think he/she probably did the computation by hand and get this formula. – Bombyx mori May 27 '14 at 08:14
  • I have to say this is amazing work (if it is correct). – Bombyx mori May 27 '14 at 08:15
  • Yes, it is amazing. It avoids the messy calculation of convolution of two Gamma random variables. Thanks for your efforts. I will investigate it further. – hengxin May 27 '14 at 08:19
  • My only bit of misunderstanding is I do not see where he used the fact $S_{i},A_{i}$ are random variables instead of scalars. Perhaps this is done in the last step using the identity $E(X)=E(E(X|Y))$ or something similar. Otherwise computing it involves some kind of simple double integral. I assume he/she must have calculated it. The only way to verify it is to compute this once ourselves. – Bombyx mori May 27 '14 at 08:25
  • A comment comes from the author of the amazing solution at MO site. You can make comments at the MO site if you will. This Math.SE site is always "complaining about" the comments which go longer and longer. – hengxin May 28 '14 at 02:33
0

see eqn 3.383(1) of book tables of integrals Gradshteyn and Ryzhik, it reads

$$\int\limits_0^u x^{\nu-1}(u-x)^{\mu-1}e^{\beta x}\,\mathrm{d}x= B(\mu,\nu)\overline{u}^{\mu+\nu-1} {}_1F_1(\nu;\mu+\nu;\beta u)\\ [\operatorname{Re}\mu>0,\,\operatorname{Re}\nu>0].$$

Alexey Burdin
  • 5,299
amit patel
  • 17