MOMENT() don't give the rigth numbers for Variance, Skewness & Kurtosis in some degenerate cases.

[alaingdl]

When all values are equal, Variance should be equal to 0.0 and Skewness and Kurtosis must be equal to NaN

GDL> moment(replicate(11., 1000))
       11.000000       0.0000000             NaN             NaN
GDL> moment(replicate(2.1, 1000))
       2.0999856   2.0484114E-10      0.99848902      -2.0020103
GDL> moment(replicate(1.1d, 1000))
       1.1000000000000087   7.5066155958530961E-29     -0.99850037506254774      -2.0019989999999828

IDL> moment(replicate(11., 1000))
       11.000000       0.0000000             NaN             NaN
IDL> moment(replicate(2.1, 1000))
       2.0999856       0.0000000             NaN             NaN
IDL> moment(replicate(1.1d, 1000))
       1.1000000000000087       0.0000000000000000                      NaN                      NaN

We have a similar trouble with Variance and STDDEV() but I will create another bug report.
And I realize we are missing some systematic tests on that in the testsuite !

Just to illustrate the (classical) subtlety of the problem :

GDL> moment(replicate(2.1d, 10000))
       2.0999999999998802   1.4378427840436983E-26      0.99985000375017119      -2.0001999899997842
GDL> moment(replicate(2.d, 10000))
       2.0000000000000000       0.0000000000000000                      NaN                      NaN

[GillesDuvert]

@alaingdl I do not believe this is necessary (and idem for STDDEV). people have to live with that, it is normal floating-point precision loss.

For the record, I note that, using two different ways of computing the mean, a simple (but expensive in terms of add/mults) running mean and the (unknown) IDL algorithm, one sees transition points (depending of the initial value), indicative of an a posteriori correction:

IDL> n=1086 & a=replicate(1.2,n) & b=a[0] & for i=1,n-1 do b=(temporary(b)*i+a[i])/(i+1) & print,b,moment(a)
      1.20000      1.20000      0.00000          NaN          NaN
IDL> n=1087 & a=replicate(1.2,n) & b=a[0] & for i=1,n-1 do b=(temporary(b)*i+a[i])/(i+1) & print,b,moment(a) 
      1.19999      1.20000 -8.17844e-19         -NaN         -NaN
% Program caused arithmetic error: Floating illegal operand
IDL> n=1088 & a=replicate(1.2,n) & b=a[0] & for i=1,n-1 do b=(temporary(b)*i+a[i])/(i+1) & print,b,moment(a)
      1.19999      1.20000      0.00000          NaN          NaN

(...)

IDL> n=1117 & a=replicate(1.2,n) & b=a[0] & for i=1,n-1 do b=(temporary(b)*i+a[i])/(i+1) & print,b,moment(a)
      1.19999      1.20000      0.00000          NaN          NaN
IDL> n=1118 & a=replicate(1.2,n) & b=a[0] & for i=1,n-1 do b=(temporary(b)*i+a[i])/(i+1) & print,b,moment(a)
      1.19999      1.20000 -1.59029e-18         -NaN         -NaN
% Program caused arithmetic error: Floating illegal operand
IDL> n=1119 & a=replicate(1.2,n) & b=a[0] & for i=1,n-1 do b=(temporary(b)*i+a[i])/(i+1) & print,b,moment(a)
      1.19999      1.20000      0.00000          NaN          NaN

(...)

IDL> n=2003 & a=replicate(1.2,n) & b=a[0] & for i=1,n-1 do b=(temporary(b)*i+a[i])/(i+1) & print,b,moment(a)
      1.19998      1.19998      0.00000          NaN          NaN
IDL> n=2004 & a=replicate(1.2,n) & b=a[0] & for i=1,n-1 do b=(temporary(b)*i+a[i])/(i+1) & print,b,moment(a)
      1.19998      1.19998 -1.14084e-14         -NaN         -NaN
% Program caused arithmetic error: Floating illegal operand

Apparently there is a Floating illegal operand warning when variance is not returned as 0, then no more Floating illegal operand, and variance is 0 again, until another transition etc. Each time when skewness passes from NaN to -NaN. May be a clue.