# Difference of Gamma random variables

Given two independent random variables $X\sim \mathrm{Gamma}(\alpha_X,\beta_X)$ and $Y\sim \mathrm{Gamma}(\alpha_Y,\beta_Y)$, what is the distribution of the difference, i.e. $D=X-Y$?

If the result is not well-known, how would I go about deriving the result?

• I think may be relevant: stats.stackexchange.com/q/2035/7071 – Dimitriy V. Masterov Jan 23 '13 at 21:14
• Unfortunately not relevant, that post considers the weighted sum of Gamma random variables where the weights are strictly positive. In my case the weights would be +1 and -1 respectively. – FBC Jan 23 '13 at 21:17
• The Moschopoulos paper claims that the method can be extended to linear combinations, but you are right that the rescaling seems to be restricted to weights greater than 0. I stand corrected. – Dimitriy V. Masterov Jan 23 '13 at 21:41
• There's little hope of deriving anything simple or in closed form unless the two scale factors are the same. – whuber Jan 23 '13 at 21:41
• Just a small remark: for the special case of exponentially distributed rvs with the same parameter the result is Laplace (en.wikipedia.org/wiki/Laplace_distribution). – Ric Mar 28 '13 at 9:57

I will outline how the problem can be approached and state what I think the end result will be for the special case when the shape parameters are integers, but not fill in the details.

• First, note that $X-Y$ takes on values in $(-\infty,\infty)$ and so $f_{X-Y}(z)$ has support $(-\infty,\infty)$.

• Second, from the standard results that the density of the sum of two independent continuous random variables is the convolution of their densities, that is, $$f_{X+Y}(z) = \int_{-\infty}^\infty f_X(x)f_Y(z-x)\,\mathrm dx$$ and that the density of the random variable $-Y$ is $f_{-Y}(\alpha) = f_Y(-\alpha)$, deduce that $$f_{X-Y}(z) = f_{X+(-Y)}(z) = \int_{-\infty}^\infty f_X(x)f_{-Y}(z-x)\,\mathrm dx = \int_{-\infty}^\infty f_X(x)f_Y(x-z)\,\mathrm dx.$$

• Third, for non-negative random variables $X$ and $Y$, note that the above expression simplifies to $$f_{X-Y}(z) = \begin{cases} \int_0^\infty f_X(x)f_Y(x-z)\,\mathrm dx, & z < 0,\\ \int_{0}^\infty f_X(y+z)f_Y(y)\,\mathrm dy, & z > 0. \end{cases}$$

• Finally, using parametrization $\Gamma(s,\lambda)$ to mean a random variable with density $\lambda\frac{(\lambda x)^{s-1}}{\Gamma(s)}\exp(-\lambda x)\mathbf 1_{x>0}(x)$, and with $X \sim \Gamma(s,\lambda)$ and $Y \sim \Gamma(t,\mu)$ random variables, we have for $z > 0$ that \begin{align*}f_{X-Y}(z) &= \int_{0}^\infty \lambda\frac{(\lambda (y+z))^{s-1}}{\Gamma(s)}\exp(-\lambda (y+z)) \mu\frac{(\mu y)^{t-1}}{\Gamma(t)}\exp(-\mu y)\,\mathrm dy\\ &= \exp(-\lambda z) \int_0^\infty p(y,z)\exp(-(\lambda+\mu)y)\,\mathrm dy.\tag{1} \end{align*} Similarly, for $z < 0$, \begin{align*}f_{X-Y}(z) &= \int_{0}^\infty \lambda\frac{(\lambda x)^{s-1}}{\Gamma(s)}\exp(-\lambda x) \mu\frac{(\mu (x-z))^{t-1}}{\Gamma(t)}\exp(-\mu (x-z))\,\mathrm dx\\ &= \exp(\mu z) \int_0^\infty q(x,z)\exp(-(\lambda+\mu)x)\,\mathrm dx.\tag{2} \end{align*}

These integrals are not easy to evaluate but for the special case $s = t$, Gradshteyn and Ryzhik, Tables of Integrals, Series, and Products, Section 3.383, lists the value of $$\int_0^\infty x^{s-1}(x+\beta)^{s-1}\exp(-\nu x)\,\mathrm dx$$ in terms of polynomial, exponential and Bessel functions of $\beta$ and this can be used to write down explicit expressions for $f_{X-Y}(z)$.

From here on, we assume that $s$ and $t$ are integers so that $p(y,z)$ is a polynomial in $y$ and $z$ of degree $(s+t-2, s-1)$ and $q(x,z)$ is a polynomial in $x$ and $z$ of degree $(s+t-2,t-1)$.

• For $z > 0$, the integral $(1)$ is the sum of $s$ Gamma integrals with respect to $y$ with coefficients $1, z, z^2, \ldots z^{s-1}$. It follows that the density of $X-Y$ is proportional to a mixture density of $\Gamma(1,\lambda), \Gamma(2,\lambda), \cdots, \Gamma(s,\lambda)$ random variables for $z > 0$. Note that this result will hold even if $t$ is not an integer.

• Similarly, for $z < 0$, the density of $X-Y$ is proportional to a mixture density of $\Gamma(1,\mu), \Gamma(2,\mu), \cdots, \Gamma(t,\mu)$ random variables flipped over, that is, it will have terms such as $(\mu|z|)^{k-1}\exp(\mu z)$ instead of the usual $(\mu z)^{k-1}\exp(-\mu z)$. Also, this result will hold even if $s$ is not an integer.

• +1: Having looked at this problem before, I find this answer fascinating. – Neil G Mar 28 '13 at 5:57
• I'm going to accept this answer even though there appears to be no closed form solution. It's as close as it gets, thanks! – FBC Jan 9 '14 at 16:26
• I love the reasoning here, but I'm wondering if there is any measure where the second step breaks, I.e., $f_{-Y}(\alpha) ≠ f_{Y}(-\alpha)$? – mpacer Sep 29 '15 at 6:53
• @mpacer No, $f_{-Y}(\alpha) = f_{Y}(-\alpha)$ always holds. It is a general result that does not require any assumptions (normality, Gamma-eity, positive RV etc). For the special case of a positive random variable (that is, $P\{Y > 0\} = 1$), $-Y$ is a negative random variable that takes on values less than $0$ with probability $1$. – Dilip Sarwate Sep 29 '15 at 14:16
• @mpacer If $Y$ is a positive random variable with density $f_Y(\alpha)$, then it is not true that $f_Y(\alpha)$ is undefined for $\alpha<0$. In fact, $f_Y(\alpha)$ is defined as having value $0$ for $\alpha<0$. Thus, $f_{-Y}(\alpha)=f_Y(\alpha)=0$ for all positive numbers $\alpha$, and the density of $Y$ is the density of $Y$ "flipped over" with respect to the origin (or vertical axis if you prefer.) I am not "interpreting" the $-$ operator differently, it is you who is demanding an "appropriate" notion of $-$ that will support your idea that the domain of $f_Y$ is $\mathbb R^+$ only – Dilip Sarwate Sep 30 '15 at 20:48

To my knowledge the distribution of the difference of two independent gamma r.v.’s was first studied by Mathai in 1993. He derived a closed form solution. I will not reproduce his work here. Instead I will point you to the original source. The closed form solution can be found on page 241 as theorem 2.1 in his paper On non-central generalized Laplacianness of quadratic forms in normal variables.