help-octave
[Top][All Lists]

## Re: nonlinear least square fitting

 From: Francesco Potorti` Subject: Re: nonlinear least square fitting Date: Mon, 8 Feb 1999 22:37:33 +0100 (CET)

```This is becoming a FAQ.  I'm  posting the latest version of leasqr.m that I
found around, tweaked by me for Octave.  I got the original source from the
matlab site.  I'd like to ask the  authors if they wish to publish it under
the GNU  public license, or  leave it  in the public  domain, but I  am not
acquainted in these  issues, so I'm not really sure about  what to ask.  If
anyone wants  to do that, Octave may  at last have a  nolinear least square
fitting function.

--
Francesco Potortì (researcher)         Voice:    +39-050-593 203 (op. 211)
Computer Networks Group                Fax:      +39-050-904052
CNUCE-CNR, Via Santa Maria 36          Email:    address@hidden
56126 Pisa - Italy                     Web:      http://fly.cnuce.cnr.it/

=========================== dfdp.m =============================
function prt=dfdp(x,f,p,dp,func)
% numerical partial derivatives (Jacobian) df/dp for use with leasqr
% --------INPUT VARIABLES---------
% x=vec or matrix of indep var(used as arg to func) x=[x0 x1 ....]
% f=func(x,p) vector initialsed by user before each call to dfdp
% p= vec of current parameter values
% dp= fractional increment of p for numerical derivatives
%      dp(j)>0 central differences calculated
%      dp(j)<0 one sided differences calculated
%      dp(j)=0 sets corresponding partials to zero; i.e. holds p(j) fixed
% func=string naming the function (.m) file
%      e.g. to calc Jacobian for function expsum prt=dfdp(x,f,p,dp,'expsum')
%----------OUTPUT VARIABLES-------
% prt= Jacobian Matrix prt(i,j)=df(i)/dp(j)
%================================
m=length(x);n=length(p);      %dimensions
ps=p; prt=zeros(m,n);del=zeros(n,1);       % initialise Jacobian to Zero
for j=1:n
del(j)=dp(j) .*p(j);    %cal delx=fract(dp)*param value(p)
if p(j)==0
del(j)=dp(j);     %if param=0 delx=fraction
end
p(j)=ps(j) + del(j);
if del(j)~=0, f1=feval(func,x,p);
if dp(j) < 0, prt(:,j)=(f1-f)./del(j);
else
p(j)=ps(j)- del(j);
prt(:,j)=(f1-feval(func,x,p))./(2 .*del(j));
end
end
p(j)=ps(j);     %restore p(j)
end
return

======================== leasqr.m ==============================
function [f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2]= ...
leasqr(x,y,pin,F,stol,niter,wt,dp,dFdp,options)
%function[f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2]=
%                   leasqr(x,y,pin,F,{stol,niter,wt,dp,dFdp,options})
%
% Version 3.beta
%  {}= optional parameters
% Levenberg-Marquardt nonlinear regression of f(x,p) to y(x), where:
% x=vec or mat of indep variables, 1 row/observation: x=[x0 x1....xm]
% y=vec of obs values, same no. of rows as x.
% wt=vec(dim=length(x)) of statistical weights.  These should be set
%   to be proportional to (sqrt of var(y))^-1; (That is, the covariance
%   matrix of the data is assumed to be proportional to diagonal with diagonal
%   equal to (wt.^2)^-1.  The constant of proportionality will be estimated.),
%   default=ones(length(y),1).
% pin=vector of initial parameters to be adjusted by leasqr.
% dp=fractional incr of p for numerical partials,default= .001*ones(size(pin))
%   dp(j)>0 means central differences.
%   dp(j)<0 means one-sided differences.
% Note: dp(j)=0 holds p(j) fixed i.e. leasqr wont change initial guess: pin(j)
% F=name of function in quotes,of the form y=f(x,p)
% dFdp=name of partials M-file in quotes default is prt=dfdp(x,f,p,dp,F)
% stol=scalar tolerances on fractional improvement in ss,default stol=.0001
% niter=scalar max no. of iterations, default = 20
% options=matrix of n rows (same number of rows as pin) containing
%   column 1: desired fractional precision in parameter estimates.
%     Iterations are terminated if change in parameter vector (chg) on two
%     consecutive iterations is less than their corresponding elements
%     in options(:,1).  [ie. all(abs(chg*current parm est) < options(:,1))
%      on two consecutive iterations.], default = zeros().
%   column 2: maximum fractional step change in parameter vector.
%     Fractional change in elements of parameter vector is constrained to be
%     at most options(:,2) between sucessive iterations.
%     [ie. abs(chg(i))=abs(min([chg(i) options(i,2)*current param estimate])).],
%     default = Inf*ones().
%
%          OUTPUT VARIABLES
% f=vec function values computed in function func.
% p=vec trial or final parameters. i.e, the solution.
% kvg=scalar: =1 if convergence, =0 otherwise.
% iter=scalar no. of interations used.
% corp= correlation matrix for parameters
% covp= covariance matrix of the parameters
% covr = diag(covariance matrix of the residuals)
% stdresid= standardized residuals
% Z= matrix that defines confidence region
% r2= coefficient of multiple determination

% All Zero guesses not acceptable
% Richard I. Shrager (301)-496-1122
% Modified by A.Jutan (519)-679-2111
% Modified by Ray Muzic 14-Jul-1992
%       1) add maxstep feature for limiting changes in parameter estimates
%          at each step.
%       2) remove forced columnization of x (x=x(:)) at beginning. x could be
%          a matrix with the ith row of containing values of the
%          independent variables at the ith observation.
%       4) add optional return arguments covp, stdresid, chi2
%       5) revise estimates of corp, stdev
% Modified by Ray Muzic 11-Oct-1992
%       1) revise estimate of Vy.  remove chi2, add Z as return values
% Modified by Ray Muzic 7-Jan-1994
%       1) Replace ones(x) with a construct that is compatible with versions
%          newer and older than v 4.1.
%       2) Added global declaration of verbose (needed for newer than v4.x)
%       3) Replace return value var, the variance of the residuals with covr,
%          the covariance matrix of the residuals.
%       4) Introduce options as 10th input argument.  Include
%          convergence criteria and maxstep in it.
%       5) Correct calculation of xtx which affects coveraince estimate.
%       6) Eliminate stdev (estimate of standard deviation of parameter
%          estimates) from the return values.  The covp is a much more
%          meaningful expression of precision because it specifies a confidence
%          region in contrast to a confidence interval..  If needed, however,
%          stdev may be calculated as stdev=sqrt(diag(covp)).
%       7) Change the order of the return values to a more logical order.
%       8) Change to more efficent algorithm of Bard for selecting epsL.
%       9) Tighten up memory usage by making use of sparse matrices (if
%          MATLAB version >= 4.0) in computation of covp, corp, stdresid.
% Modified by Sean Brennan 17-May-1994
%          verbose is now a vector:
%          verbose(1) controls output of results
%          verbose(2) controls plotting intermediate results
%
% References:
% Bard, Nonlinear Parameter Estimation, Academic Press, 1974.
% Draper and Smith, Applied Regression Analysis, John Wiley and Sons, 1981.
%
%set default args

% argument processing
%

plotcmd='plot(x(:,1),y,''+'',x(:,1),f); shg';
%if (sscanf(version,'%f') >= 4),
vernum= sscanf(version,'%f');
if vernum(1) >= 4,
global verbose
plotcmd='plot(x(:,1),y,''+'',x(:,1),f); figure(gcf)';
end;
if (exist('OCTAVE_VERSION'))
global verbose
end;

if(exist('verbose')~=1), %If verbose undefined, print nothing
verbose(1)=0;    %This will not tell them the results
verbose(2)=0;    %This will not replot each loop
end;

if (nargin <= 8), dFdp='dfdp'; end;
if (nargin <= 7), dp=.001*(pin*0+1); end; %DT
if (nargin <= 6), wt=ones(length(y),1); end;    % SMB modification
if (nargin <= 5), niter=20; end;
if (nargin == 4), stol=.0001; end;
%

y=y(:); wt=wt(:); pin=pin(:); dp=dp(:); %change all vectors to columns
% check data vectors- same length?
m=length(y); n=length(pin); p=pin;[m1,m2]=size(x);
if m1~=m ,error('input(x)/output(y) data must have same number of rows ') ,end;

if (nargin <= 9),
options=[zeros(n,1) Inf*ones(n,1)];
nor = n; noc = 2;
else
[nor noc]=size(options);
if (nor ~= n),
error('options and parameter matrices must have same number of rows'),
end;
if (noc ~= 2),
options=[options(noc,1) Inf*ones(noc,1)];
end;
end;
pprec=options(:,1);
maxstep=options(:,2);
%

% set up for iterations
%
f=feval(F,x,p); fbest=f; pbest=p;
r=wt.*(y-f);
sbest=r'*r;
nrm=zeros(n,1);
chgprev=Inf*ones(n,1);
kvg=0;
epsLlast=1;
epstab=[.1 1 1e2 1e4 1e6];

% do iterations
%
for iter=1:niter,
pprev=pbest;
prt=feval(dFdp,x,fbest,pprev,dp,F);
r=wt.*(y-fbest);
sprev=sbest;
sgoal=(1-stol)*sprev;
for j=1:n,
if dp(j)==0,
nrm(j)=0;
else
prt(:,j)=wt.*prt(:,j);
nrm(j)=prt(:,j)'*prt(:,j);
if nrm(j)>0,
nrm(j)=1/sqrt(nrm(j));
end;
end
prt(:,j)=nrm(j)*prt(:,j);
end;
% above loop could ? be replaced by:
% prt=prt.*wt(:,ones(1,n));
% nrm=dp./sqrt(diag(prt'*prt));
% prt=prt.*nrm(:,ones(1,m))';
[prt,s,v]=svd(prt,0);
s=diag(s);
g=prt'*r;
for jjj=1:length(epstab),
epsL = max(epsLlast*epstab(jjj),1e-7);
se=sqrt((s.*s)+epsL);
gse=g./se;
chg=((v*gse).*nrm);
%   check the change constraints and apply as necessary
ochg=chg;
for iii=1:n,
if (maxstep(iii)==Inf), break; end;
chg(iii)=max(chg(iii),-abs(maxstep(iii)*pprev(iii)));
chg(iii)=min(chg(iii),abs(maxstep(iii)*pprev(iii)));
end;
if (verbose(1) & any(ochg ~= chg)),
disp(['Change in parameter(s): ' ...
sprintf('%d ',find(ochg ~= chg)) 'were constrained']);
end;
aprec=abs(pprec.*pbest);       %---
% ss=scalar sum of squares=sum((wt.*(y-f))^2).
if (any(abs(chg) > 0.1*aprec)),%---  % only worth evaluating function if
p=chg+pprev;                       % there is some non-miniscule change
f=feval(F,x,p);
r=wt.*(y-f);
ss=r'*r;
if ss<sbest,
pbest=p;
fbest=f;
sbest=ss;
end;
if ss<=sgoal,
break;
end;
end;                          %---
end;
epsLlast = epsL;
if (verbose(2)),
eval(plotcmd);
end;
if ss<eps,
break;
end
aprec=abs(pprec.*pbest);
%  [aprec chg chgprev]
if (all(abs(chg) < aprec) & all(abs(chgprev) < aprec)),
kvg=1;
if (verbose(1)),
fprintf('Parameter changes converged to specified precision\n');
end;
break;
else
chgprev=chg;
end;
if ss>sgoal,
break;
end;
end;

% set return values
%
p=pbest;
f=fbest;
ss=sbest;
kvg=((sbest>sgoal)|(sbest<=eps)|kvg);
if kvg ~= 1 , disp(' CONVERGENCE NOT ACHIEVED! '), end;

% CALC VARIANCE COV MATRIX AND CORRELATION MATRIX OF PARAMETERS
% re-evaluate the Jacobian at optimal values
jac=feval(dFdp,x,f,p,dp,F);
msk = dp ~= 0;
n = sum(msk);           % reduce n to equal number of estimated parameters
jac = jac(:, msk);      % use only fitted parameters

%% following section is Ray Muzic's estimate for covariance and correlation
%% assuming covariance of data is a diagonal matrix proportional to
%% diag(1/wt.^2).
%% cov matrix of data est. from Bard Eq. 7-5-13, and Row 1 Table 5.1

if vernum(1) >= 4,
Q=sparse(1:m,1:m,(0*wt+1)./(wt.^2));  % save memory
Qinv=inv(Q);
else
Qinv=diag(wt.*wt);
Q=diag((0*wt+1)./(wt.^2));
end;
resid=y-f;                                    %un-weighted residuals
covr=resid'*Qinv*resid*Q/(m-n);                 %covariance of residuals
Vy=1/(1-n/m)*covr;  % Eq. 7-13-22, Bard         %covariance of the data

jtgjinv=inv(jac'*Qinv*jac);                     %argument of inv may be singular
covp=jtgjinv*jac'*Qinv*Vy*Qinv*jac*jtgjinv; % Eq. 7-5-13, Bard %cov of parm est
d=sqrt(abs(diag(covp)));
corp=covp./(d*d');

covr=diag(covr);                 % convert returned values to compact storage
stdresid=resid./sqrt(diag(Vy));  % compute then convert for compact storage
Z=((m-n)*jac'*Qinv*jac)/(n*resid'*Qinv*resid);

%%% alt. est. of cov. mat. of parm.:(Delforge, Circulation, 82:1494-1504, 1990
%%disp('Alternate estimate of cov. of param. est.')
%%acovp=resid'*Qinv*resid/(m-n)*jtgjinv

%Calculate R^2 (Ref Draper & Smith p.46)
%
r=corrcoef(y,f);
if (exist('OCTAVE_VERSION'))
r2=r^2;
else
r2=r(1,2).^2;
end

% if someone has asked for it, let them have it
%
if (verbose(2)), eval(plotcmd); end,
if (verbose(1)),
disp(' Least Squares Estimates of Parameters')
disp(p')
disp(' Correlation matrix of parameters estimated')
disp(corp)
disp(' Covariance matrix of Residuals' )
disp(covr)
disp(' Correlation Coefficient R^2')
disp(r2)
sprintf(' 95%% conf region: F(0.05)(%.0f,%.0f)>=
delta_pvec''*Z*delta_pvec',n,m-n)
Z
%   runs test according to Bard. p 201.
n1 = sum((f-y) < 0);
n2 = sum((f-y) > 0);
nrun=sum(abs(diff((f-y)<0)))+1;
if ((n1>10)&(n2>10)), % sufficent data for test?
zed=(nrun-(2*n1*n2/(n1+n2)+1)+0.5)/(2*n1*n2*(2*n1*n2-n1-n2)...
/((n1+n2)^2*(n1+n2-1)));
if (zed < 0),
prob = erfc(-zed/sqrt(2))/2*100;
disp([num2str(prob) '% chance of fewer than ' num2str(nrun) ' runs.']);
else,
prob = erfc(zed/sqrt(2))/2*100;
disp([num2str(prob) '% chance of greater than ' num2str(nrun) ' runs.']);
end;
end;
end;

% A modified version of Levenberg-Marquardt
% Non-Linear Regression program previously submitted by R.Schrager.
% This version corrects an error in that version and also provides
% an easier to use version with automatic numerical calculation of
% the Jacobian Matrix. In addition, this version calculates statistics
% such as correlation, etc....
%
% Version 3 Notes
% Errors in the original version submitted by Shrager (now called version 1)
% and the improved version of Jutan (now called version 2) have been corrected.
% Additional features, statistical tests, and documentation have also been
% included along with an example of usage.  BEWARE: Some the the input and
% output arguments were changed from the previous version.
%