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?

  • $\begingroup$ I think may be relevant: stats.stackexchange.com/q/2035/7071 $\endgroup$
    – dimitriy
    Jan 23, 2013 at 21:14
  • 4
    $\begingroup$ 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. $\endgroup$
    – FBC
    Jan 23, 2013 at 21:17
  • $\begingroup$ 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. $\endgroup$
    – dimitriy
    Jan 23, 2013 at 21:41
  • 3
    $\begingroup$ 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). $\endgroup$
    – Richi W
    Mar 28, 2013 at 9:57
  • 1
    $\begingroup$ maybe this helps: math.kit.edu/stoch/~klar/seite/veroeffentlichungen/media/… $\endgroup$
    – bonanza
    Jun 12, 2015 at 8:02

3 Answers 3


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.

  • 2
    $\begingroup$ +1: Having looked at this problem before, I find this answer fascinating. $\endgroup$
    – Neil G
    Mar 28, 2013 at 5:57
  • 1
    $\begingroup$ 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! $\endgroup$
    – FBC
    Jan 9, 2014 at 16:26
  • $\begingroup$ 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)$? $\endgroup$
    – mpacer
    Sep 29, 2015 at 6:53
  • $\begingroup$ @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$. $\endgroup$ Sep 29, 2015 at 14:16
  • 1
    $\begingroup$ @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 $\endgroup$ Sep 30, 2015 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.


the difference of two independent or correlated Gamma random variables are special cases of McKay distribution. The exact and complete answer can be find in:

Sum and difference of two squared correlated Nakagami variates in connection with the McKay distribution Holm, H., Alouini, M.-S.

IEEE Transactions on Communications, 2004 Vol. 52; Iss. 8


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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