Source code for moment.pro:

;$Id: moment.pro,v 1.22 2001/01/15 22:28:08 scottm Exp $
;
; Copyright (c) 1994-2001, Research Systems, Inc.  All rights reserved.
;       Unauthorized reproduction prohibited.
;+
; NAME:
;       MOMENT
;
; PURPOSE:
;       This function computes the mean, variance, skewness and kurtosis
;       of an N-element vector. IF the vector contains N identical elements, 
;       the mean and variance are computed; the skewness and kurtosis are 
;       not defined and are returned as IEEE NaN (Not a Number). Optionally, 
;       the mean absolute deviation and standard deviation may also be 
;       computed. The returned result is a 4-element vector containing the
;       mean, variance, skewness and kurtosis of the input vector X.
;
; CATEGORY:
;       Statistics.
;
; CALLING SEQUENCE:
;       Result = Moment(X)
;
; INPUTS:
;       X:      An N-element vector of type integer, float or double.
;
; KEYWORD PARAMETERS:
;       DOUBLE: IF set to a non-zero value, computations are done in
;               double precision arithmetic.
;
;       MDEV:   Use this keyword to specify a named variable which returns
;               the mean absolute deviation of X.
;
;       SDEV:   Use this keyword to specify a named variable which returns
;               the standard deviation of X.
;
;       MAXMOMENT:
;               Use this keyword to limit the number of moments:
;               Maxmoment = 1  Calculate only the mean.
;               Maxmoment = 2  Calculate the mean and variance.
;               Maxmoment = 3  Calculate the mean, variance, and skewness.
;               Maxmoment = 4  Calculate the mean, variance, skewness,
;                              and kurtosis (the default).
;
;       NAN:    Treat NaN elements as missing data.
;               (Due to problems with IEEE support on some platforms,
;                infinite data will be treated as missing as well. )
;
; EXAMPLE:
;       Define the N-element vector of sample data.
;         x = [65, 63, 67, 64, 68, 62, 70, 66, 68, 67, 69, 71, 66, 65, 70]
;       Compute the mean, variance, skewness and kurtosis.
;         result = moment(x)
;       The result should be the 4-element vector: 
;       [66.7333, 7.06667, -0.0942851, -1.18258]
;  
;
; PROCEDURE:
;       MOMENT computes the first four "moments" about the mean of an N-element
;       vector of sample data. The computational formulas are given in the IDL 
;       Reference Guide. 
;
; REFERENCE:
;       APPLIED STATISTICS (third edition)
;       J. Neter, W. Wasserman, G.A. Whitmore
;       ISBN 0-205-10328-6
;
; MODIFICATION HISTORY:
;       Written by:  GGS, RSI, August 1994
;       Modified:    GGS, RSI, September 1995
;                    Added DOUBLE keyword. 
;                    Added checking for N identical elements. 
;                    Added support for IEEE NaN (Not a Number).
;                    Modified variance formula.
;       Modified:    GGS, RSI, April 1996
;                    Modified keyword checking and use of double precision. 
;                    GSL, RSI, August 1997
;                    Added Maxmoment keyword.
;       Modified:    Wed Jan 28 13:28:07 1998, Scott Lett, RSI Added
;                    NAN keyword.
;-
FUNCTION Moment, X, Double = Double, Mdev = Mdev, Sdev = Sdev, $
                 Maxmoment = Maxmoment, NaN = nan

ON_ERROR, 2

if keyword_set( nan ) then begin ;If NaN set, remove NaNs and recurse.
    whereNotNaN = where( finite(X), nanCount)
    if nanCount gt 0 then begin
        return, moment( X[whereNotNan], Double = Double, Mdev = Mdev, $
                        Sdev = Sdev, Maxmoment = Maxmoment )
    endif
endif

TypeX = SIZE(X)

IF NOT keyword_set( Maxmoment ) THEN Maxmoment = 4

IF Maxmoment gt 1 and N_ELEMENTS(x) lt 2 THEN $ ;Check length.
  MESSAGE, "X array must contain 2 OR more elements."

  ;If the DOUBLE keyword is not set then the internal precision and
  ;result are identical to the type of input.
IF N_ELEMENTS(Double) EQ 0 THEN $
  Double = (TypeX[TypeX[0]+1] EQ 5 OR TypeX[TypeX[0]+1] EQ 9)

nX = TypeX[TypeX[0]+2]
Mean = TOTAL(X, Double = Double) / nX

Var  = !VALUES.F_NAN
Skew = !VALUES.F_NAN
Kurt = !VALUES.F_NAN

IF Maxmoment GT 1 THEN BEGIN    ; Calculate higher moments.
    Resid = X - Mean
;   Var = TOTAL(Resid^2, Double = Double) / (nX-1.0);Simple formula

; Numerically-stable "two-pass" formula, which offers less
; round-off error. Page 613, Numerical Recipes in C.
    Var = (TOTAL(Resid^2, Double = Double) - $
           (TOTAL(Resid, Double = Double)^2)/nX)/(nX-1.0)

;Mean absolute deviation (returned through the Mdev keyword).
    if arg_present(Mdev) then $
      Mdev = TOTAL(ABS(Resid), Double = Double) / nX
    
; Standard deviation (returned through the Sdev keyword).
    Sdev = SQRT(Var)

    IF Sdev NE 0 THEN BEGIN	;Skew & kurtosis defined
        IF Maxmoment GT 2 THEN $
          Skew = TOTAL(Resid^3, Double = Double) / (nX * Sdev ^ 3)
; The "-3" term makes the kurtosis value zero for normal distributions.
; Positive values of the kurtosis (lepto-kurtic) indicate pointed or
; peaked distributions; Negative values (platy-kurtic) indicate flat-
; tened or non-peaked distributions.
        IF Maxmoment GT 3 THEN $
          Kurt = TOTAL(Resid^4, Double = Double) / (nX * Sdev ^ 4) - 3.0
    ENDIF
ENDIF
RETURN, [Mean, Var, Skew, Kurt]
END

Return to the shapelets web page or the code help menu.