f15-3

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−bxi)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−bxiof two random variables xiandyi, Var(yi−a−bxi)=Var(yi)+b2Var(xi)=σ2 yi+b2σ2 xi≡1/wi (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−bxi)/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 15.3Straight-LineDatawithErrorsinBothCoordinates 661Sample 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 = 1σaAB σb 0b as r Figure 15.3.1. Standard errors for the parameters aandb. The point Bcan be found by varying the slope bwhile simultaneously minimizing the intercept a. This gives the standard error σb, and also the value s. The standard error σacan then be found by the geometric relation σ2 a=s2+r2. function to minimize with respect to b, while using equation (15.3.4) at each stage to ensure that the minimum with respect to bis also minimized with respect to a. Because of the finite error bars on the xi’s, the minimum χ2as a function of bwill befinite, though usually large, when bequals infinity (line of in finite slope). The angle θ≡arctan bisthus moresuitableasaparametrization ofslope than bitself. Thevalue of χ2 will then be periodic in θwith period π(not2π!). If any data points have very small σy’s but moderate or large σx’s, then it is also possible to have a maximum in χ2near zero slope, θ≈0. In that case, there can sometimes be two χ2minima, one at positive slope and the other at negative. Only one of these is the correct global minimum. It is therefore importantto have a good starting guess for b(orθ). Our strategy, implemented below, is to scale the y i’s so as to have variance equal to the xi’s, then to do a conventional (as in §15.2) linear fit with weights derived from the (scaled) sum σ2 yi+σ2 xi. This yields a good starting guess for bif the data are even plausibly related to a straight-line model. Finding the standard errors σaandσbon the parameters aandbis more complicated. We will see in §15.6 that, in appropriate circumstances, the standard errors in aandbare the respective projections onto the aandbaxes of the “confidence region boundary ”where χ2 takes on a value one greater than its minimum, ∆χ2=1. In the linear case of §15.2, these projections follow from the Taylor series expansion ∆χ2≈1 2/bracketleftbigg∂2χ2 ∂a2(∆a)2+∂2χ2 ∂b2(∆b)2/bracketrightbigg +∂2χ2 ∂a∂b∆a∆b (15.3.5 ) Becauseofthepresentnonlinearityin b,however, analyticformulasforthesecondderivatives arequiteunwieldy; moreimportant,thelowest-ordertermfrequently givesapoorapproxima-tion to ∆χ 2. Our strategy is therefore to find the roots of ∆χ2=1numerically, by adjusting thevalueoftheslope bawayfromtheminimum. Intheprogrambelowthegeneralroot finder zbrentis used. Itmay occur that there are no roots at all —for example, if all error bars are so large that all the data points are compatible with each other. It is important, therefore, tomake some effort at bracketing a putative root before re fining it (cf. §9.1). Because ais minimized at each stage of varying b, successful numerical root- finding leads to a value of ∆athat minimizes χ 2for the value of ∆bthat gives ∆χ2=1. This (see Figure 15.3.1) directly gives the tangent projection of the con fidence region onto the baxis, 662 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).and thus σb. It does not, however, give the tangent projection of the con fidence region onto theaaxis. In the figure, we have found the point labeled B;t ofindσawe need to find the point A. Geometry to the rescue: To the extent that the con fidence region is approximated by an ellipse, then you can prove (see figure) that σ2 a=r2+s2. The value of sis known from having found the point B. The value of rfollows from equations (15.3.2) and (15.3.3) applied at the χ2minimum (point Oin thefigure), giving r2=1/slashBigg/summationdisplay iwi (15.3.6 ) Actually, since bcan go through in finity, this whole procedure makes more sense in (a, θ)space than in (a, b)space. That is in fact how the following program works. Since it is conventional, however, to return standard errors for aandb, not aandθ,w efinally use the relation σb=σθ/cos2θ (15.3.7 ) Wecautionthatif banditsstandarderrorarebothlarge,sothatthecon fidence regionactually includes in finiteslope,then thestandard error σbisnot verymeaningful. Thefunction chixy isnormallycalledonlybytheroutine fitexy. However,ifyouwant,youcanyourselfexplore the confidence region by making repeated calls to chixy(whose argument is an angle θ, not a slope b), after a single initializing call to fitexy. Afinal caution, repeated from §15.0, is that if the goodness-of- fit is not acceptable (returned probability is too small), the standard errors σaandσbare surely not believable. In direcircumstances, you might tryscaling allyour xandyerror barsby a constant factor until the probability is acceptable (0.5, say), to get more plausible values for σaandσb. SUBROUTINE fitexy(x,y,ndat,sigx,sigy,a,b,siga,sigb,chi2,q) INTEGER ndat,NMAX REAL x(ndat),y(ndat),sigx(ndat),sigy(ndat),a,b,siga,sigb,chi2, * q,POTN,PI,BIG,ACC PARAMETER (NMAX=1000,POTN=1.571000,BIG=1.e30,PI=3.14159265, * ACC=1.e-3) C USES avevar,brent,chixy,fit,gammq,mnbrak,zbrent Straight-line fit to input data x(1:ndat) andy(1:ndat) with errors in both xandy,t h e respective standard deviations being the input quantities sigx(1:ndat) andsigy(1:ndat) . Output quantities are aandbsuch that y=a+bxminimizes χ2, whose value is returned aschi2 .T h e χ2probability is returned as q, a small value indicating a poor fit (sometimes indicating underestimated errors). Standard errors on aandbare returned as siga and sigb . These are not meaningful if either (i) the fit is poor, or (ii) bis so large that the data are consistent with a vertical (infinite b) line. If siga andsigb are returned as BIG, then the data are consistent with allvalues of b. INTEGER j,nnREAL xx(NMAX),yy(NMAX),sx(NMAX),sy(NMAX),ww(NMAX),swap,amx,amn * ,varx,vary,aa,offs,ang(6),ch(6),scale,bmn,bmx,d1,d2 * ,r2,dum1,dum2,dum3,dum4,dum5,brent,chixy,gammq,zbrent COMMON /fitxyc/ xx,yy,sx,sy,ww,aa,offs,nn EXTERNAL chixy if (ndat.gt.NMAX) pause ’NMAX too small in fitexy’call avevar(x,ndat,dum1,varx) Find the xandyvariances, and scale the data into the common block for communication with the function chixy .call avevar(y,ndat,dum1,vary) scale=sqrt(varx/vary) nn=ndatdo 11j=1,ndat xx(j)=x(j) yy(j)=y(j)*scalesx(j)=sigx(j)sy(j)=sigy(j)*scale ww(j)=sqrt(sx(j)**2+sy(j)**2) Use both xandyweights in first trial fit. enddo 11 call fit(xx,yy,nn,ww,1,dum1,b,dum2,dum3,dum4,dum5) Trial fit for b. offs=0. 15.3Straight-LineDatawithErrors inBothCoordinates 663Sample 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).ang(1)=0. Construct several angles for reference points. ang(2)=atan(b) Make ban angle. ang(4)=0. ang(5)=ang(2)ang(6)=POTN do 12j=4,6 ch(j)=chixy(ang(j)) enddo 12 call mnbrak(ang(1),ang(2),ang(3),ch(1),ch(2),ch(3),chixy) Bracket the χ2min- imum and then locate it with brent . chi2=brent(ang(1),ang(2),ang(3),chixy,ACC,b) chi2=chixy(b)a=aaq=gammq(0.5*(nn-2),0.5*chi2) Compute χ 2probability. r2=0. do13j=1,nn Save the inverse sum of weights at the mini- mum. r2=r2+ww(j) enddo 13 r2=1./r2 bmx=BIG Now, find standard errors for bas points where ∆χ2=1. bmn=BIG offs=chi2+1. do14j=1,6 Go through saved values to bracket the desired roots. Note periodicity in slope angles. if (ch(j).gt.offs) then d1=mod(abs(ang(j)-b),PI) d2=PI-d1if(ang(j).lt.b)then swap=d1 d1=d2 d2=swap endif if (d1.lt.bmx) bmx=d1 if (d2.lt.bmn) bmn=d2 endif enddo 14 if (bmx.lt. BIG) then Call zbrent to find the roots. bmx=zbrent(chixy,b,b+bmx,ACC)-bamx=aa-abmn=zbrent(chixy,b,b-bmn,ACC)-b amn=aa-a sigb=sqrt(0.5*(bmx**2+bmn**2))/(scale*cos(b)**2)siga=sqrt(0.5*(amx**2+amn**2)+r2)/scale Error in ahas additional piece r2. else sigb=BIGsiga=BIG endif a=a/scale Unscale the answers. b=tan(b)/scalereturn END FUNCTION chixy(bang) REAL chixy,bang,BIG INTEGER NMAXPARAMETER (NMAX=1000,BIG=1.E30) Captive function of fitexy , returns the value of (χ2−offs)for the slope b=tan(bang) . Scaled data and offs are communicated via the common block /fitxyc/ . INTEGER nn,jREAL xx(NMAX),yy(NMAX),sx(NMAX),sy(NMAX),ww(NMAX),aa,offs, * avex,avey,sumw,b COMMON /fitxyc/ xx,yy,sx,sy,ww,aa,offs,nnb=tan(bang) avex=0. 664 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).avey=0. sumw=0. do11j=1,nn ww(j)=(b*sx(j))**2+sy(j)**2if(ww(j).lt.1./BIG) then ww(j)=BIG else ww(j)=1./ww(j) endif sumw=sumw+ww(j) avex=avex+ww(j)*xx(j)avey=avey+ww(j)*yy(j) enddo 11 avex=avex/sumwavey=avey/sumwaa=avey-b*avex chixy=-offs do 12j=1,nn chixy=chixy+ww(j)*(yy(j)-aa-b*xx(j))**2 enddo 12 returnEND Be aware that the literature on the seemingly straightforward subject of this section is generally confusing and sometimes plain wrong. Deming ’s[1]early treatment is sound, but its reliance on Taylor expansions gives inaccurate error estimates. References [2-4]are reliable, more recent, general treatments with critiques of earlier work. York [5]and Reed [6] usefully discuss the simple case of a straight line as treated here, but the latter paper hassome errors, corrected in [7]. All this commotion has attracted the Bayesians [8-10], who have still different points of view. CITED REFERENCES AND FURTHER READING: Deming,W.E.1943, StatisticalAdjustmentofData (NewYork:Wiley),reprinted1964(NewYork: Dover). [1] Jefferys, W.H. 1980, Astronomical Journal , vol. 85, pp. 177–181; see also vol. 95, p. 1299 (1988). [2] Jefferys, W.H. 1981, Astronomical Journal , vol. 86, pp. 149–155; see also vol. 95, p. 1300 (1988). [3] Lybanon, M. 1984, American Journal of Physics , vol. 52, pp. 22–26. [4] York, D. 1966, Canadian Journal of Physics , vol. 44, pp. 1079–1086. [5] Reed, B.C. 1989, American Journal of Physics , vol. 57, pp. 642–646; see also vol. 58, p. 189, and vol. 58, p. 1209. [6] Reed, B.C. 1992, American Journal of Physics , vol. 60, pp. 59–62. [7] Zellner, A. 1971, An Introduction to Bayesian Inference in Econometrics (New York: Wiley); reprinted 1987 (Malabar, FL: R. E. Krieger Pub. Co.). [8] Gull,S.F.1989,in MaximumEntropyandBayesianMethods ,J.Skilling,ed.(Boston:Kluwer).[9] Jaynes, E.T. 1991, in Maximum-Entropy and Bayesian Methods, Proc. 10th Int. Workshop , W.T. Grandy, Jr., and L.H. Schick, eds. (Boston: Kluwer). [10] Macdonald, J.R., and Thompson, W.J. 1992, American Journalof Physics , vol. 60, pp. 66–73. 15.4 GeneralLinearLeastSquares 665Sample 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).15.4 General Linear Least Squares An immediate generalization of §15.2 is to fit a set of data points (xi,y i)to a model that is not just a linear combinationof 1andx(namely a+bx), but rather a linear combination of anyMspecified functions of x. For example, the functions could be 1,x ,x2,...,xM−1, in which case their general linear combination, y(x)=a1+a2x+a3x2+··· +aMxM−1(15.4.1 ) is a polynomial of degree M−1. Or, the functions could be sines and cosines, in which case their general linear combination is a harmonic series. The general form of this kind of model is y(x)=M/summationdisplay k=1akXk(x)( 15.4.2 ) where X1(x),...,X M(x)are arbitrary fixed functions of x, called the basis functions. Note that the functions Xk(x)can be wildly nonlinear functions of x. In this discussion “linear”refers onlyto the model ’s dependenceon its parameters ak. For these linear models we generalize the discussion of the previous section by defining a merit function χ2=N/summationdisplay i=1/bracketleftBigg yi−/summationtextM k=1akXk(xi) σi/bracketrightBigg2 (15.4.3 ) As before, σiis the measurement error (standard deviation) of the ith data point, presumed to be known. If the measurement errors are not known, they may all (as discussed at the end of §15.1) be set to the constant value σ=1. Once again,we will pick as best parametersthose that minimize χ2. Thereare severaldifferenttechniquesavailablefor findingthisminimum. Twoareparticularly useful, and we will discuss both in this section. To introduce them and elucidate their relationship, we need some notation. LetAbe a matrix whose N×Mcomponents are constructed from the M basis functionsevaluatedat the Nabscissas xi, andfromthe Nmeasurementerrors σi, by the prescription Aij=Xj(xi) σi(15.4.4 ) Thematrix Aiscalledthe designmatrix ofthefittingproblem. Noticethatingeneral Ahas more rows than columns, N≥M, since there must be more data points than modelparameterstobesolvedfor. (Youcan fitastraightlinetotwopoints,butnota verymeaningfulquintic!) ThedesignmatrixisshownschematicallyinFigure15.4.1. Also define a vector bof length Nby bi=yi σi(15.4.5 ) and denote the Mvector whose components are the parameters to be fitted, a1,...,a M,b ya.