Single-pass algorithm for kurtosis

Question

Here is a simple test I've run on MATLAB to check the validity of a single pass (online) algorithm for computing $3$rd moment and $4$th moment.

randn('state',0);

num2 = 0;
num1 = 0;
delta = 0;  
M1 = 0;
M2 = 0;
M3 = 0;
M4 = 0;

Xvec = zeros(1, 100000);

for j = 1:length(Xvec)

    X = 3*randn(1);

    % Single pass algorithm   
    num2  = num1;
    num1  = num1 + 1;
    delta = X - M1;
    delta_n  = delta / num1;
    delta_n2 = delta_n * delta_n;
    term1 = delta *delta_n * num2;
    M1 = M1 + delta_n;
    M4 = M4 + term1 * delta_n2 * (num1 * num1 - 3.0*num1 + 3.0) + ...
    6.0 * delta_n2 * M2 - 4.0 * delta_n * M3;
    M3 = M3 + term1 * delta_n * (num1 - 2.0) - 3.0 * delta_n * M2;
    M2 = M2 + term1;

    Xvec(j) = X;
end
% Quantities obtained from Single Pass
avg = M1;
variance = M2 / (num1 - 1.0);
kurtosis = (num1 * M4)/(M2 * M2) - 3.0;
skewness = sqrt(num1) * M3 / (M2^1.5);

% Reference Quantities 
avg1     = sum /length(Xvec);
moment_2 = moment(Xvec, 2);
moment_3 = moment(Xvec, 3);
moment_4 = moment(Xvec, 4);

The online algorithm provides a correct mean (avg = avg1) and variance (variance = moment_2). However, the values for kurtosis and skewness obtained from the online algorithm are way off from the actual $3$rd and $4$th moments.

What might be going wrong?

Source

Answer 1

From what I can tell, you've changed the notation from the Wikipedia article in a way that is quite confusing. You've replaced $m$ with $M_1$ which is not consistent with $M_2$, $M_3$ and $M_4$. $M_k$ is defined as $\sum_i (x_i - \bar{x})^k$ whereas $m$ is the running mean.

In case this is where the source of the confusion lies: $M_k$ is not the $k$-th central moment. $\frac{M_k}{n}$ is.

Another possible source of confusion: from your question you seem to imply that the kurtosis and skewness should equal the third and fourth central moments which they do not - as you can see from the very calculations you are performing to calculate them:

kurtosis = (num1 * M4)/(M2 * M2) - 3.0;
skewness = sqrt(num1) * M3 / (M2^1.5);

I do not know exactly how the moment() function is implemented in Matlab but assuming it calculates central moments then you should see that your value for kurtosis is equal to moment_4/moment_2^2 - 3.0