Table of Contents

Name

ini_moment, accum_moment, decr_moment, reset_moment, show_moment, n_moment, sum_moment, mean_moment, sigma_moment, skewness_moment, kurtosis_moment, mad_moment, mard_moment, robust_moment, min_moment, max_moment - various (moving) moment and minmax routines

Synopsis


#include <moment.h>
void ini_moment(m, mom, ndat)void accum_moment(m, x, w)void decr_moment(m,
x, w)void reset_moment(m)
real show_moment(m, mom)int n_moment(m)real sum_moment(m)real
mean_moment(m)real median_moment(m)real sigma_moment(m)real skewness_moment(m)real
kurtosis_moment(m)real h3_moment(m)real h4_moment(m)real mad_moment(m)real
mard_moment(m)
real min_moment(m);real max_moment(m);
void compute_robust_moment(m);int
n_robust_moment(m);real mean_robust_moment(m);real median_robust_moment(m);real
sigma_robust_moment(m);
Moment *m;int mom, ndat;real x, w;

Description

moment is a set of functions to compute the moments of a set of real values, but also keep track of the datamin and datamax. The routines are written in C but in an object-oriented fashion as to try and keep the use from the internal data-structures. The Moment can (and should) be treated as an opaque datastructure. All information should be obtained through the interface routines described below.

If only a datamin/max is needed, setting mom<0 can be used to prevent the more expensive moment calculations.

Moving or running averages (or moments) can be done by supplying ndat>0 to ini_moment. It will keep a memory of the last ndat data values and the moments now become running moments.

Note that the median_moment can only be used in x (the weights are ignored) and moving moment where ndat>0.

mean_moment returns the mean value, where sigma_moment returns the square root of the variance

robust_moment’s can also be computed. These are based on an algorithm to remove the outliers. It first computes the distance between the 1st and 3rd quartile, then removes all points further than 1.5 times that distance from the 1st and 3rd quartile. After calling compute_robust_moment, the mean, median and sigma can be retrieved as those computations are done in a temporary Moment structure. Caveat: for large datasets and/or moving moments, robust_moment’s can become expensive because of the need for sorting to get the quartiles.

mad_moment computes the Median Absolute Deviation (MAD), arguably a better measure for the Standard Deviation. As with the robust moments, it needs to keep a copy of the data available. MAD is formally RMS/1.4826. Related is mard_moment, the Mean Absolute Relative Difference (MARD).

Moments

A note on the h3 and h4 moments, somewhat peculiar to astronomy. See S2.4 in van der Marel & Franx (1993)

Examples

The following code computes a weighted mean and dispersion of a set of points:
    real x[100], w[100];
    int  i,n=100;
    Moment m;
    ...
    ini_moment(&m,2,0);      /* up to 2nd order moment, and using no circular
buffer */
    for (i=0; i<n; i)
    accum_moment(&m,x[i],w[i]);
    printf("Mean: %g   Dispersion: %g\n",
    mean_moment(&m), sigma_moment(&m));
    compute_robust_moment(&m);
    printf("Robust Mean: %g   Dispersion: %g\n",
    robust_mean_moment(&m), robust_sigma_moment(&m));

Moment Structure

A simple data structure (referred to as moment in the above SYNOPSIS) is used to communicate between different routines:
typedef struct { 
    int mom;
    int n;
    real *sum;
    real datamin, datamax;
    int ndat;
    int idat;
    real *dat;
    real *sum;
} Moment;
from the standard include file moment.h.

Bugs

When decr_moment is used, the data min/max is not correct. Only with ndat>0 for moving moments can it be recomputed correctly.

See Also

grid(3NEMO)
http://apophenia.info

Author

Peter Teuben

Files


~/src/kernel/misc    moment.c

Update History


30-oct-93    Created       PJT
8-nov-93    fixed init bug     PJT
13-jun-95    added decr_moment    PJT
2-feb-05    added moving moments    PJT
2-mar-11    added h3,h4    PJT
24-apr-13    documented robust statistics    PJT
16-jan-14    added MAD    PJT
11-jun-14    clarified MAD and MARD (the old MAD was really MARD)    PJT


Table of Contents