Mean absolute difference for the gamma distribution A wikipedia entry states that the mean absolute difference for the $\Gamma(k,\theta)$ distribution is $k\theta(4I_{0.5}(k+1,k)-2)$ where $I_z(x,y)$ is the regularized incomplete beta function, equal to the CDF of a beta distribution. However, R tells me that 4*pbeta(0.5,k+1,k)-2 is negative, e.g., for $k=2$ it is $-.75$, which is impossible for the mean absolute difference. Does anyone have the correct formula (or see something I'm doing wrong)?
 A: There's a simple typo.  To prove it, let's derive the result.
The density function for a $\Gamma(k)$ variable is
$$f_k(x) = \frac{1}{\Gamma{(k)}} x^{k-1}e^{-x}.$$
($\theta$ is a scale parameter and since $|X-Y|$ is directly proportional to the scale factor, we may set $\theta=1$ and at the end multiply the result by $\theta.$)
Thus the density function for two independent such variables $(X,Y)$ is
$$f_{k,k}(x,y) = f_k(x)f_k(y).$$
Evaluate $E[|X-Y|]$ by changing variables to
$$(t,u) = \left(\frac{x}{x+y}, x+y\right)$$
so that
$$x = tu,\ y = (1-t)u,\ |x-y| = |2t-1|u,\ \mathrm{d}x\mathrm{d}y = u\,\mathrm{d}t\mathrm{d}u$$ and the integration is over the region $0\lt u\lt\infty,$ $0\lt t\lt 1,$ whence
$$\begin{aligned}
E[|X-Y|] &= \frac{1}{\Gamma(k)^2} \iint_0^\infty |x-y|x^{k-1}y^{k-1}e^{-(x+y)}\,\mathrm{d}x\mathrm{d}y\\
&=\frac{1}{\Gamma{(k)}^2} \int_0^\infty u^{2k}e^{-u}\mathrm{d}u \int_0^1 |2t-1|t^{k-1}(1-t)^{k-1}\,\mathrm{d}t.
\end{aligned}\tag{*}$$
The $u$ integral equals $\Gamma(2k+1)$ (by definition) while the $t$ integral can evaluated by observing its integrand is symmetric around $t=1/2,$ giving
$$\begin{aligned}
\int_0^1 |2t-1|t^{k-1}(1-t)^{k-1}\,\mathrm{d}t &= 2\int_0^{1/2} (1-2t)t^{k-1}(1-t)^{k-1}\,\mathrm{d}t\\
&= 2\int_0^{1/2} t^{k-1}(1-t)^{k-1}\,\mathrm{d}t - 4\int_0^{1/2} t^k(1-t)^{k-1}\,\mathrm{d}t.
\end{aligned}$$
The last integral is an unnormalized cumulative Beta functions with parameters $(k+1,k).$  The integral from which it is subtracted is symmetric around $1/2,$ whence its value is $B(k,k)/2.$  Writing $F(x,a,b)$ for the Beta$(a,b)$ CDF allows us to express the expectation $(*)$ as

$$\begin{aligned}
E[|X-Y|] &= \frac{\Gamma(2k+1)}{\Gamma(k)^2}\left(\frac{B(k,k)}{2} - B(k+1,k)F(1/2,k+1,k)\right) \\
&= 4k\left(\frac{1}{2}-B(1/2,k+1,k)\right) \\
&= 4k\left(B(1/2,k,k+1)-\frac{1}{2}\right),
\end{aligned}$$

where the last equality follows from the symmetry between any $\text{Beta}(a,b)$ and $\text{Beta}(b,a)$ distribution.
This result is the negative of that in the Wikipedia entry.  Evidently either the Wikipedia entry subtracted the terms in the wrong order or it swapped the parameters $k$ and $k+1.$
One nice thing about integral formulas is we can often check them with numerical integration, as in the right hand plot of this figure that compares numerical to analytical results for various values of $k$ between $0$ and $4.$

The R code to perform these checks is reasonably short.
par(mfrow=c(1,2))
#
# Analytical formula.
#
e <- function(k) 4 * k * (pbeta(1/2, k, k+1) - 1/2)
curve(e(x), 0, 5, n=1001, lwd=2, ylab="Value",
       main=expression(4*k*group("(", Beta(1/2,k,k+1)-1/2, ")")))
#
# Numerical comparison.
#
f <- Vectorize(function(k, ...) {
  p <- function(t) dgamma(t, k)
  g <- Vectorize(function(y) integrate(function(x) p(y)*p(x)*abs(x-y), 0, Inf, ...)$value)
  integrate(g, 0, Inf, ...)$value
}, "k")

k <- seq(0, 2, length.out=21)^2
plot(f(k), e(k), type="b", cex=0.7, asp=1,
     xlab="Numerical", ylab="Analytical", main="Comparison")
abline(0:1, col="Gray")
par(mfrow=c(1,1))

A: Comment:
The table in your Wikipedia link seems to be correct for
the distribution $\mathsf{Gamma}(4, 1).$ It has $2.1875$ and
by simulation in R, I get $2.188.$ With 10 million iterations,
I would expect about 3 place accuracy.
set.seed(2021)
x = rgamma(10^7, 4, 1)
y = rgamma(10^7, 4, 1)
mean(abs(x-y))
[1] 2.188033

If you're going to try using R, for further exploration,
you should note that the second parameter in R is the rate $\lambda$,
not the scale $\theta = 1/\lambda.$ (Makes no difference for $\lambda = \theta = 1.)$
Addendum: For shape 2, my simulation also matches Wikipedia's
numerical value $1.5.$
set.seed(2021)
x = rgamma(10^7, 2, 1)
y = rgamma(10^7, 2, 1)
mean(abs(x-y))
[1] 1.500328

Arbitrarily modifying Wikipedia's formula to get positive rather than
negative values from rbeta, I get the same numerical results
for shape $k$ and rate $1$ as in the Wikipedia formula as below.
Not saying this is now the correct formula; just that it gives
the same numerical results as in the Wikipedia formula for
the four available cases (two of which also match simulated values).
k = 1:4
k*(2 - 4*pbeta(.5, k+1, k))
[1] 1.0000 1.5000 1.8750 2.1875

