精华区文章阅读

发信人: songbirds (songbird), 信区: Matlab
标  题: Levenberg-Marquardt算法
发信站: 哈工大紫丁香 (Sun May  2 11:44:43 2004), 站内信件

声明：借用

试用的时候注意文件换行的问题,这主要是网页显示导致的.

Levenberg-Marquardt 法非线性回归函数

function [f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2]=leasqr(x,y,pin,func,stol,
niter,wt,dp,dfdp,options)
%function[f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2]=
% leasqr(x,y,pin,{func,stol,niter,wt,dp,dfdp,options})
%
% Version 3.beta
% 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=1 or length(x)) of statistical weights. These should be set
% to be proportional to (sqrts of var(y))^-1; (That is, the covaraince
% 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=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)

% func=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,func)

% 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

% {}= optional parameters
% ss=scalar sum of squares=sum-over-i(wt(i)*(y(i)-f(i)))^2.

% 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.
% 3) add verbose option
% 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.
%
% Refrences:
% 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,''o'',x(:,1),f); shg';
if (sscanf(version,'%f') >= 4),
  global verbose
  plotcmd='plot(x(:,1),y,''o'',x(:,1),f); figure(gcf)';
end;

if(exist('verbose')~=1), verbose=1; end;
if (nargin <= 8), dfdp='dfdp'; end;
if (nargin <= 7), dp=.001*(pin*0+1); end; %DT
if (nargin == 6), wt=1.0; end;
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(func,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,func);

  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;
  [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 & any(ochg ~= chg)),
      disp(['Change in parameter(s): ' ...
         sprintf('%d ',find(ochg ~= chg)) 'were constrained']);
    end;
    aprec=abs(pprec.*pbest); %---
    if (any(abs(chg) > 0.1*aprec)),%--- % only worth evaluating function if

      p=chg+pprev; % there is some non-miniscule change
      f=feval(func,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),
    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),
      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,func);
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

Qinv=diag(wt.*wt);
Q=diag((0*wt+1)./(wt.^2));
%[nrw ncw]=size(wt);
%Q=ones(nrw,ncw)./wt; Q=diag(Q.*Q);
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
covr=diag(covr); %for compact storage
Z=((m-n)*jac'*Qinv*jac)/(n*resid'*Qinv*resid);
stdresid=resid./sqrt(diag(Vy));

jtgjinv=inv(jac'*Qinv*jac);
covp=jtgjinv*jac'*Qinv*Vy*Qinv*jac*jtgjinv; % Eq. 7-5-13, Bard %cov of parm es
t
for k=1:n,
  for j=k:n,
    corp(k,j)=covp(k,j)/sqrt(abs(covp(k,k)*covp(j,j)));
    corp(j,k)=corp(k,j);
  end;
end;

%%% 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);
r2=r(1,2).^2;

% if someone has asked for it, let them have it
%
if (verbose),
  eval(plotcmd);
  disp(' Least Squares Estimates of Parameters')
  disp(p')
  disp(' Correlation matrix of parameters estimated')
  disp(corp)
  disp('Covariance matriix 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
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, statisitcal 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.
%
% Ray Muzic rfm2@ds2.uh.cwru.edu
% Arthur Jutan jutan@charon.engga.uwo.ca

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
还需要另外两个函数，根据自己的情况改写

function y=leasqrfunc(t,p)
    c0=41;
    delta=0.006;
    v=0.000035;
    A=0.075^2*3.14;
     L=A/v; % Loading factor
    Q=0.0001/60;
    N=100/60/35; %s^-1
    Kma=p(1);
    Dm=p(2);
    h=0.001;
Bim=h*delta/Dm;
alfa=N*delta^2/Dm;
beta=L*delta;
xx=sqrt(alfa+beta*Bim);
% xx is not the root of the equation, but when x approaches xx, the denominato
r equals zero, which will cause problem when using function fzero
xxx=round(xx/pi) ;
syms x
T=x*tan(x)-(alfa-x^2)/[Kma*beta+(alfa-x^2)*Kma*Bim^-1];
    for i=1:200
       if(i==1)
           x0=[1e-10,pi/2-1e-10];
       else
           x0=[pi/2+1e-10,3*pi/2-1e-10]+(i-2)*pi;
           if (i==xxx+1)
               x0=[(xxx-1/2)*pi+1e-10,(xx-1e-10)];
           end
           if (i==xxx+2)
               x0=[(xx+1e-10),(xxx+1/2)*pi-1e-10];
           end
       end
       x(i)=fzero(inline(T),x0); % the egienvalues are correct
   end
    x=numeric(x);
          sum=0;
      for i=1:200
      An=[Kma*beta+(alfa-x(i)^2)*Kma*Bim^(-1)+2]*x(i)^2*cos(x(i))+x(i)*sin(x(i
))*[Kma*beta+(alfa-3*x(i)^2)*Kma*Bim^(-1)+alfa-x(i)^2];
      ca=2*c0*beta*x(i)*sin(x(i))*exp(-Dm*delta^(-2)*x(i)^2*t*3600)/An;
      sum=ca+sum;
    end
y=sum;
return

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

function y = leasqrdfdp(x,f,p,dp,func)
%here x stands for time. p : parameters. dp, func has been defined in
%leaqrexamp.m
epslon=1e-5;
dfdp1=(leasqrfunc(x,[p(1)+epslon*p(1)/2;p(2)])-leasqrfunc(x,[p(1)-epslon*p(1)/
2;p(2)]))/(epslon*p(1));
dfdp2=(leasqrfunc(x,[p(1);p(2)+epslon*p(2)/2])-leasqrfunc(x,[p(1);p(2)-epslon*
p(2)/2]))/(epslon*p(2));
y=[dfdp1,dfdp2];

＝＝＝＝＝＝
以及主调用函数

% leasqr example/test
clc
clear
global verbose
verbose=0;
t=[2.67,7.67,23.08,29.17,47.83,71.67,86.83,97.17,108.67,122.17,169.83,203.17,2
26.25,250,274.75,298.5]';
pin=[1000; 3e-13];
data=[0.004414916,0.003813583,0.002522687,0.002809648,0.002480928,0.002152208,
0.001721766,0.001579259,0.001488343,0.001375831,0.001227828,0.001363707,0.0011
72298,0.0011379,0.001103502,0.001104185]';
%wt1=(1+0*t)./sqrt(data); % = 1 /sqrt of variances of data
wt1=1;
% if (sscanf(version,'%f') < 4),
% rand('normal');
% data=data+0.05*rand(16,1)./wt1; % 1/wt1 = sqrt of var = standard deviation

% else
% data=data+0.05*randn(16,1)./wt1;
% end;
options=[[0.0001; 0.0001] [0.8; 0.8]];
[f1,p1,kvg1,iter1,corp1,covp1,covr1,stdresid1,Z1,r21]=...
  leasqr(t,data,pin,'leasqrfunc',1e-6,500,wt1,[1;1],...
  'leasqrdfdp',options);
%function [f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2]=leasqr(x,y,pin,func,stol
,niter,wt,dp,dfdp,options)

--
http://forum.jnnc.com/UploadFile/200351314291199524.gif
http://aswater.hit.edu.cn/images/baby3.gif
http://chinaxue.y365.com/tupian/qunwu.gif
※ 来源:．哈工大紫丁香 bbs.hit.edu.cn [FROM: 210.46.77.90]

Matlab 版 (精华区)