non_lin_fit

NLINFIT Nonlinear least-squares data fitting by the Gauss-Newton method.
   NLINFIT(X,Y,FUN,BETA0) estimates the coefficients of a nonlinear
   function using least squares.  Y is a vector of response (dependent
   variable) values.  Typically, X is a design matrix of predictor
   (independent variable) values, with one row for each value in Y.
   However, X may be any array that FUN is prepared to accept.  FUN is
   a function that accepts two arguments, a coefficient vector and the
   array X, and returns a vector of fitted Y values.  BETA0 is a vector
   containing initial values for the coefficients.

   [BETA,R,J] = NLINFIT(X,Y,FUN,BETA0) returns the fitted coefficients
   BETA, the residuals R, and the Jacobian J.  You can use these outputs
   with NLPREDCI to produce error estimates on predictions, and with
   NLPARCI to produce error estimates on the estimated coefficients.

0001 function [beta,r,J] = non_lin_fit(X,y,model,beta0)
0002 %NLINFIT Nonlinear least-squares data fitting by the Gauss-Newton method.
0003 %   NLINFIT(X,Y,FUN,BETA0) estimates the coefficients of a nonlinear
0004 %   function using least squares.  Y is a vector of response (dependent
0005 %   variable) values.  Typically, X is a design matrix of predictor
0006 %   (independent variable) values, with one row for each value in Y.
0007 %   However, X may be any array that FUN is prepared to accept.  FUN is
0008 %   a function that accepts two arguments, a coefficient vector and the
0009 %   array X, and returns a vector of fitted Y values.  BETA0 is a vector
0010 %   containing initial values for the coefficients.
0011 %
0012 %   [BETA,R,J] = NLINFIT(X,Y,FUN,BETA0) returns the fitted coefficients
0013 %   BETA, the residuals R, and the Jacobian J.  You can use these outputs
0014 %   with NLPREDCI to produce error estimates on predictions, and with
0015 %   NLPARCI to produce error estimates on the estimated coefficients.
0016 %
0017 if (nargin<4), error('non_lin_fit requires four arguments.'); end
0018 
0019 if min(size(y)) ~= 1
0020    error('Requires a vector second input argument.');
0021 end
0022 y = y(:);
0023 
0024 if size(X,1) == 1 % turn a row vector into a column vector.
0025    X = X(:);
0026 end
0027 
0028 wasnan = (isnan(y) | any(isnan(X),2));
0029 if (any(wasnan))
0030    y(wasnan) = [];
0031    X(wasnan,:) = [];
0032 end
0033 n = length(y);
0034 
0035 p = length(beta0);
0036 beta0 = beta0(:);
0037 
0038 if any(size(feval(model,beta0,X)) ~= size(y))
0039    error('FUN should return a column vector of the same length as Y.');
0040 end
0041 
0042 J = zeros(n,p);
0043 beta = beta0;
0044 betanew = beta + 1;
0045 maxiter = 100;
0046 iter = 0;
0047 betatol = 1.0E-4;
0048 rtol = 1.0E-4;
0049 sse = 1;
0050 sseold = sse;
0051 seps = sqrt(eps);
0052 zbeta = zeros(size(beta));
0053 s10 = sqrt(10);
0054 eyep = eye(p);
0055 zerosp = zeros(p,1);
0056 
0057 while (norm((betanew-beta)./(beta+seps)) > betatol | abs(sseold-sse)/(sse+seps) > rtol) & iter < maxiter
0058    if iter > 0, 
0059       beta = betanew;
0060    end
0061 
0062    iter = iter + 1;
0063    yfit = feval(model,beta,X);
0064    r = y - yfit;
0065    sseold = r'*r;
0066 
0067    for k = 1:p,
0068       delta = zbeta;
0069       if (beta(k) == 0)
0070          nb = sqrt(norm(beta));
0071          delta(k) = seps * (nb + (nb==0));
0072       else
0073          delta(k) = seps*beta(k);
0074       end
0075       yplus = feval(model,beta+delta,X);
0076       J(:,k) = (yplus - yfit)/delta(k);
0077    end
0078 
0079    Jplus = [J;(1.0E-2)*eyep];
0080    rplus = [r;zerosp];
0081 
0082    % Levenberg-Marquardt type adjustment
0083    % Gauss-Newton step -> J\r
0084    % LM step -> inv(J'*J+constant*eye(p))*J'*r
0085    step = Jplus\rplus;
0086    
0087    betanew = beta + step;
0088    yfitnew = feval(model,betanew,X);
0089    rnew = y - yfitnew;
0090    sse = rnew'*rnew;
0091    iter1 = 0;
0092    while sse > sseold & iter1 < 12
0093       step = step/s10;
0094       betanew = beta + step;
0095       yfitnew = feval(model,betanew,X);
0096       rnew = y - yfitnew;
0097       sse = rnew'*rnew;
0098       iter1 = iter1 + 1;
0099    end
0100 end
0101 if iter == maxiter
0102    disp(['non_lin_fit did NOT converge in ', int2str(iter),' iterations. Returning results from last iteration.']);
0103    else
0104    disp(['non_lin_fit converges in ', int2str(iter),' iterations.']);
0105 end