f15-5

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 modification: 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 cantranslatefromthefictitious 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 define 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. 676 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).Howisthisproblemdifferentfromthegeneralnonlinearfunctionminimization problem already dealt with in Chapter 10? Superficially, not at all: Sufficientlyclose to the minimum, we expect the χ 2function to be well approximated by a quadratic form, which we can write as χ2(a)≈γ−d·a+1 2a·D·a (15.5.1 ) wheredis an M-vector and Dis an M×Mmatrix. (Compare equation 10.6.1.) If the approximation is a good one, we know how to jump from the current trial parameters acurto the minimizing ones aminin a single leap, namely amin =acur +D−1·/bracketleftbig −∇χ2(acur)/bracketrightbig (15.5.2 ) (Compare equation 10.7.4.) On the other hand, (15.5.1) might be a poor local approximation to the shape of the function that we are trying to minimize at acur. In that case, about all we can do is take a step down the gradient, as in the steepest descent method ( §10.6). In other words, anext =acur−constant ×∇χ2(acur)( 15.5.3 ) where the constant is small enough not to exhaust the downhill direction. To use (15.5.2) or (15.5.3), we must be able to compute the gradient of the χ2 functionatanysetofparameters a. Touse(15.5.2)wealsoneedthematrix D,which is thesecondderivativematrix(Hessian matrix)of the χ2meritfunction,at any a. Now, this is the crucial difference from Chapter 10: There, we had no way of directly evaluating the Hessian matrix. We were given only the ability to evaluate the function to be minimized and (in some cases) its gradient. Therefore, we had to resort to iterative methods not justbecause our function was nonlinear, but also in order to build up information about the Hessian matrix. Sections 10.7 and 10.6concernedthemselveswithtwodifferenttechniquesforbuildingupthisinformation. Here, life is much simpler. We knowexactly the form of χ 2, since it is based on a modelfunctionthat we ourselveshavespecified. Thereforethe Hessian matrixis known to us. Thus we are free to use (15.5.2) whenever we care to do so. The only reason to use (15.5.3) will be failure of (15.5.2) to improve the fit, signaling failure of (15.5.1) as a good local approximation. Calculation ofthe Gradientand Hessian The model to be fitted is y=y(x;a)( 15.5.4 ) and the χ2merit function is χ2(a)=N/summationdisplay i=1/bracketleftbiggyi−y(xi;a) σi/bracketrightbigg2 (15.5.5 ) 15.5NonlinearModels 677Sample 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).The gradient of χ2with respect to the parameters a, which will be zero at the χ2 minimum, has components ∂χ2 ∂a k=−2N/summationdisplay i=1[yi−y(xi;a)] σ2 i∂y(xi;a) ∂a kk=1,2,...,M (15.5.6 ) Taking an additional partial derivative gives ∂2χ2 ∂a k∂a l=2N/summationdisplay i=11 σ2 i/bracketleftbigg∂y(xi;a) ∂a k∂y(xi;a) ∂a l−[yi−y(xi;a)]∂2y(xi;a) ∂a l∂a k/bracketrightbigg (15.5.7 ) It is conventional to remove the factors of 2 by defining βk≡−1 2∂χ2 ∂a kαkl≡1 2∂2χ2 ∂a k∂a l(15.5.8 ) making [α]=1 2Din equation (15.5.2), in terms of which that equation can be rewritten as the set of linear equations M/summationdisplay l=1αklδa l=βk (15.5.9 ) This set is solved for the increments δa lthat, added to the current approximation, givethe nextapproximation. Inthe contextof least-squares,thematrix [α], equalto one-half times the Hessian matrix, is usually called the curvature matrix . Equation (15.5.3), the steepest descent formula, translates to δa l=constant ×βl (15.5.10 ) Notethatthecomponents αkloftheHessianmatrix(15.5.7)dependbothonthe first derivatives and on the second derivatives of the basis functions with respect to their parameters. Some treatments proceed to ignore the second derivative without comment. We will ignore it also, but only aftera few comments. Secondderivativesoccurbecausethegradient(15.5.6)alreadyhasadependence on∂y/∂a k,sothenextderivativesimplymustcontaintermsinvolving ∂2y/∂a l∂a k. The second derivative term can be dismissed when it is zero (as in the linear case of equation 15.4.8), or small enough to be negligible when compared to the term involvingthe first derivative. It also has an additionalpossibility of beingignorablysmall in practice: The term multiplying the second derivative in equation (15.5.7) is[y i−y(xi;a)]. For a successful model, this term should just be the random measurement error of each point. This error can have either sign, and should ingeneralbe uncorrelatedwith themodel. Therefore,the secondderivativeterms tend to cancel out when summed over i. Inclusionofthesecond-derivativetermcaninfactbedestabilizingifthemodel fits badly or is contaminated by outlier points that are unlikely to be offset by 678 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).compensating points of opposite sign. From this point on, we will always use as the definition of αklthe formula αkl=N/summationdisplay i=11 σ2 i/bracketleftbigg∂y(xi;a) ∂a k∂y(xi;a) ∂a l/bracketrightbigg (15.5.11 ) This expression more closely resembles its linear cousin (15.4.8). You should understand that minor (or even major) fiddling with [α]has no effect at all on what final set of parameters ais reached, but affects only the iterative route that is taken in getting there. The condition at the χ2minimum, that βk=0for all k, is independent of how [α]is defined. Levenberg-Marquardt Method Marquardt [1]hasputforthanelegantmethod,relatedtoanearliersuggestionof Levenberg,forvaryingsmoothlybetweentheextremesoftheinverse-Hessianmethod (15.5.9)andthesteepestdescentmethod(15.5.10). Thelattermethodisusedfarfrom the minimum, switching continuouslyto the formeras the minimumis approached.ThisLevenberg-Marquardtmethod (alsocalled Marquardtmethod )worksverywell in practice and has become the standard of nonlinearleast-squares routines. The method is based on two elementary, but important, insights. Consider the “constant” in equation (15.5.10). What should it be, even in order of magnitude? What sets its scale? There is no informationabout the answer in the gradient. That tells only the slope, not how far that slope extends. Marquardt’s first insight is thatthe components of the Hessian matrix, even if they are not usable in any precise fashion,give someinformationabouttheorder-of-magnitudescale of theproblem. The quantity χ 2is nondimensional,i.e., is a pure number; this is evident from its definition (15.5.5). On the other hand, βkhas the dimensions of 1/a k, which may well be dimensional,i.e., have units like cm−1, or kilowatt-hours,or whatever. (In fact, each component of βkcan have different dimensions!) The constant of proportionalitybetween βkandδa kmustthereforehavethedimensionsof a2 k. Scan thecomponentsof [α]andyouseethatthereisonlyoneobviousquantitywiththese dimensions, and that is 1/α kk, the reciprocalof the diagonal element. So that must set the scale of the constant. But that scale might itself be too big. So let’s divide theconstantbysome(nondimensional)fudgefactor λ,withthepossibilityofsetting λ/greatermuch1to cut down the step. In other words, replace equation (15.5.10)by δa l=1 λα llβlor λα llδa l=βl (15.5.12 ) It is necessary that αllbe positive, but this is guaranteed by definition (15.5.11) — another reason for adopting that equation. Marquardt’s second insight is that equations (15.5.12) and (15.5.9) can be combined if we define a new matrix α/primeby the following prescription α/prime jj≡αjj(1 +λ) α/prime jk≡αjk (j/negationslash=k)(15.5.13 ) 15.5NonlinearModels 679Sample 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 then replace both (15.5.12) and (15.5.9) by M/summationdisplay l=1α/prime klδa l=βk (15.5.14 ) When λis very large, the matrix α/primeis forced into being diagonally dominant ,s o equation (15.5.14) goes over to be identical to (15.5.12). On the other hand, as λ approaches zero, equation (15.5.14) goes over to (15.5.9). Given an initial guess for the set of fitted parameters a, the recommended Marquardt recipe is as follows: •Compute χ2(a). •Pick a modest value for λ, say λ=0.001. •(†) Solvethe linear equations(15.5.14)for δaandevaluate χ2(a+δa). •Ifχ2(a+δa)≥χ2(a),increase λby a factor of 10 (or any other substantial factor) and go back to ( †). •Ifχ2(a+δa)<χ2(a),decrease λby a factor of 10, update the trial solutiona←a+δa, and go back to ( †). Alsonecessaryisaconditionforstopping. Iteratingtoconvergence(tomachine accuracy or to the roundoff limit) is generally wasteful and unnecessary since the minimum is at best only a statistical estimate of the parameters a. As we will see in§15.6, a change in the parameters that changes χ2by an amount /lessmuch 1isnever statistically meaningful. Furthermore, it is not uncommon to find the parameters wandering around near the minimum in a flat valley of complicated topography. The rea- son is that Marquardt’smethodgeneralizesthemethodofnormalequations( §15.4), hence has the same problem as that method with regard to near-degeneracy of the minimum. Outright failure by a zero pivot is possible, but unlikely. More often, a small pivot will generate a large correction which is then rejected, the value ofλbeing then increased. For sufficiently large λthe matrix [α /prime]is positive definite and can have no small pivots. Thus the method does tend to stay away from zero pivots, but at the cost of a tendency to wander around doing steepest descent in very un-steep degenerate valleys. These considerations suggest that, in practice, one might as well stop iterating on the first or second occasion that χ2decreases by a negligible amount, say either less than 0.01absolutely or (in case roundoff prevents that being reached) some fractional amount like 10−3. Don’t stop after a step where χ2increases: That only shows that λhas not yet adjusted itself optimally. Once the acceptable minimum has been found, one wants to set λ=0and compute the matrix [C]≡[α]−1(15.5.15 ) which, as before, is the estimated covariance matrix of the standard errors in the fitted parameters a(see next section). The following pair of subroutines encodes Marquardt’s method for nonlinear parameter estimation. Much of the organization matches that used in lfitof §15.4. Inparticularthe array ia(1:ma) must be inputwithcomponentsoneor zero corresponding to whether the respective parameter values a(1:ma) are to be fitted for or held fixed at their input values, respectively. 680 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).The routine mrqminperforms one iteration of Marquardt’s method. It is first called (once) with alamda <0, which signals the routine to initialize. alamdais returned on the first and all subsequent calls as the suggested value of λfor the next iteration; aandchisqare always returned as the best parameters found so far and their χ2. When convergence is deemed satisfactory, set alamdato zero before a final call. The matrices alphaandcovar(which were used as workspace in all previous calls) will then be set to the curvature and covariance matrices for the converged parameter values. The arguments alpha,a, and chisqmust not be modified between calls, nor should alamdabe, except to set it to zero for the final call. When an uphill step is taken, chisqandaare returned with their input (best) values, but alamdais returned with an increased value. Theroutine mrqmincallstheroutine mrqcofforthecomputationofthematrix [α](equation 15.5.11) and vector β(equations 15.5.6 and 15.5.8). In turn mrqcof calls theuser-suppliedroutine funcs(x,a,y,dyda) ,whichforinputvalues x≡xi anda≡areturns the model function y≡y(xi;a)and the vector of derivatives dyda≡∂y/∂a k. SUBROUTINE mrqmin(x,y,sig,ndata,a,ia,ma,covar,alpha,nca, * chisq,funcs,alamda) INTEGER ma,nca,ndata,ia(ma),MMAXREAL alamda,chisq,funcs,a(ma),alpha(nca,nca),covar(nca,nca), * sig(ndata),x(ndata),y(ndata) PARAMETER (MMAX=20) Set to largest number of fit parameters. C USES covsrt,gaussj,mrqcof Levenberg-Marquardt method, attempting to reduce the value χ2of a fit between a set of data points x(1:ndata) ,y(1:ndata) with individual standard deviations sig(1:ndata) , and a nonlinear function dependent on macoefficients a(1:ma) . The input array ia(1:ma) indicates 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 current best-fit values for the parameters a(1:ma) ,a n d χ2=chisq.T h ea r - rayscovar(1:nca,1:nca) ,alpha(1:nca,1:nca) with physical dimension nca(≥the number of fitted parameters) are used as working space during most iterations. Supply a subroutine funcs(x,a,yfit,dyda,ma) that evaluates the fitting function yfit,a n di t s derivatives dydawith respect to the fitting parameters aatx. On the first call provide an initial guess for the parameters a,a n ds e t alamda<0 for initialization (which then sets alamda=.001 ). If a step succeeds chisqbecomes smaller and alamda decreases by a factor of 10. If a step fails alamda grows by a factor of 10. You must call this routine repeatedly until convergence is achieved. Then, make one final call with alamda=0 ,s o thatcovar(1:ma,1:ma) returns the covariance matrix, and alphathe curvature matrix. (Parameters held fixed will return zero covariances.) INTEGER j,k,l,mfitREAL ochisq,atry(MMAX),beta(MMAX),da(MMAX) SAVE ochisq,atry,beta,da,mfit if(alamda.lt.0.)then Initialization. mfit=0do 11j=1,ma if (ia(j).ne.0) mfit=mfit+1 enddo 11 alamda=0.001 call mrqcof(x,y,sig,ndata,a,ia,ma,alpha,beta,nca,chisq,funcs) ochisq=chisqdo 12j=1,ma atry(j)=a(j) enddo 12 endif do14j=1,mfit Alter linearized fitting matrix, by augmenting diagonal elements. do13k=1,mfit 15.5NonlinearModels 681Sample 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).covar(j,k)=alpha(j,k) enddo 13 covar(j,j)=alpha(j,j)*(1.+alamda)da(j)=beta(j) enddo 14 call gaussj(covar,mfit,nca,da,1,1) Matrix solution. if(alamda.eq.0.)then Once converged, evaluate covariance matrix. call covsrt(covar,nca,ma,ia,mfit)call covsrt(alpha,nca,ma,ia,mfit) Spread out alphato its full size too. return endifj=0do 15l=1,ma Did the trial succeed? if(ia(l).ne.0) then j=j+1atry(l)=a(l)+da(j) endif enddo 15 call mrqcof(x,y,sig,ndata,atry,ia,ma,covar,da,nca,chisq,funcs)if(chisq.lt.ochisq)then Success, accept the new solution. alamda=0.1*alamda ochisq=chisqdo 17j=1,mfit do16k=1,mfit alpha(j,k)=covar(j,k) enddo 16 beta(j)=da(j) enddo 17 do18l=1,ma a(l)=atry(l) enddo 18 else Failure, increase alamda and return. alamda=10.*alamdachisq=ochisq endif returnEND Notice the use of the routine covsrtfrom§15.4. This is merely for rearranging the covariancematrix covarinto the order of all maparameters. The aboveroutine also makes use of SUBROUTINE mrqcof(x,y,sig,ndata,a,ia,ma,alpha,beta,nalp, * chisq,funcs) INTEGER ma,nalp,ndata,ia(ma),MMAX REAL chisq,a(ma),alpha(nalp,nalp),beta(ma),sig(ndata),x(ndata), * y(ndata) EXTERNAL funcs PARAMETER (MMAX=20) Used by mrqmin to evaluate the linearized fitting matrix alpha, and vector betaas in (15.5.8), and calculate χ2. INTEGER mfit,i,j,k,l,m REAL dy,sig2i,wt,ymod,dyda(MMAX)mfit=0do 11j=1,ma if (ia(j).ne.0) mfit=mfit+1 enddo 11 do13j=1,mfit Initialize (symmetric) alpha,beta. do12k=1,j 682 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).alpha(j,k)=0. enddo 12 beta(j)=0. enddo 13 chisq=0. do16i=1,ndata Summation loop over all data. call funcs(x(i),a,ymod,dyda,ma)sig2i=1./(sig(i)*sig(i))dy=y(i)-ymod j=0 do 15l=1,ma if(ia(l).ne.0) then j=j+1 wt=dyda(l)*sig2i k=0do 14m=1,l if(ia(m).ne.0) then k=k+1alpha(j,k)=alpha(j,k)+wt*dyda(m) endif enddo 14 beta(j)=beta(j)+dy*wt endif enddo 15 chisq=chisq+dy*dy*sig2i And find χ2. enddo 16 do18j=2,mfit Fill in the symmetric side. do17k=1,j-1 alpha(k,j)=alpha(j,k) enddo 17 enddo 18 returnEND Example The following subroutine fgaussis an example of a user-supplied subroutine funcs. Used with the above routine mrqmin(in turn using mrqcof,covsrt, and gaussj), it fits for the model y(x)=K/summationdisplay k=1Bkexp/bracketleftBigg −/parenleftbiggx−Ek Gk/parenrightbigg2/bracketrightBigg (15.5.16 ) which is a sum of KGaussians, each having a variable position, amplitude, and width. We store the parameters in the order B1,E 1,G 1,B 2,E 2,G 2,...,B K, EK,G K. 15.5NonlinearModels 683Sample 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).SUBROUTINE fgauss(x,a,y,dyda,na) INTEGER na REAL x,y,a(na),dyda(na) y(x;a)is the sum of na/3Gaussians (15.5.16). The amplitude, center, and width of the Gaussians are stored in consecutive locations of a:a(i) =Bk,a(i+1) =Ek,a(i+2) = Gk,k=1, ...,na/3. INTEGER iREAL arg,ex,facy=0. do 11i=1,na-1,3 arg=(x-a(i+1))/a(i+2)ex=exp(-arg**2)fac=a(i)*ex*2.*arg y=y+a(i)*ex dyda(i)=exdyda(i+1)=fac/a(i+2) dyda(i+2)=fac*arg/a(i+2) enddo 11 returnEND More Advanced Methods for NonlinearLeast Squares The Levenberg-Marquardt algorithm can be implemented as a model-trust region method for minimization (see §9.7 and ref. [2]) applied to the special case of a least squares function. A code of this kind due to Mor ´e[3]can be found in MINPACK [4]. Another algorithm for nonlinear least-squares keeps the second- derivativeterm we droppedin the Levenberg-Marquardtmethod wheneverit would be better to do so. These methods are called “full Newton-type” methods and are reputed to be more robust than Levenberg-Marquardt, but more complex. One implementation is the code NL2SOL [5]. CITED REFERENCES AND FURTHER READING: Bevington, P.R. 1969, Data Reduction and Error Analysis for the Physical Sciences (New York: McGraw-Hill), Chapter 11. Marquardt, D.W. 1963, Journal of the Society for Industrial and Applied Mathematics , vol. 11, pp. 431–441. [1] Jacobs,D.A.H.(ed.)1977, TheStateoftheArtinNumericalAnalysis (London:AcademicPress), Chapter III.2 (by J.E. Dennis). Dennis,J.E., andSchnabel,R.B. 1983, NumericalMethods forUnconstrained Optimizationand Nonlinear Equations (Englewood Cliffs, NJ: Prentice-Hall). [2] Mor´e, J.J. 1977, in Numerical Analysis , Lecture Notes in Mathematics, vol. 630, G.A. Watson, ed. (Berlin: Springer-Verlag), pp. 105–116. [3] Mor´e,J.J.,Garbow,B.S.,andHillstrom,K.E.1980, UserGuideforMINPACK-1 ,ArgonneNational Laboratory Report ANL-80-74. [4] Dennis, J.E., Gay, D.M, and Welsch, R.E. 1981, ACM Transactions on Mathematical Software , vol. 7, pp. 348–368; op. cit., pp. 369–383. [5]. 684 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).15.6 Confidence Limits on Estimated Model Parameters Severaltimesalreadyinthischapterwehavemadestatementsaboutthestandard errors, or uncertainties, in a set of Mestimated parameters a. We have given some formulas for computing standard deviations or variances of individual parameters (equations 15.2.9, 15.4.15, 15.4.19), as well as some formulas for covariances between pairs of parameters (equation 15.2.10; remark following equation 15.4.15;equation 15.4.20; equation 15.5.15). In this section, we want to be more explicit regarding the precise meaning of these quantitative uncertainties, and to give further information about how quantitative confidence limits on fitted parameters can be estimated. The subject can get somewhat technical, and even somewhat confusing, so we will try to makeprecise statements, even when they must be offered without proof. Figure 15.6.1 shows the conceptual scheme of an experiment that “measures” a set of parameters. There is some underlying true set of parameters a truethat are known to Mother Nature but hidden from the experimenter. These true parameters arestatistically realized,alongwithrandommeasurementerrors,asameasureddata set,whichwewillsymbolizeas D(0). Thedataset D(0)isknowntotheexperimenter. He or she fits the data to a model by χ2minimizationor some other technique, and obtainsmeasured,i.e.,fitted, valuesfortheparameters,whichweheredenote a(0). Because measurement errors have a random component, D(0)is not a unique realization of the true parameters atrue. Rather, there are infinitely many other realizations of the true parameters as “hypothetical data sets” each of which could have been the one measured, but happened not to be. Let us symbolize these byD(1),D(2),.... Each one, had it been realized, would have given a slightly different set of fitted parameters, a(1),a(2),..., respectively. These parameter sets a(i)therefore occur with some probability distribution in the M-dimensional space of all possible parametersets a. The actual measured set a(0)is one memberdrawn from this distribution. Even more interesting than the probability distribution of a(i)would be the distribution of the difference a(i)−atrue. This distribution differs from the former onebyatranslationthatputsMotherNature’struevalueattheorigin. Ifweknew this distribution, we would know everything that there is to know about the quantitative uncertainties in our experimental measurement a(0). So the name of the game is to find some way of estimating or approximating theprobabilitydistributionof a(i)−atruewithoutknowing atrueandwithouthaving available to us an infinite universe of hypothetical data sets. Monte CarloSimulationof Synthetic Data Sets Although the measured parameter set a(0)is not the true one, let us consider a fictitious world in which it wasthe true one. Since we hope that our measured parameters are not toowrong, we hope that that fictitious world is not too different from the actual world with parameters atrue. In particular, let us hope — no, let us assume— that the shape of the probability distribution a(i)−a(0)in the fictitious worldisthesame,orverynearlythesame,astheshapeoftheprobabilitydistribution