#include <moment.h> void ini_moment(m, mom, ndat)void slow_moment(m, slow)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 sratio_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 min_robust_moment(m);real max_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;
If only a datamin/max is needed, setting mom<0 can be used to prevent the more expensive moment calculations. See also slow_moment.
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. This is also the only way the slow_moment (more accurate for large values) will work.
Note that the median_moment can only be used in x (the weights are ignored) and moving moment where ndat>0.
slow_moment can be set with slow=TRUE to force re-computing higher order moments based on the value of the mean. This only will work when ndat>0.
mean_moment returns the mean value, where sigma_moment returns the square root of the variance
sratio_moment returns a measure of how much signal there is compared to the noise, by compting (SumP+SumN)/(SumP-SumN). For pure noise this should be near 0, and for pure signal, where SumN=0, near 1.
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).
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",
mean_robust_moment(&m), sigma_robust_moment(&m));
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.
http://apophenia.info
~/src/kernel/misc moment.c
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 12-jul-20 added min/max for robust moment PJT 14-nov-21 added sratio PJT 18-jan-25 added slow option PJT