f14-1

604 Chapter14. StatisticalDescriptionofDataSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Intheothercategory,model-dependentstatistics, welumpthewholesubjectof fitting data to a theory, parameter estimation, least-squares fits, and so on. Thosesubjects are introduced in Chapter 15. Section14.1dealswithso-called measuresofcentraltendency ,themomentsof adistribution,themedianandmode. In §14.2we learntotest whetherdifferentdata sets are drawn from distributions with different values of these measures of central tendency. Thisleadsnaturally,in §14.3,tothemoregeneralquestionofwhethertwo distributions can be shown to be (significantly) different. In§14.4– §14.7, we deal with measures of association for two distributions. We want to determine whether two variables are “correlated” or “dependent” onone another. If they are, we want to characterize the degree of correlation in some simple ways. The distinction between parametric and nonparametric (rank) methods is emphasized. Section 14.8 introduces the concept of data smoothing, and discusses the particular case of Savitzky-Golay smoothing filters. This chapter draws mathematically on the material on special functions that was presented in Chapter 6, especially §6.1–§6.4. You may wish, at this point, to review those sections. CITED REFERENCES AND FURTHER READING: Bevington, P.R. 1969, Data Reduction and Error Analysis for the Physical Sciences (New York: McGraw-Hill). Stuart, A., andOrd, J.K. 1987, Kendall’sAdvanced Theory of Statistics , 5th ed. (London:Griffin and Co.) [previous eds. published as Kendall, M., and Stuart, A., The Advanced Theory of Statistics ]. Norusis,M.J.1982, SPSSIntroductoryGuide:BasicStatisticsandOperations ;and1985, SPSS- X Advanced Statistics Guide (New York: McGraw-Hill). Dunn, O.J., andClark, V.A. 1974, AppliedStatistics: Analysis of Variance andRegression (New York: Wiley). 14.1 Moments of a Distribution: Mean, Variance, Skewness, and So Forth Whenasetofvalueshasasufficientlystrongcentraltendency,thatis,atendency to cluster around some particular value, then it may be useful to characterize the set by a few numbers that are related to its moments, the sums of integer powers of the values. Best known is the meanof the values x1,...,x N, x=1 NN/summationdisplay j=1xj (14.1.1 ) which estimates the value around which central clustering occurs. Note the use of anoverbartodenotethemean;anglebracketsareanequallycommonnotation,e.g., /angbracketleftx/angbracketright. You should be aware that the mean is not the only available estimator of this 14.1MomentsofaDistribution: Mean,Variance,Skewness 605Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).quantity, nor is it necessarily the best one. For values drawn from a probability distribution with very broad “tails,” the mean may convergepoorly, or not at all, asthe number of sampled points is increased. Alternative estimators, the medianand themode, are mentioned at the end of this section. Having characterized a distribution’s central value, one conventionally next characterizes its “width” or “variability” around that value. Here again, more than one measure is available. Most common is the variance, Var(x 1...x N)=1 N−1N/summationdisplay j=1(xj−x)2(14.1.2 ) or its square root, the standard deviation , σ(x1...x N)=/radicalbig Var(x1...x N)( 14.1.3 ) Equation (14.1.2) estimates the mean squared deviation of xfrom its mean value. There is a long story about why the denominator of (14.1.2) is N−1instead of N. If you have never heard that story, you may consult any good statistics text. Here we will be content to note that the N−1shouldbe changed to Nif you are ever in the situation of measuring the variance of a distribution whose mean xis known a priorirather than being estimated from the data. (We might also comment that if the differencebetween NandN−1ever matters to you, then you are probably up to no good anyway — e.g., trying to substantiate a questionable hypothesis with marginal data.) As the mean depends on the first moment of the data, so do the variance and standard deviation depend on the second moment. It is not uncommon, in real life, to be dealing with a distribution whose second moment does not exist (i.e., is infinite). In this case, the variance or standard deviation is useless as a measure of the data’s width around its central value: The values obtained from equations(14.1.2) or (14.1.3) will not converge with increased numbers of points, nor show anyconsistencyfromdatasettodataset drawnfromthesamedistribution. Thiscan occurevenwhenthewidthofthepeaklooks,byeye,perfectlyfinite. Amorerobustestimatorofthewidthisthe averagede viationormeanabsolutedeviation ,definedby ADev (x 1...x N)=1 NN/summationdisplay j=1|xj−x| (14.1.4 ) One often substitutes the sample median xmedforxin equation (14.1.4). For any fixed sample, the median in fact minimizes the mean absolute deviation. Statisticians have historically sniffed at the use of (14.1.4) instead of (14.1.2), since the absolute value brackets in (14.1.4) are “nonanalytic” and make theorem- provingdifficult. In recent years, however,the fashion has changed,and the subject ofrobust estimation (meaning, estimation for broad distributions with significant numbers of “outlier” points) has become a popular and important one. Higher moments, or statistics involving higher powers of the input data, are almost always less robust than lower moments or statistics that involve only linear sums or (the lowest moment of all) counting. 606 Chapter14. StatisticalDescriptionofDataSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).(b) (a)Skewness negative positivepositive (leptokurtic)negative (platykurtic)Kurtosis Figure 14.1.1. Distributions whose third and fourth moments are signi ficantly different from a normal (Gaussian) distribution. (a) Skewness or third moment. (b) Kurtosis or fourth moment. That being the case, the skewness orthird moment , and the kurtosisorfourth momentshould be used with caution or, better yet, not at all. Theskewnesscharacterizesthedegreeofasymmetryofadistributionaroundits mean. While the mean, standard deviation, and average deviation are dimensional quantities, that is, have the same units as the measured quantities xj, the skewness is conventionally de fined in such a way as to make it nondimensional . It is a pure numberthat characterizesonlythe shapeofthedistribution. Theusualde finitionis Skew (x1...x N)=1 NN/summationdisplay j=1/bracketleftbiggxj−x σ/bracketrightbigg3 (14.1.5 ) where σ=σ(x1...x N)is the distribution ’sstandarddeviation(14.1.3). A positive value of skewness signi fies a distribution with an asymmetric tail extending out towards more positive x; a negativevalue signi fies a distribution whose tail extends out towards more negative x(see Figure 14.1.1). Of course, any set of Nmeasured values is likely to give a nonzero value for (14.1.5),eveniftheunderlyingdistributionisinfactsymmetrical(haszeroskewness). For (14.1.5) to be meaningful, we need to have some idea of itsstandard deviation as an estimator of the skewness of the underlying distribution. Unfortunately, that dependsonthe shapeof theunderlyingdistribution,and rathercritically onits tails! For the idealized case of a normal(Gaussian) distribution,the standard deviationof(14.1.5) is approximately/radicalbig 15/Nwhen xis the true mean, and/radicalbig 6/Nwhen it is estimated by the sample mean, (14.1.1). In real life it is good practice to believe in skewnesses only when they are several or many times as large as this. The kurtosis is also a nondimensional quantity. It measures the relative peakedness or flatness of a distribution. Relative to what? A normal distribution, what else! A distribution with positive kurtosis is termed leptokurtic ; the outline of the Matterhorn is an example. A distribution with negative kurtosis is termed platykurtic ; the outline of a loaf of bread is an example. (See Figure 14.1.1.) And, as you no doubt expect, an in-between distribution is termed mesokurtic . The conventional de finition of the kurtosis is Kurt (x1...x N)=  1 NN/summationdisplay j=1/bracketleftbiggxj−x σ/bracketrightbigg4  −3( 14.1.6 ) 14.1MomentsofaDistribution: Mean,Variance,Skewness 607Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).where the −3term makes the value zero for a normal distribution. Thestandarddeviationof(14.1.6)asanestimatorofthekurtosisofanunderlying normal distribution is/radicalbig 96/Nwhen σis the true standard deviation, and/radicalbig 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 standarddeviation of (14.1.6) as an estimator is effectively in finite. Calculationofthequantitiesde finedinthissectionis perfectlystraightforward. Many textbooks use the binomial theorem to expand out the de finitions into sums of various powers of the data, e.g., the familiar Var(x 1...x N)=1 N−1  N/summationdisplay j=1x2 j −Nx2 ≈x2−x2(14.1.7 ) butthiscanmagnifytheroundofferrorbyalargefactorandisgenerallyunjusti fiable 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...x N)by Var(x1...x N)=1 N−1  N/summationdisplay j=1(xj−x)2−1 N N/summationdisplay j=1(xj−x) 2  (14.1.8 ) The second sum would be zero if xwere exact, but otherwise it does a good job of correcting the roundoff error in the first term. SUBROUTINE moment(data,n,ave,adev,sdev,var,skew,curt) INTEGER nREAL adev,ave,curt,sdev,skew,var,data(n) G i v e na na r r a yo f data(1:n) , this routine returns its mean ave , average deviation adev , standard deviation sdev ,v a r i a n c e var ,s k e w n e s s skew ,a n dk u r t o s i s curt . INTEGER jREAL p,s,ep if(n.le.1)pause ’n must be at least 2 in moment’ s=0. First pass to get the mean. do 11j=1,n s=s+data(j) enddo 11 ave=s/n adev=0. Second pass to get the first (absolute), second, third, and fourth moments of the deviation from the mean. var=0. skew=0.curt=0.ep=0. do 12j=1,n s=data(j)-aveep=ep+s adev=adev+abs(s) p=s*svar=var+pp=p*s skew=skew+p p=p*scurt=curt+p enddo 12 608 Chapter14. StatisticalDescriptionofDataSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).adev=adev/n Put the pieces together according to the conventional definitions. var=(var-ep**2/n)/(n-1) Corrected two-pass formula. sdev=sqrt(var) if(var.ne.0.)then skew=skew/(n*sdev**3) curt=curt/(n*var**2)-3. else pause ’no skew or kurtosis when zero variance in moment’ endif return END Semi-Invariants The mean and variance of independent random variables are additive: If xand yare drawn independently from two, possibly different, probability distributions, then (x+y)= x+yVar(x+y)=Var(x)+Var(x)( 14.1.9 ) Higher moments are not, in general, additive. However, certain combinations of them, calledsemi-invariants , are in fact additive. If the centered moments of a distribution are denoted M k, M k≡/angbracketleftBig (xi−x)k/angbracketrightBig (14.1.10 ) so that, e.g., M2=Var(x), then the first few semi-invariants, denoted Ikare given by I2=M2 I3=M3 I4=M4−3M2 2 I5=M5−10M2M3 I6=M6−15M2M4−10M2 3+3 0 M3 2(14.1.11 ) Noticethattheskewness andkurtosis,equations(14.1.5)and(14.1.6)aresimplepowers of the semi-invariants, Skew (x)= I3/I3/2 2Kurt (x)= I4/I2 2 (14.1.12 ) A Gaussian distribution has all its semi-invariants higher than I2equal to zero. A Poisson distribution has all of its semi-invariants equal to its mean. For more details, see [2]. Median and Mode The median of a probability distribution function p(x)is the value xmedfor which larger and smaller values of xare equally probable: /integraldisplayxmed −∞p(x)dx=1 2=/integraldisplay∞ xmedp(x)dx (14.1.13 ) The median of a distribution is estimated from a sample of values x1,..., xNbyfinding that value xiwhich has equal numbers of values above it and below it. Of course, this is not possible when Nis even. In that case it is conventional to estimate the median as the mean of the unique twocentral values. If the values xjj=1,...,Nare sorted into ascending (or, for that matter, descending) order, then the formula for the median is xmed =/braceleftbiggx(N+1) /2,N odd 1 2(xN/2+x(N/2)+1 ),Neven(14.1.14 ) 14.2DoTwoDistributionsHavetheSameMeansorVariances? 609Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).If a distribution has a strong central tendency, so that most of its area is under a single peak, then the median is an estimator of the central value. It is a morerobust estimator than the mean is: The median fails as an estimator only if the area in the tails is large, while the mean fails if the first moment of the tails is large; it is easy to construct examples where the first moment of the tails is large even though their area is negligible. Tofind the median of a set of values, one can proceed by sorting the set and thenapplying(14.1.14). Thisisaprocessoforder NlogN. Youmightrightlythink that this is wasteful, since it yields much more information than just the median (e.g., the upper and lower quartile points, the deciles, etc.). In fact, we saw in§8.5 that the element x (N+1) /2can be located in of order Noperations. Consult that section for routines. Themodeof a probability distribution function p(x)is the value of xwhere it takesonamaximumvalue. Themodeisusefulprimarilywhenthereisasingle,sharp maximum, in which case it estimates the central value. Occasionally, a distribution will bebimodal, with two relative maxima; then one may wish to know the two modes individually. Note that, in such cases, the mean and median are not very useful, since they will give only a “compromise ”value between the two peaks. CITED REFERENCES AND FURTHER READING: Bevington, P.R. 1969, Data Reduction and Error Analysis for the Physical Sciences (New York: McGraw-Hill), Chapter 2. Stuart, A., andOrd, J.K. 1987, Kendall’sAdvanced Theory of Statistics , 5th ed. (London:Griffin and Co.) [previous eds. published as Kendall, M., and Stuart, A., The Advanced Theory of Statistics ], vol. 1, §10.15 Norusis,M.J.1982, SPSSIntroductoryGuide:BasicStatisticsandOperations ;and1985, SPSS- X Advanced Statistics Guide (New York: McGraw-Hill). Chan,T.F.,Golub,G.H.,andLeVeque,R.J.1983, AmericanStatistician ,vol.37,pp.242–247.[1] Cram´er, H. 1946, Mathematical Methods of Statistics (Princeton: Princeton University Press), §15.10. [2] 14.2 Do Two Distributions Have the Same Means or Variances? Not uncommonly we want to know whether two distributions have the same mean. For example, a first set of measured values may have been gathered before some event,a secondset after it. We want to knowwhetherthe event,a “treatment ” or a“change in a control parameter, ”made a difference. Ourfirst thoughtis to ask “howmanystandarddeviations ”onesamplemeanis from the other. That numbermay in fact be a useful thing to know. It does relate to the strength or “importance ”of a difference of means if that difference is genuine . However, by itself, it says nothing about whether the difference isgenuine, that is, statistically signi ficant. A difference of means can be very small compared to the standard deviation, and yet very signi ficant, if the number of data points is large. Conversely, a difference may be moderately large but not signi ficant, if the data