Expected value and variance of the absolute value of difference of gamma distributed variables Consider two independent and identically distributed random variables $X$ and $Y$ such that
$$X, Y \sim \Gamma\left(\frac{1}{2}, 2^{-n+1}\right), $$
where $n$ is a constant.
I am trying to calculate
\begin{equation}
\text{E}[|X - Y|], ~~\text{and}~~ \text{Var}[|X - Y|].
\end{equation}
Is there any simple way of doing this without having to compute the double integrals? Or, is there at least any way to simplify the double integrals? The double integrals are too complicated to solve on Mathematica.
 A: Here's a sketch of a solution.
First, ignore the scale factor of $\sigma = 2^{-n+1}$ because it can be restored at the end.  Because the characteristic function of a Gamma$(1/2)$ variable like $X$ or $Y$ is
$$\psi_X(t) = \psi_Y(t) = (1 - it)^{-1/2},$$
the characteristic function of $X-Y$ (a version of the Fourier Transform of its density) is
$$\psi_{X-Y}(t) = \psi_X(t)\psi_Y(-t) = (1+t^2)^{-1/2}.$$
The density function of $X - Y$ therefore is recovered by applying the inverse Fourier transform, giving
$$f_{X-Y}(x) = \widehat{\psi_{X-Y}}(t) = \frac{1}{\pi}K_0(|x|),$$
a multiple of the Bessel $K_0$ function.
Because all events of the form $|X-Y|\le t$ for $t\ge 0$ correspond to events $-t \le X-Y \le t$ and the symmetry (and continuity) of $f_{X-Y}$ show the events $-t \le X-Y\le 0$ and $0 \le X-Y \le t$ must have equal probabilities, the density of $|X-Y|$ is twice the density of $X-Y$ for all positive values (and otherwise is zero, of course).  Thus
$$f_{|X-Y|}(x) = \frac{2}{\pi}K_0(|x|).$$
The mean therefore is
$$E\left[\,|X-Y|\,\right] = \int_0^\infty x f_{|X-Y|}(x)\,\mathrm{d}x =  \int_0^\infty x \frac{2}{\pi}K_0(|x|)\,\mathrm{d}x = \frac{2}{\pi}.$$
The raw second moment similarly is
$$E\left[\,|X-Y|^2\,\right] = \int_0^\infty x^2 \frac{2}{\pi}K_0(|x|)\,\mathrm{d}x = 1.$$
(We could have obtained this earlier from $\psi_{X-Y}$ because the second moment of $X-Y$ is that of $|X-Y|$ and
$$\psi_{X-Y}(t) = (1+t^2)^{-1/2} = 1 + \binom{-1/2}{1}(t^2)^1 + \cdots = 1 - \frac{t^2}{2} + \cdots$$
and the second moment is, as always, $2!$ times the coefficient of $(it)^2.$)
Therefore
$$\operatorname{Var}(|X-Y|) = 1 - \left(\frac{2}{\pi}\right)^2 \approx 0.5947153.$$
Finally, when $X$ and $Y$ are multiplied by the same constant $\sigma,$ $|X-Y|$ is also multiplied by $\sigma,$ whence its mean is $(2/\pi)\sigma$ and its variance is $\left(1-(2/\pi)^2\right)\sigma^2.$
Results of a simulation of $|X-Y|$ compared to the theoretical density (shown in red):

This R code conducted the simulation and produced the graphic.
n <- 1e6
z <- abs(rgamma(n, 1/2) - rgamma(n, 1/2))

hist(z, xlim=c(0,5), breaks=200, freq=FALSE, col="Gray", main="Realizations of |X-Y|",
     xlab=expression(abs(X-Y)))
curve(besselK(x, 0) * 2/pi, add=TRUE, col="Red", lwd=2, n=1501)

(cbind(Mean=c(Simulated=mean(z), Theoretical=2/pi),
       Variance=c(var(z), 1 - (2/pi)^2)))

