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

• Hm. I played around a bit, and you seem to be right. The Wikipedia page unfortunately does not give a source for that table - I suspect the formula in the table is simply wrong. Of course, we can assume without loss of generality that $\theta=1$, but even so, WolframAlpha exceeds its standard computation time for the relevant integral. Do you have a Pro account? Commented Feb 11, 2021 at 7:24
• @BruceET: no, the beta is quite correct, at least based on that Wikipedia page, which claims that one should use the regularized incomplete beta function - which is just the CDF of the beta distribution. Commented Feb 11, 2021 at 7:27
• @StephanKolassa. OK, will erase my comment. But I can't get pbeta to work. Maybe the formula in the linked table is wrong, but my simulation matches the numerical value for the one instance I tried. Commented Feb 11, 2021 at 7:44

## 3 Answers

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

• I anticipate the Wikipedia page will soon be updated to reflect this result. ;-)
– whuber
Commented Feb 11, 2021 at 18:24

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)


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

• For experimenting with rgamma(), it's easiest to name the parameters and use rgamma(10^7, shape=4,scale=1). And I don't quite see how you find that the Wikipedia table gets it right. 4*1*(4*pbeta(0.5,4+1,4)-2) gives me -2.1875, which is the negative of your simulation result. (A relationship that does not hold for a shape of 2.) Commented Feb 11, 2021 at 7:52
• I meant that my simulation matched Wikipedia's numerical value. It's the Wikipedia formula I doubt. I got the same result from pbeta you got; just the negative of the correct answer. (Notice your R code using pbeta has a factor of $k$ that Wikipedia has, but OP does not have.) Commented Feb 11, 2021 at 8:16