f15-4

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 severaldifferenttechniquesavailableforfindingthisminimum. Twoareparticularlyuseful, 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 x i, 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. (Youcanfitastraightlinetotwopoints,butnotaverymeaningfulquintic!) ThedesignmatrixisshownschematicallyinFigure15.4.1. Also define a vector bof length Nby b i=yi σi(15.4.5 ) and denote the Mvector whose components are the parameters to be fitted, a1,...,a M,b ya. 666 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).X1(x1) σ1x1 X2(x1) σ1. . .XM(x1) σ1X1()X2(). . .XM() X1(x2) σ2x2 X2(x2) σ2. . .XM(x2) σ2 ......... ...... ...... X1(xN) σNxN X2(xN) σN. . .XM(xN) σNdata pointsbasis functions Figure 15.4.1. Design matrix for the least-squares fit of a linear combination of Mbasis functions to N data points. The matrix elements involve the basis functions evaluated at the values of the independentvariableatwhichmeasurementsaremade,andthestandarddeviationsofthemeasureddependentvariable. The measured values of the dependent variable do not enter the design matrix. Solutionby Use ofthe NormalEquations The minimum of (15.4.3)occurs where the derivative of χ2with respect to all Mparameters akvanishes. Specializing equation (15.1.7) to the case of the model (15.4.2), this condition yields the Mequations 0=N/summationdisplay i=11 σ2 i yi−M/summationdisplay j=1ajXj(xi) Xk(xi) k=1,...,M (15.4.6 ) Interchangingtheorderofsummations,wecanwrite(15.4.6)asthematrixequation M/summationdisplay j=1αkjaj=βk (15.4.7 ) where αkj=N/summationdisplay i=1Xj(xi)Xk(xi) σ2 iorequivalently [α]=AT·A (15.4.8 ) anM×Mmatrix, and βk=N/summationdisplay i=1yiXk(xi) σ2 iorequivalently [β]=AT·b (15.4.9 ) 15.4 GeneralLinearLeastSquares 667Sample 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).a vector of length M. The equations (15.4.6)or (15.4.7)are called the normalequations of the least- squares problem. Theycan be solved forthe vectorof parameters aby the standard methods of Chapter 2, notably LUdecomposition and backsubstitution, Choleksy decomposition, or Gauss-Jordan elimination. In matrix form, the normal equationscan be written as either [α]·a=[β]oras/parenleftbig A T·A/parenrightbig ·a=AT·b (15.4.10 ) The inverse matrix Cjk≡[α]−1 jkis closely related to the probable (or, more precisely, standard) uncertainties of the estimated parameters a. To estimate these uncertainties, consider that aj=M/summationdisplay k=1[α]−1 jkβk=M/summationdisplay k=1Cjk/bracketleftBiggN/summationdisplay i=1yiXk(xi) σ2 i/bracketrightBigg (15.4.11 ) andthatthevarianceassociatedwiththeestimate ajcanbefoundasin(15.2.7)from σ2(aj)=N/summationdisplay i=1σ2 i/parenleftbigg∂a j ∂y i/parenrightbigg2 (15.4.12 ) Note that αjkis independent of yi, so that ∂a j ∂y i=M/summationdisplay k=1CjkXk(xi)/σ2 i (15.4.13 ) Consequently, we find that σ2(aj)=M/summationdisplay k=1M/summationdisplay l=1CjkCjl/bracketleftBiggN/summationdisplay i=1Xk(xi)Xl(xi) σ2 i/bracketrightBigg (15.4.14 ) Thefinal term in brackets is just the matrix [α]. Since this is the matrix inverse of [C], (15.4.14) reduces immediately to σ2(aj)=Cjj (15.4.15 ) In other words, the diagonal elements of [C]are the variances (squared uncertainties) of the fitted parameters a. It should not surprise you to learn that the off-diagonalelements Cjkare the covariances between ajandak(cf. 15.2.10); but we shall defer discussion of these to §15.6. We will now givea routinethat implements the aboveformulasfor the general linear least-squares problem, by the method of normal equations. Since we wish to computenotonly the solutionvector abut also the covariancematrix [C], it is most convenienttouseGauss-Jordanelimination(routine gaussjof§2.1)to performthe linearalgebra. Theoperationcount,inthisapplication,isnolargerthanthatfor LU decomposition. If you have no need for the covariance matrix, however, you can save a factor of 3 on the linear algebra by switching to LUdecomposition, without 668 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).computation of the matrix inverse. In theory, since AT·Ais positive de finite, Cholesky decomposition is the most ef ficient way to solve the normal equations. However, in practice most of the computing time is spent in looping over the data to form the equations, and Gauss-Jordan is quite adequate. We need to warn you that the solution of a least-squaresproblemdirectlyfrom the normal equations is rather susceptible to roundoff error. An alternative, and preferred, technique involves QRdecomposition ( §2.10,§11.3, and §11.6) of the designmatrix A. Thisisessentiallywhatwedidattheendof §15.2forfittingdatato astraightline,butwithoutinvokingall themachineryof QRtoderivethenecessary formulas. Later in this section, we will discuss otherdif ficulties in the least-squares problem,forwhichthecureis singularvaluedecomposition (SVD),ofwhichwegive animplementation. ItturnsoutthatSVDalso fixestheroundoffproblem,soitisour recommendedtechniqueforallbut “easy”least-squaresproblems. Itisfortheseeasy problemsthatthefollowingroutine,whichsolvesthe normalequations,is intended. The routine below introduces one bookkeeping trick that is quite useful in practical work. Frequently it is a matter of “art”to decide which parameters ak in a model should be fit from the data set, and which should be held constant at fixed values, for example values predicted by a theory or measured in a previous experiment. One wants, therefore, to have a convenient means for “freezing” and“unfreezing ”the parameters ak. In the following routine the total number of parameters akis denoted ma(called Mabove). As input to the routine, you supply an array ia(1:ma) , whose components are either zero or nonzero (e.g., 1). Zeros indicatethatyouwant thecorrespondingelementsof theparametervector a(1:ma) to be held fixed at their input values. Nonzeros indicate parameters that should be fitted for. On output, any frozen parameters will have their variances, and all their covariances, set to zero in the covariance matrix. SUBROUTINE lfit(x,y,sig,ndat,a,ia,ma,covar,npc,chisq,funcs) INTEGER ma,ia(ma),npc,ndat,MMAX REAL chisq,a(ma),covar(npc,npc),sig(ndat),x(ndat),y(ndat) EXTERNAL funcsPARAMETER (MMAX=50) Set to the maximum number of coefficients ma. C USES covsrt,gaussj Given a set of data points x(1:ndat) ,y(1:ndat) with individual standard deviations sig(1:ndat) ,u s eχ2minimization to fit for some or all of the coefficients a(1:ma) of a function that depends linearly on a,y=/summationtext iai×afunc i(x). The input array ia(1:ma) in- dicates by nonzero entries those components of athat should be fitted for, and by zero entries those components that should be held fixed at their input values. The program returns values fora(1:ma) ,χ2=chisq, and the covariance matrix covar(1:ma,1:ma) . (Parameters held fixed will return zero covariances.) npcis the physical dimension of covar(npc,npc) in the calling routine. The user supplies a subroutine funcs(x,afunc,ma) that returns themabasis functions evaluated at x=xin the array afunc. INTEGER i,j,k,l,m,mfitREAL sig2i,sum,wt,ym,afunc(MMAX),beta(MMAX) mfit=0 do 11j=1,ma if(ia(j).ne.0) mfit=mfit+1 enddo 11 if(mfit.eq.0) pause ’lfit: no parameters to be fitted’ do13j=1,mfit Initialize the (symmetric) matrix. do12k=1,mfit covar(j,k)=0. enddo 12 beta(j)=0. enddo 13 15.4 GeneralLinearLeastSquares 669Sample 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).do17i=1,ndat Loop over data to accumulate coefficients of the normal equations. call funcs(x(i),afunc,ma) ym=y(i) if(mfit.lt.ma) then Subtract off dependences on known pieces of the fitting function. do14j=1,ma if(ia(j).eq.0) ym=ym-a(j)*afunc(j) enddo 14 endifsig2i=1./sig(i)**2 j=0 do 16l=1,ma if (ia(l).ne.0) then j=j+1 wt=afunc(l)*sig2i k=0do 15m=1,l if (ia(m).ne.0) then k=k+1covar(j,k)=covar(j,k)+wt*afunc(m) endif enddo 15 beta(j)=beta(j)+ym*wt endif enddo 16 enddo 17 do19j=2,mfit Fill in above the diagonal from symmetry. do18k=1,j-1 covar(k,j)=covar(j,k) enddo 18 enddo 19 call gaussj(covar,mfit,npc,beta,1,1) Matrix solution. j=0do 21l=1,ma if(ia(l).ne.0) then j=j+1 a(l)=beta(j) Partition solution to appropriate coefficients a. endif enddo 21 chisq=0. Evaluate χ2of the fit. do23i=1,ndat call funcs(x(i),afunc,ma) sum=0. do22j=1,ma sum=sum+a(j)*afunc(j) enddo 22 chisq=chisq+((y(i)-sum)/sig(i))**2 enddo 23 call covsrt(covar,npc,ma,ia,mfit) Sort covariance matrix to true order of fitting return coefficients. END That last call to a subroutine covsrtis only for the purpose of spreading the covariances back into the full ma×macovariance matrix, in the proper rows and columns and with zero variances and covariances set for variables which were held frozen. The subroutine covsrtis as follows. SUBROUTINE covsrt(covar,npc,ma,ia,mfit) INTEGER ma,mfit,npc,ia(ma)REAL covar(npc,npc) Expand in storage the covariance matrix covar, so as to take into account parameters that are being held fixed. (For the latter, return zero covariances.) INTEGER i,j,k REAL swap 670 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).do12i=mfit+1,ma do11j=1,i covar(i,j)=0. covar(j,i)=0. enddo 11 enddo 12 k=mfitdo 15j=ma,1,-1 if(ia(j).ne.0)then do13i=1,ma swap=covar(i,k)covar(i,k)=covar(i,j)covar(i,j)=swap enddo 13 do14i=1,ma swap=covar(k,i) covar(k,i)=covar(j,i) covar(j,i)=swap enddo 14 k=k-1 endif enddo 15 return END Solutionby Use ofSingular Value Decomposition In some applications, the normal equations are perfectly adequate for linear least-squaresproblems. However,inmanycasesthenormalequationsareveryclose to singular. A zero pivot element may be encountered during the solution of thelinear equations (e.g., in gaussj), in which case you get no solution at all. Or a very small pivot may occur, in which case you typically get fitted parameters a k with verylarge magnitudesthat are delicately (andunstably)balanced to cancel out almost precisely when the fitted function is evaluated. Why does this commonly occur? The reason is that, more often than experi- menterswouldliketo admit,datadonot clearlydistinguishbetweentwo ormoreof the basis functions provided. If two such functions, or two different combinations of functions, happen to fit the data about equally well —or equally badly —then the matrix [α], unable to distinguish between them, neatly folds up its tent and becomessingular. Thereisacertainmathematicalironyinthefactthatleast-squares problems are bothoverdetermined (number of data points greater than number of parameters) andunderdetermined (ambiguous combinations of parameters exist); but that is how it frequently is. The ambiguities can be extremely hard to notice a prioriin complicated problems. Enter singularvalue decomposition(SVD). This would be a goodtime for you to review the material in §2.6, which we will not repeat here. In the case of an overdetermined system, SVD produces a solution that is the best approximation in the least-squares sense, cf. equation (2.6.10). That is exactly what we want. In the case of an underdeterminedsystem, SVD producesa solution whose values (for us,thea k’s) are smallest in the least-squares sense, cf. equation (2.6.8). That is also whatwewant: Whensomecombinationofbasisfunctionsisirrelevanttothe fit,that combination will be driven down to a small, innocuous, value, rather than pushed up to delicately canceling in finities. 15.4 GeneralLinearLeastSquares 671Sample 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).In terms of the design matrix A(equation 15.4.4) and the vector b(equation 15.4.5), minimization of χ2in (15.4.3) can be written as findathatminimizes χ2=|A·a−b|2(15.4.16 ) Comparingtoequation(2.6.9),weseethatthisispreciselytheproblemthatroutines svdcmpandsvbksbaredesignedtosolve. Thesolution,whichisgivenbyequation (2.6.12), can be rewritten as follows: If UandVenter the SVD decomposition ofAaccording to equation (2.6.1), as computed by svdcmp, then let the vectors U(i)i=1,...,Mdenote the columnsofU(each one a vector of length N); and let the vectors V(i);i=1,...,Mdenote the columns ofV(each one a vector of length M). Then the solution (2.6.12) of the least-squares problem (15.4.16) can be written as a=M/summationdisplay i=1/parenleftbiggU(i)·b wi/parenrightbigg V(i) (15.4.17 ) where the wiare, as in §2.6, the singular values returned by svdcmp. Equation (15.4.17) says that the fitted parameters aare linear combinations of thecolumnsof V,withcoef ficientsobtainedbyformingdotproductsofthecolumns ofUwith theweighteddatavector(15.4.5). Thoughit is beyondourscopeto prove here,itturnsoutthatthestandard(loosely, “probable”)errorsinthe fittedparameters are also linear combinations of the columns of V. In fact, equation (15.4.17) can be written in a form displaying these errors as a=/bracketleftBiggM/summationdisplay i=1/parenleftbiggU(i)·b wi/parenrightbigg V(i)/bracketrightBigg ±1 w1V(1)±···±1 wMV(M) (15.4.18 ) Here each ±is followed by a standard deviation. The amazing fact is that, decomposed in this fashion, the standard deviations are all mutually independent (uncorrelated). Therefore they can be added together in root-mean-square fashion.What is going on is that the vectors V (i)are the principal axes of the error ellipsoid of thefitted parameters a(see§15.6). It follows that the variance in the estimate of a parameter ajis givenby σ2(aj)=M/summationdisplay i=11 w2 i[V(i)]2 j=M/summationdisplay i=1/parenleftbiggVji wi/parenrightbigg2 (15.4.19 ) whose result should be identical with (15.4.14). As before, you should not be surprised at the formula for the covariances, here given without proof, Cov (aj,a k)=M/summationdisplay i=1/parenleftbiggVjiVki w2 i/parenrightbigg (15.4.20 ) We introduced this subsection by noting that the normal equations can fail by encountering a zero pivot. We have not yet, however, mentioned how SVD overcomes this problem. The answer is: If any singular value wiis zero, its 672 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).reciprocal in equation (15.4.18) should be set to zero, not in finity. (Compare the discussion preceding equation 2.6.7.) This corresponds to adding to the fitted parameters aazeromultiple, rather than some random large multiple, of any linear combinationofbasisfunctionsthataredegenerateinthe fit. Itisa goodthingtodo! Moreover, if a singular value wiis nonzero but very small, you should also defineitsreciprocal to be zero, since its apparent value is probably an artifact of roundoff error, not a meaningful number. A plausible answer to the question “how small is small? ”is to edit in this fashion all singular values whose ratio to the largest singular value is less than Ntimes the machine precision /epsilon1. (You might argue for√ N, or a constant, instead of Nas the multiple; that starts getting into hardware-dependent questions.) There is another reason for editing even additional singular values, ones large enough that roundoff error is not a question. Singular value decomposition allowsyou to identify linear combinations of variables that just happen not to contribute much to reducing the χ 2of your data set. Editing these can sometimes reduce the probableerroronyourcoef ficientsquitesigni ficantly,whileincreasingtheminimum χ2only negligibly. We will learn more about identifying and treating such cases in§15.6. In the following routine, the point at which this kind of editing would occur is indicated. Generallyspeaking,werecommendthatyoualwaysuseSVDtechniquesinstead ofusingthenormalequations. SVD ’sonlysigni ficantdisadvantageisthatitrequires an extra array of size N×Mto store the whole design matrix. This storage is overwritten by the matrix U. Storage is also required for the M×Mmatrix V, but this is instead of the same-sized matrix for the coef ficients of the normal equations. SVD can be signi ficantly slower than solving the normal equations; however, its great advantage, that it (theoretically) cannot fail , more than makes up for the speed disadvantage. Intheroutinethatfollows,thematrices u,vandthevector wareinputasworking space. npandmparetheirvariousphysicaldimensions.Thelogicaldimensionsofthe problemare ndatadata points by mabasis functions(and fitted parameters). If you care only about the values aof thefitted parameters, then u,v,wcontain no useful informationonoutput. Ifyouwantprobableerrorsforthe fittedparameters,readon. SUBROUTINE svdfit(x,y,sig,ndata,a,ma,u,v,w,mp,np, * chisq,funcs) INTEGER ma,mp,ndata,np,NMAX,MMAX REAL chisq,a(ma),sig(ndata),u(mp,np),v(np,np),w(np), * x(ndata),y(ndata),TOL EXTERNAL funcs PARAMETER (NMAX=1000,MMAX=50,TOL=1.e-5) NMAXis the maximum expected value of ndata;MMAXthe maximum expected for ma;t h e default TOLvalue is appropriate for single precision and variables scaled to be of order unity. C USES svbksb,svdcmp Given a set of data points x(1:ndata) ,y(1:ndata) with individual standard deviations sig(1:ndata) ,u s eχ2minimization to determine the macoefficients aof the fitting func- tiony=/summationtext iai×afunc i(x). Here we solve the fitting equations using singular value decom- position of the ndata bymamatrix, as in §2.6. Arrays u(1:mp,1:np),v(1:np,1:np), w(1:np) provide workspace on input; on output they define the singular value decomposi- tion, and can be used to obtain the covariance matrix. mp,np are the physical dimensions of the matrices u,v,w, as indicated above. It is necessary that mp≥ndata,np≥ma.T h e program returns values for the mafit parameters a,a n d χ2,chisq. The user supplies a subroutine funcs(x,afunc,ma) that returns the mabasis functions evaluated at x=x in the array afunc. INTEGER i,j 15.4 GeneralLinearLeastSquares 673Sample 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).REAL sum,thresh,tmp,wmax,afunc(MMAX),b(NMAX) do12i=1,ndata Accumulate coefficients of the fitting ma- trix. call funcs(x(i),afunc,ma) tmp=1./sig(i)do 11j=1,ma u(i,j)=afunc(j)*tmp enddo 11 b(i)=y(i)*tmp enddo 12 call svdcmp(u,ndata,ma,mp,np,w,v) Singular value decomposition. wmax=0. Edit the singular values, given TOLfrom the parameter statement, between here ... do13j=1,ma if(w(j).gt.wmax)wmax=w(j) enddo 13 thresh=TOL*wmax do14j=1,ma if(w(j).lt.thresh)w(j)=0. enddo 14 ...and here. call svbksb(u,w,v,ndata,ma,mp,np,b,a)chisq=0. Evaluate chi-square. do 16i=1,ndata call funcs(x(i),afunc,ma)sum=0. do 15j=1,ma sum=sum+a(j)*afunc(j) enddo 15 chisq=chisq+((y(i)-sum)/sig(i))**2 enddo 16 return END Feeding the matrix vand vector woutput by the above program into the following short routine, you easily obtain variances and covariances of the fitted parameters a. The square roots of the variances are the standard deviations of thefitted parameters. The routine straightforwardly implements equation (15.4.20) above, with the convention that singular values equal to zero are recognized as having been edited out of the fit. SUBROUTINE svdvar(v,ma,np,w,cvm,ncvm) INTEGER ma,ncvm,np,MMAXREAL cvm(ncvm,ncvm),v(np,np),w(np) PARAMETER (MMAX=20) Set to the maximum number of fit parameters. To evaluate the covariance matrix cvmof the fit for maparameters obtained by svdfit , call this routine with matrices v,was returned from svdfit .np,ncvm give the physical dimensions of v,w,cvm as indicated. INTEGER i,j,k REAL sum,wti(MMAX)do 11i=1,ma wti(i)=0. if(w(i).ne.0.) wti(i)=1./(w(i)*w(i)) enddo 11 do14i=1,ma Sum contributions to covariance matrix (15.4.20). do13j=1,i sum=0.do 12k=1,ma sum=sum+v(i,k)*v(j,k)*wti(k) enddo 12 cvm(i,j)=sumcvm(j,i)=sum enddo 13 enddo 14 returnEND 674 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).Examples Be aware that some apparently nonlinear problems can be expressed so that theyare linear. For example,an exponentialmodelwith two parameters aandb, y(x)=aexp(−bx)( 15.4.21 ) canberewrittenas log[y(x)] =c−bx (15.4.22 ) which is linear in its parameters candb. (Of course you must be aware that such transformations do not exactly take Gaussian errors into Gaussian errors.) Also watch out for “non-parameters, ”as in y(x)=aexp(−bx+d)( 15.4.23 ) Heretheparameters aanddare,infact,indistinguishable. Thisisagoodexampleof wherethenormalequationswillbeexactlysingular,andwhereSVDwill findazero singular value. SVD will then make a “least-squares ”choice for setting a balance between aandd(or, rather, their equivalents in the linear model derived by taking the logarithms). However —and this is true whenever SVD returns a zero singular value—you are better advised to figure out analytically where the degeneracy is amongyourbasis functions,and thenmake appropriatedeletionsin the basis set. Herearetwoexamplesforuser-suppliedroutines funcs. Thefirstoneistrivial andfits a general polynomial to a set of data: SUBROUTINE fpoly(x,p,np) INTEGER npREAL x,p(np) Fitting routine for a polynomial of degree np-1, with npcoefficients. INTEGER jp(1)=1.do 11j=2,np p(j)=p(j-1)*x enddo 11 return END Thesecondexampleisslightlylesstrivial. Itisusedto fitLegendrepolynomials up to some order nl-1through a data set. SUBROUTINE fleg(x,pl,nl) INTEGER nl REAL x,pl(nl) Fitting routine for an expansion with nlLegendre polynomials pl, evaluated using the recurrence relation as in §5.5. INTEGER j REAL d,f1,f2,twoxpl(1)=1.pl(2)=x if(nl.gt.2) then twox=2.*xf2=x d=1. 15.5NonlinearModels 675Sample 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).do11j=3,nl f1=d f2=f2+twox d=d+1.pl(j)=(f2*pl(j-1)-f1*pl(j-2))/d enddo 11 endif returnEND MultidimensionalFits Ifyouaremeasuringasinglevariable yasafunctionofmorethanonevariable —say,avectorofvariables x,thenyourbasisfunctionswillbefunctionsofavector, X1(x),...,X M(x). The χ2merit function is now χ2=N/summationdisplay i=1/bracketleftBigg yi−/summationtextM k=1akXk(xi) σi/bracketrightBigg2 (15.4.24 ) All of the preceding discussion goes through unchanged, with xreplaced by x.I n fact, if youare willing to tolerate a bit of programminghack, youcan use the above programs without any modi fication: In both lfitandsvdfit, the only use made ofthearrayelements x(i)isthateachelementisinturnpassedtotheuser-supplied routine funcs, which duly returns the values of the basis functions at that point. If you set x(i)=ibefore calling lfitorsvdfit, and independently provide funcs withthetruevectorvaluesofyourdatapoints(e.g.,ina COMMONblock),then funcs cantranslatefromthe fictitious x(i)’stotheactualdatapointsbeforedoingitswork. CITED REFERENCES AND FURTHER READING: Bevington, P.R. 1969, Data Reduction and Error Analysis for the Physical Sciences (New York: McGraw-Hill), Chapters 8–9. Lawson, C.L., and Hanson, R. 1974, Solving Least Squares Problems (Englewood Cliffs, NJ: Prentice-Hall). Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for Mathematical Computations (Englewood Cliffs, NJ: Prentice-Hall), Chapter 9. 15.5 Nonlinear Models We now consider fitting when the model depends nonlinearly on the set of M unknownparameters ak,k=1,2,...,M. We usethesameapproachas inprevious sections, namely to de fine aχ2merit function and determine best- fit parameters by its minimization. With nonlinear dependences, however, the minimization must proceed iteratively. Given trial values for the parameters, we develop a procedure that improves the trial solution. The procedure is then repeated until χ2stops (or effectively stops) decreasing.