Sum of two gamma/Erlang random variables $\Gamma(m,\lambda)$ and $\Gamma(n, \mu)$ with integer numbers $m \neq n, \lambda \neq \mu$

Question

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.

Answer 1

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.

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

Answer 3

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].$$