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 V. Masterov Jan 23 '13 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 '13 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 V. Masterov Jan 23 '13 at 21:41
  • $\begingroup$ There's little hope of deriving anything simple or in closed form unless the two scale factors are the same. $\endgroup$ – whuber Jan 23 '13 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$ – 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.

| cite | improve this answer | |
  • 2
    $\begingroup$ +1: Having looked at this problem before, I find this answer fascinating. $\endgroup$ – Neil G Mar 28 '13 at 5:57
  • $\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 '14 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 '15 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$ – Dilip Sarwate Sep 29 '15 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$ – 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.

| cite | improve this answer | |

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.