Produced using MATLAB software simple nonlinear LS Page 1 of 10 simple nonlinear LS m This MATLAB function uses only standard MATLAB routines to implement the nonlinear fitting of parameters using the Marquardt method The errors for each response measurement are assumed to be independently distributed according to a normal distribution with a mean of zero and a common variance Once the parameter vector minimizing the residual sum of squarred errors is calculated approximate confidence intervalues are generated using the linearized design matrix around the minimum The normal equations are solved using singular value decomposition INPUT y the vector of observed responses calc yhat the name of a function that takes as input a given estimate of the model parameters and returns as output the vector of model predictions theta guess the initial guess of the vector of parameters alpha CI determines the width of the confidence intervals a value of 0 05 returns a 95 confidence interval Options a data structure that provides information to override the default options of the solver routine Fields are with default values in parentheses atol convergence tolerance of gradient norm 1e 8 max iter max of Gauss Newton iterations 100 max line max of weak line steps 10 tau LM the initial gues of the positive value added to the diagonal elements of the approximate Hessian in the Levenberg Marquardt method 1e 3 search angle max the maximum angle in degrees allowed between the search direction and the direction of steepest descent 80 verbose if 0 don t print to screen results If 1 print final results If 1 print results of each iteration 0 Param fix a data structure of fixed parameters that may be 7 16 2002 simple nonlinear LS Page 2 of 10 passed to the function calc yhat OUTPUT theta the parameter found as a local minimum of the residual sum of squared errors theta CI a vector containing the half widths for the confidence intervals for each parameter yhat the vector of model predictions using the fitted values of the model parameters yhat CI a vector containing the half widths for the confidence intervals for each model prediction Stats a data structure containing various measures of fit of the data RSS the residual sum of squared errors sample var the sample variance sample std the sample standard deviation R2 corr the R 2 correlation coefficient that reports the quality of the fit iflag an integer that takes a value of 1 upon successful completion of the routine K Beers MIT ChE 12 8 2001 function theta theta CI yhat yhat CI Stats iflag simple nonlinear LS y calc yhat theta guess alpha CI Options Param fix iflag 0 We extract the number of the measurements and the number of parameters to be fitted num exp length y num param length theta guess From the initial guess we now calculate the initial model predicions From this we then calculate the residual error vector theta theta guess yhat iflag2 feval calc yhat theta Param fix if iflag2 0 7 16 2002 simple nonlinear LS Page 3 of 10 iflag 1 error calc yhat returned error int2str iflag2 end residual y yhat RSS dot residual residual We now use finite differences to estimate the linearized design matrix X X iflag2 approx X FD theta yhat calc yhat Param fix if iflag2 0 iflag 2 error approx X FD returned error int2str iflag2 end We now set the parameters that direct the optimization algorithm if exist Options Options null 0 end atol 1e 8 if isfield Options atol atol Options atol end max iter 100 if isfield Options max iter max iter Options max iter end max line 10 if isfield Options max line max line Options max line end tau LM 1e 3 if isfield Options tau LM tau LM Options tau LM end search angle max 80 if isfield Options angle max search angle max Options search angle max end search cos min cos search angle max 180 pi verbose 0 if isfield Options verbose verbose Options verbose 7 16 2002 simple nonlinear LS Page 4 of 10 if verbose 0 verbose 0 end end if verbose disp disp Initial RSS num2str RSS end We now begin the iterations of the Gauss Newton method using the Levenberg Marquardt method to modify the approximate Hessian if the search direction is too far away from the direction of steepest descent if verbose disp disp Iteration RSS grad tau LM end coverged 0 for iter 1 max iter First we store the old value of the residual sum of squares to use in the descent criterion RSS old RSS We now calculate the search direction using the Levenberg Marquardt method of adding a diagonal matrix here we use Levenberg s choice of I to bias the search direction towards the steepest descent direction when the approximate Hessian behaves poorly Since we are adding a postive diagonal matrix problems with near singularity are somewhat reduced We use a QR decomposition to solve We use a while loop to keep calculating search directions until we find one that is less than the maximum angle away from the direction of steepest descent accept search 0 d X residual while accept search Q R qr X X tau LM eye num param p R Q X residual try search cos dot p d sqrt dot p p 7 16 2002 simple nonlinear LS Page 5 of 10 sqrt dot d d if search cos search cos min accept search 1 tau LM 0 1 tau LM else tau LM 10 tau LM end catch accept search 1 end end Now that we have selected the search direction we begin a weak line search to force acceptance of the descent criterion ifound descent 0 for iter line 0 max line We try the following fractional step length alpha 2 iter line We now calculate the model predictions for this new estimate of the parameter vector yhat feval calc yhat theta alpha p Param fix RSS dot y yhat y yhat If this fractional step satisfies the descent criterion accept it Otherwise try a smaller fractional step if RSS RSS old theta theta alpha p X approx X FD theta yhat calc yhat Param fix residual y yhat ifound descent 1 break end end If we have not been able to decrease the cost function increase tau LM and try again if ifound descent tau LM 10 tau LM end We print the results of the iteration to the screen if requested if verbose 1 disp int2str iter num2str RSS 7 16 2002 simple nonlinear LS num2str sqrt dot d d num2str tau LM Page 6 of 10 end If we have updated the estimate we check the magnitude of the gradient for convergence d the direction of steepest descent is merely the negative of the gradient and so we use this vector d X residual if dot d d atol atol converged 1 break end end iflag converged if verbose disp disp Final parameter vector estimate disp theta disp converged