Absolute difference of two gamma distribution

Question

Given probability density functions of two Gamma distributions $X$ and $Y$:

$f_1 (x) \propto x^{\alpha_1 -1} e^{-\beta_1 x}$

and

$f_2 (y) \propto y^{\alpha_2 -1} e^{-\beta_2 y}$.

How to compute the pdf of $|X-Y|$?

Answer 1

I will not attempt to solve this in closed form, that looks difficult, so I go directly for a numerical solution modeled on Difference of two i.i.d. lognormal random variables. Let $X$ and $Y$ be independent random variates with pdf, cdf $f,F$ and $g,G$, respectively. With $D=X-Y$ we find that the density of the absolute difference $| D |$ is $$ f_{|D|}(t) = \int_{\text{Range}(Y)} \left( f(t+y)+f(-t+y) \right)g(y)\; dy $$ (where we understand that the density $f$ is zero outside the range of $X$).

and implementing that in R gives:

    make_dabsDIFF  <-  function(f, g, rangeY=c(0, +Inf) ) {
        function(t) {
            res  <-  t
            for (tt in seq(along=t)) {
                res[tt]  <-  integrate(Vectorize(function(y)
                   if ((y-t[tt]) >= 0.0) { (f(y+t[tt]) + f(y - 
                       t[tt]))*g(y)} else {f(y+t[tt])*g(y)}), 
                       lower=rangeY[1], upper=rangeY[2])$value   
            }
            return(res)    
            }
    }

A symmetric case:

     dabsDIFF <- make_dabsDIFF(function(x) dgamma(x, 1, 1), 
                       function(x) dgamma(x, 1, 1))

and a plot:

Then a case with two gamma distributions with different parameters:

    dabsDIFF2 <- make_dabsDIFF(function(x) dgamma(x, 1, 1), 
                   function(x) dgamma(x, 0.5, 2.0))