f15-1

15.1LeastSquaresas aMaximum LikelihoodEstimator 651Sample 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).should provide (i) parameters, (ii) error estimates on the parameters, and (iii) a statistical measure of goodness-of-fit. When the third item suggests that the modelis an unlikely match to the data, then items (i) and (ii) are probably worthless. Unfortunately, many practitioners of parameter estimation never proceed beyond item(i). Theydeemafitacceptableifagraphofdataandmodel“looksgood.” Thisapproachis knownas chi-by-eye . Luckily,its practitionersget what theydeserve. CITED REFERENCES AND FURTHER READING: Bevington, P.R. 1969, Data Reduction and Error Analysis for the Physical Sciences (New York: McGraw-Hill). Brownlee, K.A. 1965, Statistical Theory and Methodology , 2nd ed. (New York: Wiley). Martin, B.R. 1971, Statistics for Physicists (New York: Academic Press). von Mises, R. 1964, Mathematical Theory of Probability and Statistics (New York: Academic Press), Chapter X. Korn, G.A., andKorn, T.M. 1968, Mathematical Handbookfor Scientists andEngineers , 2nded. (New York: McGraw-Hill), Chapters 18–19. 15.1 Least Squares as a Maximum Likelihood Estimator Supposethatwearefitting Ndatapoints (xi,yi)i=1,...,N,toamodelthat hasMadjustable parameters aj,j=1,...,M. The model predicts a functional relationship between the measured independentand dependent variables, y(x)=y(x;a1...a M)( 15.1.1 ) wherethedependenceontheparametersisindicatedexplicitlyontheright-handside. What, exactly, do we want to minimize to get fitted values for the aj’s? The first thing that comes to mind is the familiar least-squares fit, minimizeover a1...a M:N/summationdisplay i=1[yi−y(xi;a1...a M)]2(15.1.2 ) Butwheredoesthiscomefrom? Whatgeneralprinciplesisitbasedon? Theanswer to these questions takes us into the subject of maximum likelihoodestimators . Given a particular data set of xi’s and yi’s, we have the intuitive feeling that some parameter sets a1...a Mare very unlikely — those for which the model function y(x)looksnothinglike thedata—whileothersmaybeverylikely—those thatcloselyresemblethedata. Howcanwequantifythisintuitivefeeling? Howcan we select fitted parameters that are “most likely” to be correct? It is not meaningfulto ask thequestion,“What is the probabilitythat a particularset offitted parameters a 1...a Mis correct?” The reason is that there is no statistical universe of models fromwhich the parametersare drawn. Thereis just one model, the correct one, and a statistical universe of data sets that are drawn from it! 652 Chapter15. ModelingofDataSample 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).Thatbeingthecase,wecan,however,turnthequestionaround,andask,“ Given a particular set of parameters , what is the probability that this data set could have occurred?” If the yi’s take on continuous values, the probability will always be zero unless we add the phrase,“...plusor minus some fixed ∆yoneach data point.” So let’s always take this phrase as understood. If the probability of obtaining thedata set is infinitesimally small, then we can conclude that the parameters under consideration are “unlikely” to be right. Conversely, our intuition tells us that the data set should not be too improbablefor the correct choice of parameters. In other words, we identify the probability of the data given the parameters (whichis amathematicallycomputablenumber),as the likelihood oftheparameters given the data. This identification is entirely based on intuition. It has no formal mathematical basis in and of itself; as we already remarked, statistics is nota branch of mathematics! Once we make this intuitive identification, however, it is only a small further step to decide to fit for the parameters a 1...a Mprecisely by finding those values thatmaximize the likelihood defined in the above way. This form of parameter estimation is maximum likelihood estimation . We are now ready to make the connection to (15.1.2). Suppose that each data point yihas a measurement error that is independentlyrandom and distributed as a normal (Gaussian) distributionaroundthe “true” model y(x). And suppose that the standarddeviations σofthese normaldistributionsare thesame forall points. Then the probabilityof the data set is the productof the probabilitiesof each point, P∝N/productdisplay i=1/braceleftBigg exp/bracketleftBigg −1 2/parenleftbiggyi−y(xi) σ/parenrightbigg2/bracketrightBigg ∆y/bracerightBigg (15.1.3 ) Notice that there is a factor ∆yin each term in the product. Maximizing(15.1.3)is equivalentto maximizingits logarithm,or minimizingthe negativeof its logarithm, namely, /bracketleftBiggN/summationdisplay i=1[yi−y(xi)]2 2σ2/bracketrightBigg −Nlog ∆ y (15.1.4 ) Since N,σ, and ∆yare all constants, minimizing this equation is equivalent to minimizing (15.1.2). What we see is that least-squares fitting isa maximum likelihood estimation of the fitted parameters ifthe measurement errors are independent and normally distributed with constant standard deviation. Notice that we made no assumption about the linearity or nonlinearity of the model y(x;a1...)in its parameters a1...a M. Just below,we will relax ourassumptionof constant standarddeviations and obtain the very similar formulas for what is called “chi-square fitting” or “weighted least-squares fitting.” First, however, let us discuss further our very stringent assumption of a normal distribution. For a hundred years or so, mathematical statisticians have been in love with the fact that the probability distribution of the sum of a very large number of very small random deviations almost always converges to a normal distribution. (For precise statements of this central limit theorem , consult [1]or other standard works 15.1LeastSquaresas aMaximum LikelihoodEstimator 653Sample 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).on mathematical statistics.) This infatuation tended to focus interest away from the fact that, for real data, the normal distribution is often rather poorly realized, if it isrealized at all. We are often taught, rather casually, that, on average, measurements will fall within ±σof the true value 68 percent of the time, within ±2σ95 percent of the time, and within ±3σ99.7 percent of the time. Extending this, one would expect a measurement to be off by ±20σonly one time out of 2×10 88. We all know that “glitches” are much more likely than that! In some instances, the deviations from a normal distribution are easy to understand and quantify. For example, in measurements obtained by counting events, the measurement errors are usually distributed as a Poisson distribution,whose cumulative probability function was already discussed in §6.2. When the numberofcountsgoingintoonedatapointislarge,thePoissondistributionconverges towards a Gaussian. However, the convergence is not uniform when measured infractionalaccuracy. The morestandarddeviationsout on the tail of the distribution, the larger the number of counts must be before a value close to the Gaussian is realized. Thesignoftheeffectisalwaysthesame: TheGaussianpredictsthat“tail”events are much less likely than they actually (by Poisson) are. This causes such events, when theyoccur,to skew a least-squares fit much morethan they ought. Other times, the deviations from a normal distribution are not so easy to understand in detail. Experimental points are occasionally just way off. Perhaps thepowerflickeredduringapoint’smeasurement,orsomeonekickedtheapparatus,or someone wrote down a wrong number. Points like this are called outliers. They can easily turn a least-squares fit on otherwise adequate data into nonsense. Their probability of occurrence in the assumed Gaussian model is so small that themaximum likelihood estimator is willing to distort the whole curve to try to bring them, mistakenly, into line. The subject of robust statistics deals with cases where the normal or Gaussian modelisabadapproximation,orcaseswhereoutliersareimportant. Wewilldiscuss robust methods briefly in §15.7. All the sections between this one and that one assume, one way or the other, a Gaussian model for the measurement errors in the data. It it quite important that you keep the limitations of that model in mind, even as you use the very useful methods that follow from assuming it. Finally, note that our discussion of measurement errors has been limited to statistical errors, the kind that will average away if we only take enough data. Measurements are also susceptible to systematic errors that will not go away with any amount of averaging. For example, the calibration of a metal meter stick might depend on its temperature. If we take all our measurements at the same wrongtemperature, then no amount of averaging or numerical processing will correct for this unrecognized systematic error. Chi-Square Fitting We considered the chi-square statistic once before, in §14.3. Here it arises in a slightly different context. If each data point (xi,y i)has its own, known standard deviation σi, then equation (15.1.3) is modified only by putting a subscript ion the symbol σ. That subscript also propagates docilely into (15.1.4), so that the maximum likelihood 654 Chapter15. ModelingofDataSample 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).estimate of the model parameters is obtained by minimizing the quantity χ2≡N/summationdisplay i=1/parenleftbiggyi−y(xi;a1...a M) σi/parenrightbigg2 (15.1.5 ) called the “chi-square.” Towhateverextentthemeasurementerrorsactually arenormallydistributed,the quantity χ2iscorrespondinglyasumof Nsquaresofnormallydistributedquantities, eachnormalizedto unitvariance. Oncewe haveadjustedthe a1...a Mto minimize thevalueof χ2,thetermsinthesumarenotallstatisticallyindependent. Formodels that are linear in the a’s, however, it turns out that the probability distribution for different values of χ2at its minimum can nevertheless be derived analytically, and is thechi-square distribution for N−Mdegrees of freedom . We learned how to compute this probability function using the incomplete gamma function gammqin §6.2. In particular, equation (6.2.18) gives the probability Qthat the chi-square shouldexceeda particularvalue χ2by chance,where ν=N−Mis thenumberof degrees of freedom . The quantity Q, or its complement P≡1−Q, is frequently tabulated in appendices to statistics books, but we generally find it easier to use gammqandcomputeourownvalues: Q=gammq (0.5ν,0.5χ2). It is quitecommon, and usually not too wrong, to assume that the chi-square distributionholds even formodels that are not strictly linear in the a’s. This computed probabilitygives a quantitative measure for the goodness-of-fit of the model. If Qis a very small probability for some particular data set, then the apparent discrepancies are unlikely to be chance fluctuations. Much more probably either(i)themodeliswrong—canbestatisticallyrejected,or(ii)someonehasliedtoyouaboutthesizeofthemeasurementerrors σ i—theyarereallylargerthanstated. It is an important point that the chi-square probability Qdoes not directly measure the credibility of the assumption that the measurement errors are normallydistributed. It assumes they are. In most, but not all, cases, however, the effect of nonnormal errors is to create an abundance of outlier points. These decrease the probability Q,sothatwecanaddanotherpossible,thoughlessdefinitive,conclusion to the abovelist: (iii) the measurementerrors may not be normallydistributed. Possibility (iii) is fairly common, and also fairly benign. It is for this reason that reasonable experimenters are often rather tolerant of low probabilities Q.I ti s notuncommontodeemacceptableonequaltermsanymodelswith,say, Q> 0.001. This is not as sloppy as it sounds: Truly wrongmodels will often be rejected with vastly smaller values of Q,10 −18, say. However, if day-in and day-out you find yourselfacceptingmodelswith Q∼10−3, youreallyshouldtrackdownthe cause. If you happen to know the actual distribution law of your measurement errors, thenyoumightwishto MonteCarlosimulate somedatasetsdrawnfromaparticular model, cf. §7.2–§7.3. You can then subject these synthetic data sets to your actual fitting procedure, so as to determine both the probability distribution of the χ2 statistic, and also the accuracy with which your model parameters are reproduced by the fit. We discuss this further in §15.6. The technique is very general, but it can also be very expensive. Attheoppositeextreme,itsometimeshappensthattheprobability Qistoolarge, too near to 1, literally too good to be true! Nonnormal measurement errors cannot in general producethis disease, since the normal distribution is about as “compact” 15.2FittingDatatoaStraightLine 655Sample 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).as a distribution can be. Almost always, the cause of too good a chi-square fit is that the experimenter, in a “fit” of conservativism, has overestimated his or her measurement errors. Very rarely, too good a chi-square signals actual fraud, data that has been “fudged” to fit the model. A rule of thumb is that a “typical” value of χ2for a “moderately” good fit is χ2≈ν. Morepreciseisthestatementthatthe χ2statistichasamean νandastandard deviation√ 2ν, and, asymptoticallyfor large ν, becomes normallydistributed. In some cases the uncertainties associated with a set of measurements are not knownin advance,andconsiderationsrelatedto χ2fitting areused to derivea value forσ. Ifweassumethatallmeasurementshavethesamestandarddeviation, σi=σ, and that the model does fit well, then we can proceed by first assigning an arbitrary constant σto all points, next fitting for the model parameters by minimizing χ2, and finally recomputing σ2=N/summationdisplay i=1[yi−y(xi)]2/(N−M)( 15.1.6 ) Obviously, this approach prohibits an independent assessment of goodness-of-fit, a fact occasionally missed by its adherents. When, however, the measurement error is not known, this approach at least allows somekind of error bar to be assigned to the points. Ifwe takethe derivativeofequation(15.1.5)withrespectto theparameters ak, we obtain equations that must hold at the chi-square minimum, 0=N/summationdisplay i=1/parenleftbiggyi−y(xi) σ2 i/parenrightbigg/parenleftbigg∂y(xi;...a k...) ∂a k/parenrightbigg k=1,...,M (15.1.7 ) Equation(15.1.7)is, in general,a set of Mnonlinearequationsfor the Munknown ak. Various of the procedures described subsequently in this chapter derive from (15.1.7) and its specializations. CITED REFERENCES AND FURTHER READING: Bevington, P.R. 1969, Data Reduction and Error Analysis for the Physical Sciences (New York: McGraw-Hill), Chapters 1–4. von Mises, R. 1964, Mathematical Theory of Probability and Statistics (New York: Academic Press), §VI.C. [1] 15.2 Fitting Data to a Straight Line A concrete example will make the considerationsof the previous section more meaningful. We consider the problem of fitting a set of Ndata points (xi,y i)to a straight-line model y(x)=y(x;a, b)=a+bx (15.2.1 )