14.1 Moments of a Distribution:Mean.Variance,Skewness 613 The standard deviation of(14.1.6)as an estimator ofthe kurtosis ofan underlying normal distribution is v96/N when o is the true standard deviation,and v24/N when it is the sample estimate (14.1.3).However,the kurtosis depends on such a high moment that there are many real-life distributions for which the standard deviation of(14.1.6)as an estimator is effectively infinite. Calculation of the quantities defined in this section is perfectly straightforward. Many textbooks use the binomial theorem to expand out the definitions into sums of various powers of the data,e.g.,the familiar Var(1...IN) [位 (14.1.7) but this can magnify the roundofferror by a large factor and is generally unjustifiable in terms of computing speed.A clever way to minimize roundoff error,especially for large samples,is to use the corrected two-pass algorithm [11:First calculate then calculate Var(z1...xN)by 、、1 (Nort server 令 Var(x1..xN (14.1.8) 景S%s3dN America The second sum would be zero if were exact,but otherwise it does a good job of correcting the roundoff error in the first term. Programs #include <math.h> SCIENTIFIC void moment(float data[],int n,float *ave,float *adev,float *sdev, float *var,float *skew,float *curt) Given an array of data[1..n],this routine returns its mean ave,average deviation adev, standard deviation sdev,variance var,skewness skew,and kurtosis curt. void nrerror(char error._text[☐)； 1920 COMPUTING(ISBN int j; float ep=0.0,s,P; if (n <1)nrerror("n must be at least 2 in moment"); 8=0.0； First pass to get the mean. for (j=1;j<=n;j++)s +data[j]; Numerical Recipes 021 43108 来aVe世g/n: *adev=(*var)=(*skew)=(*curt)=0.0; Second pass to get the first (absolute),sec- for (j=1;j<=n;j++){ ond,third,and fourth moments of the (outside *adev +fabs(s=data[j]-(tave)); deviation from the mean. 膜 ep +s; North oftware. 米Var+=(p=s*s)； *skew +(p *s); *curt +(p *s) *adev /n; *var=(*var-ep*ep/n)/(n-1); Corrected two-pass formula. *sdev=sgrt(*var); Put the pieces together according to the con- if (*var){ ventional definitions. *skew /(n*(*var)*(*sdev)); *curt=(*curt)/(n*(*var)*(*var))-3.0; else nrerror("No skew/kurtosis when variance 0 (in moment)");14.1 Moments of a Distribution: Mean, Variance, Skewness 613 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). The standard deviation of (14.1.6) as an estimator of the kurtosis of an underlying normal distribution is 96/N when σ is the true standard deviation, and 24/N when it is the sample estimate (14.1.3). However, the kurtosis depends on such a high moment that there are many real-life distributions for which the standard deviation of (14.1.6) as an estimator is effectively infinite. Calculation of the quantities defined in this section is perfectly straightforward. Many textbooks use the binomial theorem to expand out the definitions into sums of various powers of the data, e.g., the familiar Var(x1 ...xN ) = 1 N − 1    
N j=1 x2 j   − Nx2   ≈ x2 − x2 (14.1.7) but this can magnify the roundoff error by a large factor and is generally unjustifiable in terms of computing speed. A clever way to minimize roundoff error, especially for large samples, is to use the corrected two-pass algorithm [1]: First calculate x, then calculate Var(x1 ...xN ) by Var(x1 ...xN ) = 1 N − 1   
N j=1 (xj − x) 2 − 1 N  
N j=1 (xj − x)   2    (14.1.8) The second sum would be zero if x were exact, but otherwise it does a good job of correcting the roundoff error in the first term. #include <math.h> void moment(float data[], int n, float *ave, float *adev, float *sdev, float *var, float *skew, float *curt) Given an array of data[1..n], this routine returns its mean ave, average deviation adev, standard deviation sdev, variance var, skewness skew, and kurtosis curt. { void nrerror(char error_text[]); int j; float ep=0.0,s,p; if (n <= 1) nrerror("n must be at least 2 in moment"); s=0.0; First pass to get the mean. for (j=1;j<=n;j++) s += data[j]; *ave=s/n; *adev=(*var)=(*skew)=(*curt)=0.0; Second pass to get the first (absolute), sec￾ond, third, and fourth moments of the deviation from the mean. for (j=1;j<=n;j++) { *adev += fabs(s=data[j]-(*ave)); ep += s; *var += (p=s*s); *skew += (p *= s); *curt += (p *= s); } *adev /= n; *var=(*var-ep*ep/n)/(n-1); Corrected two-pass formula. *sdev=sqrt(*var); Put the pieces together according to the con￾if (*var) { ventional definitions. *skew /= (n*(*var)*(*sdev)); *curt=(*curt)/(n*(*var)*(*var))-3.0; } else nrerror("No skew/kurtosis when variance = 0 (in moment)"); }