Mean absolute difference for the gamma distribution

Question

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)?

Answer 1

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))

Answer 2

Unfortunately, WolframAlpha will not calculate the relevant integral with my Standard account, any Pro user would be welcome to try.

That said, simulations strongly suggest that there is a simple sign error, and that the correct formula should be

$$ k\theta\big(2-4I_{\frac{1}{2}}(k+1,k)\big).$$

To see this, I first assumed without loss of generality that $\theta=1$ (it's a simple scaling factor that we can multiply in front of the result), then simulated the mean absolute difference $\text{MD}$ for $k=1, \dots, 100$, and finally regressed $\frac{\text{MD}}{k}$ on $I_{\frac{1}{2}}(k+1,k)$, obtaining parameters that were reasonably close to $2$ (for the intercept) and $-4$ (for $I_{\frac{1}{2}}(k+1,k)$).

Appending this theoretical formula to my simulations matches well.

    shape       sim     pbeta    theory
1       1  1.000461 0.2500000  1.000000
2       2  1.502236 0.3125000  1.500000
3       3  1.875163 0.3437500  1.875000
4       4  2.181705 0.3632813  2.187500
5       5  2.460114 0.3769531  2.460937
6       6  2.701749 0.3872070  2.707031
7       7  2.930308 0.3952637  2.932617
8       8  3.148350 0.4018097  3.142090
9       9  3.360636 0.4072647  3.338470
10     10  3.515429 0.4119015  3.523941
90     90 10.715282 0.4703059 10.689887
91     91 10.748322 0.4704690 10.749276
92     92 10.876108 0.4706295 10.808338
93     93 10.872729 0.4707874 10.867079
94     94 10.927393 0.4709428 10.925504
95     95 11.026546 0.4710957 10.983618
96     96 11.014793 0.4712463 11.041427
97     97 11.132948 0.4713945 11.098934
98     98 11.207437 0.4715404 11.156145
99     99 11.186310 0.4716842 11.213064
100   100 11.224114 0.4718258 11.269696

Simulation runs:

result <- data.frame(shape=1:100,sim=NA,pbeta=NA,theory=NA)
for ( ii in 1:nrow(result) ) {
    sims <- matrix(rgamma(2e5,shape=result$shape[ii],scale=1),ncol=2)
	result[ii,"sim"] <- mean(abs(sims[,1]-sims[,2]))
	result[ii,"pbeta"] <- pbeta(0.5,result$shape[ii]+1,result$shape[ii])
	result[ii,"theory"] <- result$shape[ii]*(-4*pbeta(0.5,result$shape[ii]+1,result$shape[ii])+2)
}
result[c(1:10,90:100),]
lm(sim/shape~pbeta,data=result)

Answer 3

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