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.

  • 2
    $\begingroup$ The best online reference is called Wikipedia. $\endgroup$
    – whuber
    Feb 3 at 0:22
  • 4
    $\begingroup$ Wikipedia is very handy reference, but - while it often gives an expression for moments - it usually doesn't give third absolute moments. If absolute moments are given at all, it's typically only the first. I can see some value in having these available, since, for example, they come up in the Berry-Esseen inequality. $\endgroup$
    – Glen_b
    Feb 3 at 0:27
  • 1
    $\begingroup$ @Glen_b Good point. This is a fraught exercise because there are all sorts of numerical issues involved with the automatic integration of an arbitrary distribution function. Here is a quick-and-dirty implementation which, if nothing else, will supply some rough reference values (and hints at an effective algorithm). f <- function(p = pnorm, ...) { mp <- integrate(\(x) 1 - p(x), 0, Inf, ...)$value; mm <- integrate(\(x) p(x), -Inf, 0, ...)$value; m <- mp + mm; integrate(\(x) 1 - (p(m + x^(1/3)) - p(m - x^(1/3))), 0, Inf, ...)$value } $\endgroup$
    – whuber
    Feb 3 at 0:34
  • $\begingroup$ @whuber I first went to Wikipedia as any good student does, but as the other user mentioned, Wikipedia doesn't really talk much about absolute moments. Thanks for your code, anyway, I know that the exercise is probably stupid, but as I said, I unfortunately have to do this. $\endgroup$ Feb 3 at 1:50
  • $\begingroup$ @whuber I have a question: how would I adapt your approach to the discrete case? I thought that the function I was using was correct, but it unfortunately wasn't. My problem is how to compute the series involved in the formula of the expected value. In the continuous case I know the "integrate" function that you also used, but in the discrete case I am at a loss. I mean, I thought that I could somehow truncate the series, but the index in my series is not an integer or something (as I foolishly believed initially), so what could I do? $\endgroup$ Feb 3 at 2:50

2 Answers 2


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.$


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:

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:

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.

  • $\begingroup$ thanks a lot for your answer! I think that I understand what you did, but I am still a bit confused about the use of > integrate for the discrete case. I mean, I tried this > f <- function(x){ if(x==1){ 0.7 } else{ 0 } } f <- Vectorize(f) integrate(f, 0, Inf, subdivisions=100L)$value and it returned 0. Shouldn't it have returned 0.7? Sorry for not managing to write the code in a block, I am still new here. $\endgroup$ Feb 3 at 13:48
  • $\begingroup$ I mean what I tried was to write the PMF of a Bernoulli with parameter $0.7$ and to integrate it to convince myself this works, but I think I made some mistake. $\endgroup$ Feb 3 at 13:56
  • 1
    $\begingroup$ I don't follow that question, because my code examples show how to do it for "properly defined" distribution functions. Of course if you pass that code a function p that does not implement a valid distribution function, then you will get erroneous results. The requirements are that it return predictable (not random) numbers; for very large negative values it returns 0; for large positive values it returns 1; and it never decreases when you increase its argument. (There are additional mathematical requirements that are irrelevant for computing.) $\endgroup$
    – whuber
    Feb 3 at 14:22
  • 1
    $\begingroup$ f <- \(x) pbinom(x, 1, 0.7) does it. I would guess you are reluctant to use some built-in statistical functions, but that leaves me wondering which functions you are willing to use. I will suppose only the most basic ones and propose something like H <- \(x) ifelse(x >= 0, 1, 0); f <- \(x) 0.3 * H(x) + 0.7 * H(x-1) because this shows how to generalize it to more complicated discrete distributions. Questions about how to code stuff are usually not considered on topic here unless they clearly require statistical explanations. $\endgroup$
    – whuber
    Feb 3 at 14:34
  • 2
    $\begingroup$ That's what we're for. But the labor is split up: programming questions are fielded on Stack Overflow (which has some very stats-knowledgeable users, too) while here, although we often offer programs as parts of solutions, we focus on the stats. $\endgroup$
    – whuber
    Feb 3 at 14:44

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)

  • $\begingroup$ thanks a lot! Is there any efficient way to compute these third absolute moments by hand? I haven't been taught how to do this in my classes and the only thing I could come up with was computing the PMF/CDF of $|X-\mu|^3$ in the generic case. $\endgroup$ Feb 3 at 1:52
  • $\begingroup$ Well, that's more work than you need to go to; LOTUS reduces the effort. That is, $E[g(X)] = \int_{-\infty}^{\infty} g(x) f(x) dx$, but particular cases can be reduced still further. I solved the normal case by using properties of the half-normal. I solved the exponential case by writing the integral in two parts (integral from $0$ to $\mu$ and $\mu$ to $\infty$ solved at $\mu=1$ and then generalized), then manipulated that ... ctd $\endgroup$
    – Glen_b
    Feb 3 at 10:10
  • $\begingroup$ ctd ... into an integral over the half line minus twice an integral over the first part; the first term is just a straight third central moment and the second is straightforward. The uniform reduces to evaluating a polynomial. $\endgroup$
    – Glen_b
    Feb 3 at 10:12
  • $\begingroup$ @justanamateur You can compute these values by integration: $m_3=\int |x-E(X)|^3 f(x)\,dx$, where $f(x)$ is the probability density. The non-centered moments can alternatively be computed by derivation of the "characteristic function" at $x=0$; the "characteristic function" is the Fourier transform of the probability density and you can find a list on wikipedia. $\endgroup$
    – cdalitz
    Feb 3 at 10:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.