f15-2

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 ) 656 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).This problem is often called linear regression , a terminology that originated, long ago, in the social sciences. We assume that the uncertainty σiassociated with each measurement yiis known, and that the xi’s (values of the dependent variable) are known exactly. To measure how well the model agrees with the data, we use the chi-square merit function (15.1.5), which in this case is χ2(a, b)=N/summationdisplay i=1/parenleftbiggyi−a−bx i σi/parenrightbigg2 (15.2.2 ) Ifthemeasurementerrorsarenormallydistributed,thenthismeritfunctionwillgive maximumlikelihoodparameterestimationsof aandb;iftheerrorsarenotnormally distributed,thentheestimationsarenotmaximumlikelihood,butmaystillbeuseful in a practical sense. In §15.7, we will treat the case where outlier points are so numerous as to render the χ2merit function useless. Equation (15.2.2) is minimized to determine aandb. At its minimum, derivatives of χ2(a, b)with respect to a, bvanish. 0=∂χ2 ∂a=−2N/summationdisplay i=1yi−a−bx i σ2 i 0=∂χ2 ∂b=−2N/summationdisplay i=1xi(yi−a−bx i) σ2 i(15.2.3 ) These conditions can be rewritten in a convenient form if we define the following sums: S≡N/summationdisplay i=11 σ2 iSx≡N/summationdisplay i=1xi σ2 iSy≡N/summationdisplay i=1yi σ2 i Sxx≡N/summationdisplay i=1x2 i σ2 iSxy≡N/summationdisplay i=1xiyi σ2 i(15.2.4 ) With these definitions (15.2.3) becomes aS+bS x=Sy aS x+bS xx=Sxy(15.2.5 ) The solution of these two equations in two unknowns is calculated as ∆≡SS xx−(Sx)2 a=SxxSy−SxSxy ∆ b=SS xy−SxSy ∆(15.2.6 ) Equation(15.2.6)gives the solution for the best-fit model parameters aandb. 15.2FittingDatatoaStraightLine 657Sample 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).We are not done, however. We must estimate the probable uncertainties in the estimates of aandb, since obviously the measurement errors in the data must introduce some uncertainty in the determination of those parameters. If the data are independent, then each contributes its own bit of uncertainty to the parameters. Consideration of propagation of errors shows that the variance σ2 fin the value of any function will be σ2 f=N/summationdisplay i=1σ2 i/parenleftbigg∂f ∂y i/parenrightbigg2 (15.2.7 ) For the straight line, the derivatives of aandbwith respect to yican be directly evaluated from the solution: ∂a ∂y i=Sxx−Sxxi σ2 i∆ ∂b ∂y i=Sx i−Sx σ2 i∆(15.2.8 ) Summing over the points as in (15.2.7), we get σ2 a=Sxx/∆ σ2 b=S/∆(15.2.9 ) which are the variances in the estimates of aandb, respectively. We will see in §15.6that an additionalnumberis also neededtocharacterizeproperlythe probable uncertainty of the parameter estimation. That number is the covariance ofaandb, and (as we will see below) is given by Cov(a, b)=−Sx/∆( 15.2.10 ) The coefficient of correlation between the uncertainty in aand the uncertainty inb, which is a number between −1and 1, follows from (15.2.10) (compare equation 14.5.1), rab=−Sx√SS xx(15.2.11 ) A positive value of rabindicates that the errors in aandbare likely to have the same sign, while a negative value indicates the errors are anticorrelated, likely to have opposite signs. We arestillnot done. We must estimate the goodness-of-fit of the data to the model. Absentthisestimate,wehavenottheslightestindicationthattheparameters aandbin the model have any meaning at all! The probability Qthat a value of chi-square as pooras the value (15.2.2) should occur by chance is Q=gammq/parenleftbiggN−2 2,χ2 2/parenrightbigg (15.2.12 ) 658 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).Here gammqis our routine for the incomplete gamma function Q(a, x),§6.2. If Qis larger than, say, 0.1, then the goodness-of-fit is believable. If it is larger than, say, 0.001, then the fit maybe acceptable if the errors are nonnormal or have been moderately underestimated. If Qis less than 0.001then the model and/or estimation procedure can rightly be called into question. In this latter case, turnto§15.7 to proceed further. If you do not know the individualmeasurementerrorsof the points σ i, and are proceeding (dangerously) to use equation (15.1.6) for estimating these errors, then here is the procedure for estimating the probable uncertainties of the parameters a andb: Set σi≡1in all equations through (15.2.6), and multiply σaandσb,a s obtained from equation (15.2.9), by the additional factor/radicalbig χ2/(N−2), where χ2 is computed by (15.2.2) using the fitted parameters aandb. As discussed above, this procedure is equivalent to assuming a good fit, so you get no independent goodness-of-fit probability Q. In§14.5 we promised a relation between the linear correlation coefficient r(equation 14.5.1) and a goodness-of-fit measure, χ2(equation 15.2.2). For unweighted data (all σi=1), that relation is χ2=( 1−r2)NVar (y1...y N)( 15.2.13 ) where NVar (y1...y N)≡N/summationdisplay i=1(yi−y)2(15.2.14 ) For data with varying weights σi, the above equations remain valid if the sums in equation (14.5.1) are weighted by 1/σ2 i. The following subroutine, fit, carries out exactly the operations that we have discussed. When the weights σare known in advance, the calculations exactly correspond to the formulas above. However, when weights σare unavailable, the routine assumesequal values of σfor each point and assumesa good fit, as discussed in §15.1. The formulas (15.2.6) are susceptible to roundoff error. Accordingly, we rewrite them as follows: Define ti=1 σi/parenleftbigg xi−Sx S/parenrightbigg ,i =1,2,...,N (15.2.15 ) and Stt=N/summationdisplay i=1t2 i (15.2.16 ) Then, as you can verify by direct substitution, b=1 SttN/summationdisplay i=1tiyi σi(15.2.17 ) a=Sy−Sxb S(15.2.18 ) 15.2FittingDatatoaStraightLine 659Sample 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).σ2 a=1 S/parenleftbigg 1+S2 x SS tt/parenrightbigg (15.2.19 ) σ2 b=1 Stt(15.2.20 ) Cov(a, b)=−Sx SS tt(15.2.21 ) rab=Cov(a, b) σaσb(15.2.22 ) SUBROUTINE fit(x,y,ndata,sig,mwt,a,b,siga,sigb,chi2,q) INTEGER mwt,ndataREAL a,b,chi2,q,siga,sigb,sig(ndata),x(ndata),y(ndata) C USES gammq Given a set of data points x(1:ndata) ,y(1:ndata) with individual standard deviations sig(1:ndata) , fit them to a straight line y=a+bxby minimizing χ2. Returned are a,b and their respective probable uncertainties siga andsigb , the chi-square chi2 ,a n d the goodness-of-fit probability q(that the fit would have χ2this large or larger). If mwt=0 on input, then the standard deviations are assumed to be unavailable: qis returned as 1.0 and the normalization of chi2 is to unit standard deviation on all points. INTEGER i REAL sigdat,ss,st2,sx,sxoss,sy,t,wt,gammqsx=0. Initialize sums to zero. sy=0. st2=0. b=0.if(mwt.ne.0) then Accumulate sums ... ss=0. do 11i=1,ndata ...with weights wt=1./(sig(i)**2)ss=ss+wt sx=sx+x(i)*wt sy=sy+y(i)*wt enddo 11 else do12i=1,ndata ...or without weights. sx=sx+x(i)sy=sy+y(i) enddo 12 ss=float(ndata) endifsxoss=sx/ss if(mwt.ne.0) then do 13i=1,ndata t=(x(i)-sxoss)/sig(i) st2=st2+t*t b=b+t*y(i)/sig(i) enddo 13 else do14i=1,ndata t=x(i)-sxossst2=st2+t*t b=b+t*y(i) enddo 14 endifb=b/st2 Solve for a,b,σ a,a n d σb. a=(sy-sx*b)/ss siga=sqrt((1.+sx*sx/(ss*st2))/ss)sigb=sqrt(1./st2) chi2=0. Calculate χ 2. 660 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).q=1. if(mwt.eq.0) then do15i=1,ndata chi2=chi2+(y(i)-a-b*x(i))**2 enddo 15 sigdat=sqrt(chi2/(ndata-2)) For unweighted data evaluate typical sig us- ingchi2 , and adjust the standard devia- tions.siga=siga*sigdat sigb=sigb*sigdat else do16i=1,ndata chi2=chi2+((y(i)-a-b*x(i))/sig(i))**2 enddo 16 if(ndata.gt.2) q=gammq(0.5*(ndata-2),0.5*chi2) Equation (15.2.12). endif returnEND CITED REFERENCES AND FURTHER READING: Bevington, P.R. 1969, Data Reduction and Error Analysis for the Physical Sciences (New York: McGraw-Hill), Chapter 6. 15.3 Straight-Line Data with Errors in Both Coordinates If experimental data are subject to measurement error not only in the yi’s, but also in thexi’s, then the task of fitting a straight-line model y(x)=a+bx (15.3.1 ) is considerably harder. Itis straightforward to write down the χ2meritfunction for this case, χ2(a, b )=N/summationdisplay i=1(yi−a−bx i)2 σ2 yi+b2σ2 xi(15.3.2 ) where σxiandσyiare, respectively, the xandystandard deviations for the ith point. The weighted sum of variances in the denominator of equation (15.3.2) can be understood bothas the variance in the direction of the smallest χ 2between each data point and the line with slope b, and also as the variance of the linear combination yi−a−bx iof two random variables xiandyi, Var (yi−a−bx i)=Var (yi)+b2Var (xi)=σ2 yi+b2σ2 xi≡1/w i (15.3.3 ) The sum of the square of Nrandom variables, each normalized by its variance, is thus χ2-distributed. We want to minimize equation (15.3.2) with respect to aandb. Unfortunately, the occurrence of bin the denominator of equation (15.3.2) makes the resulting equation for the slope ∂χ2/∂b =0nonlinear. However, the corresponding condition for the intercept, ∂χ2/∂a =0, is still linear and yields a=/bracketleftBigg/summationdisplay iwi(yi−bx i)/bracketrightBigg/slashBigg/summationdisplay iwi (15.3.4 ) where the wi’s are defined by equation (15.3.3). A reasonable strategy, now, is to use the machinery of Chapter 10 (e.g., the routine brent) for minimizing a general one-dimensional