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: [http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance]

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