f14-7

640 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).14.7 Do Two-DimensionalDistributions Differ? We here discuss a useful generalization of the K–S test ( §14.3) totwo-dimensional distributions. ThisgeneralizationisduetoFasanoandFranceschini [1],avariant onanearlier idea due to Peacock [2]. In a two-dimensional distribution, each data point is characterized by an (x, y)pair of values. An example near to our hearts is that each of the 19 neutrinos that were detectedfrom Supernova 1987A is characterized by a time t iand by an energy Ei(see[3]). We might wish to know whether these measured pairs (ti,E i),i=1 ...19are consistent with a theoretical model that predicts neutrino flux as a function of both time and energy — that is,a two-dimensional probability distribution in the (x, y)[here, (t, E)] plane. That would be a one-sample test. Or, given two sets of neutrino detections, from two comparable detectors, we might want to know whether they are compatible with each other, a two-sample test. In the spirit of the tried-and-true, one-dimensional K–S test, we want to range over the(x, y)plane in search of some kind of maximum cumulative difference between two two-dimensional distributions. Unfortunately, cumulative probability distribution is notwell-defined in more than one dimension! Peacock’s insight was that a good surrogate istheintegrated probability in each of four natural quadrants around a given point (x i,y i), namely the total probabilities (or fraction of data) in (x>x i,y > y i),(x<x i,y > y i), (x<x i,y < y i),(x>x i,y < y i). The two-dimensional K–S statistic Dis now taken to be the maximum difference (ranging both over data points and over quadrants) of thecorresponding integrated probabilities. When comparing two data sets, the value of Dmay depend on which data set is ranged over. In that case, define an effective Das the average of the two values obtained. If you are confused at this point about the exact definition of D, don’t fret; the accompanying computer routines amount to a precise algorithmic definition. Figure14.7.1gives afeelingforwhatisgoing on. The65trianglesand 35squares seem to have somewhat different distributions in the plane. The dotted lines are centered on thetriangle that maximizes the Dstatistic; the maximum occurs in the upper-left quadrant. That quadrant contains only 0.12 of all the triangles, but it contains 0.56 of all the squares. Thevalue of Dis thus 0.44. Is this statistically significant? Evenforfixedsamplesizes,itisunfortunately notrigorously truethatthedistributionof Dinthenullhypothesisisindependentoftheshapeofthetwo-dimensionaldistribution. Inthis respectthetwo-dimensionalK–Stestisnotasnaturalasitsone-dimensionalparent. However,extensive Monte Carlo integrations have shown that the distribution of the two-dimensionalDisvery nearly identical for even quite different distributions, as long as they have the same coefficient of correlation r, defined in the usual way by equation (14.5.1). In their paper, FasanoandFranceschinitabulateMonteCarloresultsfor(whatamountsto)thedistributionofDas a function of (of course) D, sample size N, and coefficient of correlation r. Analyzing their results, one finds that the significance levels for the two-dimensional K–S test can besummarized by the simple, though approximate, formulas, Probability (D>observed )= Q KS/parenleftbigg √ ND 1+√ 1−r2(0.25−0.75/√ N)/parenrightbigg (14.7.1 ) for the one-sample case, and the same for the two-sample case, but with N=N1N2 N1+N2. (14.7.2 ) The above formulas are accurate enough when N>∼20, and when the indicated probability (significance level) is less than (more significant than) 0.20or so. When the indicated probability is >0.20, its value may not be accurate, but the implication that the data and model (or two data sets) are not significantly different is certainly correct. Noticethat in the limit of r→1(perfect correlation), equations (14.7.1) and (14.7.2) reduce to equations (14.3.9) and (14.3.10): The two-dimensional data lie on a perfect straight line, andthe two-dimensional K–S test becomes a one-dimensional K–S test. 14.7DoTwo-DimensionalDistributionsDiffer? 641Sample 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).−3 −2 −1 012 3−3−2−10123 .11|.09.65|.26 .12|.09.12|.56 Figure 14.7.1. Two-dimensional distributions of 65 triangles and 35 squares. The two-dimensional K –S testfinds that point one of whose quadrants (shown by dotted lines) maximizes the difference between fraction of triangles and fraction of squares. Then, equation (14.7.1) indicates whether the difference isstatisticallysigni ficant,i.e.,whetherthetrianglesandsquaresmusthavedifferentunderlyingdistributions. The signi ficance level for the data in Figure 14.7.1, by the way, is about 0.001. This establishes to a near-certainty that the triangles and squares were drawn from differentdistributions. (As in fact they were.) Of course, if you do not want to rely on the Monte Carlo experiments embodied in equation (14.7.1), you can do your own: Generate a lot of synthetic data sets from yourmodel, each one with the same number of points as the real data set. Compute Dfor each synthetic data set, using the accompanying computer routines (but ignoring their calculatedprobabilities), and count what fraction of the time these synthetic D’s exceed the Dfrom the real data. That fraction is your signi ficance. Onedisadvantageofthetwo-dimensionaltests,bycomparisonwiththeirone-dimensional progenitors, is that the two-dimensional tests require of order N 2operations: Two nested loops of order Ntake the place of an Nlog Nsort. For small computers, this restricts the usefulness of the tests to Nless than several thousand. We now give computer implementations. The one-sample case is embodied in the routine ks2d1s(that is, 2-dimensions, 1-sample). This routine calls a straightforward utility routine quadctto count points in the four quadrants, and it calls a user-supplied routine quadvlthat must be capable of returning the integrated probability of an analytic model in each of four quadrants around an arbitrary (x, y)point. A trivial sample quadvlis shown; realistic quadvls can be quite complicated, often incorporating numerical quadratures over 642 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).analytic two-dimensional distributions. SUBROUTINE ks2d1s(x1,y1,n1,quadvl,d1,prob) INTEGER n1REAL d1,prob,x1(n1),y1(n1) EXTERNAL quadvl C USES pearsn,probks,quadct,quadvl Two-dimensional Kolmogorov-Smirnov test of one sample against a model. Given the x and ycoordinates of n1data points in arrays x1(1:n1) andy1(1:n1) ,a n dg i v e na user-supplied function quadvlthat exemplifies the model, this routine returns the two- dimensional K-S statistic as d1, and its significance level as prob. Small values of prob show that the sample is significantly different from the model. Note that the test is slightly distribution-dependent, so probis only an estimate. INTEGER jREAL dum,dumm,fa,fb,fc,fd,ga,gb,gc,gd,r1,rr,sqen,probks d1=0.0 do 11j=1,n1 Loop over the data points. call quadct(x1(j),y1(j),x1,y1,n1,fa,fb,fc,fd)call quadvl(x1(j),y1(j),ga,gb,gc,gd) d1=max(d1,abs(fa-ga),abs(fb-gb),abs(fc-gc),abs(fd-gd)) For both the sample and the model, the distribution is integrated in each of four quad-rants, and the maximum difference is saved. enddo 11 call pearsn(x1,y1,n1,r1,dum,dumm) Get the linear correlation coefficient r1. sqen=sqrt(float(n1))rr=sqrt(1.0-r1**2) Estimate the probability using the K-S probability function probks. prob=probks(d1*sqen/(1.0+rr*(0.25-0.75/sqen)))returnEND SUBROUTINE quadct(x,y,xx,yy,nn,fa,fb,fc,fd) INTEGER nn REAL fa,fb,fc,fd,x,y,xx(nn),yy(nn) G i v e na no r i g i n ( x,y), and an array of nnpoints with coordinates xxandyy, count how many of them are in each quadrant around the origin, and return the normalized frac- tions. Quadrants are labeled alphabetically, counterclockwise from the upper right. Used byks2d1sandks2d2s. INTEGER k,na,nb,nc,ndREAL ff na=0 nb=0nc=0nd=0 do 11k=1,nn if(yy(k).gt.y)then if(xx(k).gt.x)then na=na+1 else nb=nb+1 endif else if(xx(k).gt.x)then nd=nd+1 else nc=nc+1 endif endif enddo 11 ff=1.0/nn fa=ff*na fb=ff*nb 14.7DoTwo-DimensionalDistributionsDiffer? 643Sample 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).fc=ff*nc fd=ff*nd return END SUBROUTINE quadvl(x,y,fa,fb,fc,fd) REAL fa,fb,fc,fd,x,y This is a sample of a user-supplied routine to be used with ks2d1s. In this case, the model distribution is uniform inside the square −1<x< 1,−1<y< 1. In general this routine should return, for any point (x,y), the fraction of the total distribution in each of the four quadrants around that point. The fractions, fa,fb,fc,a n d fd, must add up to 1. Quadrants are alphabetical, counterclockwise from the upper right. REAL qa,qb,qc,qd qa=min(2.,max(0.,1.-x)) qb=min(2.,max(0.,1.-y))qc=min(2.,max(0.,x+1.))qd=min(2.,max(0.,y+1.)) fa=0.25*qa*qb fb=0.25*qb*qcfc=0.25*qc*qdfd=0.25*qd*qa return END Theroutine ks2d2sisthetwo-samplecaseofthetwo-dimensionalK –Stest. Italsocalls quadct,pearsn, and probks. Being a two-sample test, itdoes not need an analytic model. SUBROUTINE ks2d2s(x1,y1,n1,x2,y2,n2,d,prob) INTEGER n1,n2 REAL d,prob,x1(n1),x2(n2),y1(n1),y2(n2) C USES pearsn,probks,quadct Two-dimensional Kolmogorov-Smirnov test on two samples. Given the xand ycoordinates of the first sample as n1values in arrays x1(1:n1) andy1(1:n1) , and likewise for the second sample, n2valuesin arrays x2andy2, this routine returns the two-dimensional, two- sample K-S statistic as d, and its significance level as prob. Small values of probshow that the two samples are significantly different. Note that the test is slightly distribution- dependent, so probis only an estimate. INTEGER jREAL d1,d2,dum,dumm,fa,fb,fc,fd,ga,gb,gc,gd,r1,r2,rr, * sqen,probks d1=0.0do 11j=1,n1 First, use points in the first sample as origins. call quadct(x1(j),y1(j),x1,y1,n1,fa,fb,fc,fd) call quadct(x1(j),y1(j),x2,y2,n2,ga,gb,gc,gd) d1=max(d1,abs(fa-ga),abs(fb-gb),abs(fc-gc),abs(fd-gd)) enddo 11 d2=0.0do 12j=1,n2 Then, use points in the second sample as origins. call quadct(x2(j),y2(j),x1,y1,n1,fa,fb,fc,fd)call quadct(x2(j),y2(j),x2,y2,n2,ga,gb,gc,gd) d2=max(d2,abs(fa-ga),abs(fb-gb),abs(fc-gc),abs(fd-gd)) enddo 12 d=0.5*(d1+d2) Average the K-S statistics. sqen=sqrt(float(n1)*float(n2)/float(n1+n2)) call pearsn(x1,y1,n1,r1,dum,dumm) Getthelinearcorrelationcoefficientforeachsample. call pearsn(x2,y2,n2,r2,dum,dumm)rr=sqrt(1.0-0.5*(r1**2+r2**2)) Estimate the probability using the K-S probability function probks. prob=probks(d*sqen/(1.0+rr*(0.25-0.75/sqen)))return END 644 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).CITED REFERENCES AND FURTHER READING: Fasano, G. and Franceschini, A. 1987, Monthly Notices of the Royal Astronomical Society , vol. 225, pp. 155–170. [1] Peacock,J.A.1983, MonthlyNoticesoftheRoyalAstronomicalSociety ,vol.202,pp.615–627.[2] Spergel, D.N., Piran, T., Loeb, A., Goodman, J., and Bahcall, J.N. 1987, Science, vol. 237, pp. 1471–1473. [3] 14.8 Savitzky-Golay Smoothing Filters In§13.5 we learned something about the construction and application of digital filters, but little guidance was given on which particular filter to use. That, of course, depends on what you want to accomplish by filtering. One obvious use for low-pass filters is to smooth noisy data. The premise of data smoothing is that one is measuring a variable that is both slowly varying and also corrupted by random noise. Then it can sometimes be useful to replaceeach data point by some kind of local average of surrounding data points. Since nearbypointsmeasurevery nearlythesameunderlying value, averagingcanreduce thelevelofnoisewithout (much) biasing the value obtained. We must comment editorially that the smoothing of data lies in a murky area, beyond the fringe of some better posed, and therefore more highly recommended, techniques that arediscussed elsewhere in this book. If you are fitting data to a parametric model, for example (see Chapter 15), it is almost always better to use raw data than to use data that has beenpre-processed by a smoothing procedure. Another alternative to blind smoothing is so-called“optimal”or Wiener filtering, as discussed in §13.3 and more generally in §13.6. Data smoothing isprobably most justi fiedwhen itisused simply asa graphical technique, to guide theeyethrough aforestofdatapointsallwithlargeerrorbars;orasameans ofmaking initialroughestimates of simple parameters from a graph. In this section we discuss a particular type of low-pass filter, well-adapted for data smoothing, and termed variously Savitzky-Golay [1],least-squares [2],o rDISPO(Digital Smoothing Polynomial) [3]filters. Rather than having their properties de fined in the Fourier domain, and then translated to the time domain, Savitzky-Golay filters derive directly from a particular formulation of the data smoothing problem in the time domain, as we will nowsee. Savitzky-Golay filterswereinitially(andarestilloften)usedtorendervisibletherelative widths and heights of spectral lines in noisy spectrometric data. Recall thata digital filteris applied to a series of equally spaced data values f i≡f(ti), where ti≡t0+i∆for some constant sample spacing ∆and i=...−2,−1,0,1,2,.... Wehave seen ( §13.5) that the simplesttype of digital filter(thenonrecursive or finiteimpulse response filter) replaces each data value fiby a linear combination giof itself and some number of nearby neighbors, gi=nR/summationdisplay n=−nLcnfi+n (14.8.1 ) Here nLis the number of points used “to the left ”of a data point i, i.e., earlier than it, while nRis the number used to the right, i.e.,later. A so-called causalfilter would have nR=0. As a starting point for understanding Savitzky-Golay filters, consider the simplest possible averaging procedure: For some fixed nL=nR, compute each gias the average of the data points from fi−nLtofi+nR. This is sometimes called moving window averaging andcorresponds toequation(14.8.1)withconstant cn=1 /(nL+nR+1 ). Iftheunderlying function is constant, or is changing linearly with time (increasing or decreasing), then nobias is introduced into the result. Higher points at one end of the averaging interval are on