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 Answers
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.
-
$\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 likeH <- \(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_bFeb 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_bFeb 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$– cdalitzFeb 3 at 10:12
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$