Absolute third central moment for standard distributions reference

Question

I have to write an R function that computes the absolute third central moment (i.e. $\mathbb{E}[|X-\mathbb{E}[X]|^3]$) in the cases that you are given the name of the distribution or the PMF/CDF. I know how to do this, but what I want is a reference for the absolute third central moment for distributions such as Bernoulli, Binomial, Geometric, Exponential, Gamma etc. (the pretty standard ones) so that I can easily check the output of my function. I don't know of any predefined R functions that do this either.

Answer 1

It sounds like a little analysis might help you.

The absolute moment of order $k$ (I will consider only $k\gt 0$ for convenience) around a central value $a$ for a distribution function $F$ is defined to be

$$\nu_k(F;a) = \int |x - a|^k\,\mathrm{d}F(x)$$

in the sense of Lebesgue-Stieltjes or Riemann-Stieltjes integration: this applies to continuous and discrete distributions.

If we let $X$ be a random variable with distribution $F,$ this is a standard formula for an expectation

$$\nu_k(F;a) = E\left[|X-a|^k\right].$$

But the distribution function of $Y = |X-a|^k,$ which has non-negative support, is readily computed from its definition for all $y\ge 0$ as

$$\begin{aligned} F_Y(y) &= \Pr(Y\le y) = \Pr(|X-a|^k\le y) = \Pr(a - y^{1/k}\le X \le a + y^{1/k})\\ &= \Pr(X\le a + y^{1/k}) - \Pr(X \lt a - y^{1/k})\\ &= F(a + y^{1/k}) - F^{-}(a - y^{1/k}) \end{aligned}$$

where I have written

$$F^{-}(x) = \lim_{x\to 0^{-}} F(x) = \Pr(X \lt x).$$

Integration by parts shows that for any random variable $Y$ whose expectation exists (even if it's infinite),

$$E[Y] = \int_0^\infty (1 - F_Y(y))\,\mathrm d y - \int_{-\infty}^0 F_Y(y)\,\mathrm d y.\tag{*}$$

When $Y$ has non-negative support, as here, the second term vanishes.

Putting these facts together gives

$$\nu_k(F;a) = E[Y] = \int_0^\infty (1 - F_Y(y))\,\mathrm d y = \int_0^\infty 1 - F(a + y^{1/k}) + F^{-}(a - y^{1/k})\,\mathrm d y.$$

Now for any $F$ you would want to consider, you may ignore up to a countable infinity of values, which means you don't have to worry about the distinction between $F$ and $F^{-}$ in the integral. The formula to consider implementing for any absolute moment around the mean is

$$\nu_k(F) = \int_0^\infty 1 - F(\mu(F) + y^{1/k}) + F(\mu(F) - y^{1/k})\,\mathrm d y$$

where, as a preliminary, you have previously computed the expectation $\mu(F) = E[X]$ by applying formula $(*)$ to $F.$

Remarks

I warmly recommend implementing this more general function for arbitrary $k$ rather than focusing on $k=3,$ if only because it will permit you to test it with more readily-verified values like $k=1$ and $k=2.$

---

You can test your implementation with simulation. Here are examples in R. The integration is quick and dirty: it does not check the results and has only the barest provision to control the integration (via the number of subdivisions s). The function mu computes ordinary or absolute moments about a central value and is used for all the integration.

mu <- function(k = 1, p = pnorm, a = 0, abs.value = TRUE, s = 100L, ...) {
  sgn <- if(isTRUE(abs.value)) 1 else -1
  integrate(\(x) 1 - p(a + x^(1/k),...) + sgn * p(a - x^(1/k), ...), 
            0, Inf, subdivisions = s)$value
}
mu.p <- function(k = 1, p = pnorm, s = 100L, ...) {
  m <- mu(1, p, 0, FALSE, s, ...) # Mean
  mu(k, p, m, TRUE, s, ...)       # Absolute moment about the mean
}

The optional arguments ... to mu and mu.p are any parameters you might want to pass to the distribution function p. For example, here is the third absolute central moment of a Poisson$(2)$ distribution compared to a Monte-Carlo estimate:

set.seed(17)
x <- rpois(1e6, 2)
print(c(Calculated = mu.p(3, ppois, lambda = 2, s = 1e4), 
        `Monte-Carlo` = mean(abs(x - mean(x))^3)))

 Calculated Monte-Carlo 
   4.706693    4.719109

The difference is attributable to random fluctuations among the million values used in the M-C estimate.

That was a Poisson calculation. Let's try a Negative Binomial distribution (and, for fun, change $k$ from $3$ to $2$):

mu.p(2, pnbinom, size = 4, prob = 1/3)

Error in integrate(function(x) 1 - (p(a + x^(1/k), ...) - p(a - x^(1/k),  : 
  maximum number of subdivisions reached

That illustrates what you might be up against with numerical integration of arbitrary CDFs. The jumps give the integrator the willies. Let's increase the subdivisions with the s argument:

set.seed(17)
x <- rnbinom(1e6, size = 4, prob = 1/3)
print(c(Calculated = mu.p(2, pnbinom, size = 4, prob = 1/3, s = 1e3), 
        `Monte-Carlo` =  mean(abs(x - mean(x))^2)))

 Calculated Monte-Carlo 
   23.99989    24.11780

Again the agreement is fine.

---

The bigger challenge lies in implementing a numerical (or symbolic) integrator that handles a wide variety of distribution functions $F.$ Doing that well and robustly is a huge enterprise. So, if you can anticipate what $F$ is likely to be, you can optimize your code for that class of distributions.

Answer 2

Here's a few values to get started with since it might occasionally be handy to have some third absolute moments to refer to. I've made this answer community wiki -- anyone should feel free to edit to include more if they wish. I will try to add more when I get some time.

In each case $\sigma^2$ is the variance of the distribution in question (this way the specifics of the parameterization doesn't matter).

Normal: $\sqrt{\frac{\pi}{8}}\,\sigma^3$ $\quad$ ($\approx 1.596$ for standard normal)

Exponential: $(\frac{12}{e}-2)\,\sigma^3$ $\quad$ ($\approx 2.415$ for standard exponential)

Uniform (continuous): $\frac{12^{3/2}}{32}\,\sigma^3$ $\quad$ ($\frac{1}{32}=0.03125$ for standard uniform)

(computed algebraically, checked by large simulation)