donewton.m

基于李雅普诺夫的睡眠脑电动力学分析 核心算法程序。

%%%%%%%%%% start of function donewton.m %%%%%%%%%%
%
% Function to calculate the Lagrange multipliers, lambda, 
% using Newton's method
%
% Written by: Roseanna M. Neupauer
% Modification Date: April 24, 1999
% Modification Date: June 25, 1999 (added option to use additive noise,
%		multiplicative noise, or a combination of the two)
% Modification Date: February 12, 2000 (added feature to check if the initial
%		guess of lambda is a solution)
%
% [lambda,mhat,a]=...
%          donewton(funfcn,jacfcn,lamiter,tollam,tolls,...
%          g,upper,cin,noise,d,ndata,beta,...
%          nearzero,large);
% 
% Inputs
%     funfcn       name of matlab function that calculates the 
%                    F vector (Equation 2.28) funfcn='getfk'
%     jacfcn       name of matlab function that calculates the 
%                    Jacobian matrix (Equation B.13) 
%                    funfcn='getjac'
%     lamiter      maximum number of iterations in 
%                    lambda-fitting
%     tollam       tolerance for lambda-fitting
%     tolls        tolerance for golden section search
%     g            matrix of scaled kernel functions
%     upper        array containing upper limit of prior 
%                    distributions
%     t            array containing solution times
%     cin          array containing true source history
%     noise        vector of noise levels:
%			1. standard deviation of normally-distributed 
%                    	additive random noise in measurements;
%			2. standard deviation of normally-distributed 
%                    	multiplicative random noise in measurements;
%     d      array containing sampled concentrations
%     ndata        number of measured data points
%     beta         array of Lagrange multipliers
%     nearzero     value below which the asymptotic 
%                    approximation to zero is used
%     large        value above which the asymptotic 
%                    approximation to infinity is used
%
% Outputs
%     lambda       array containing the Lagrange multipliers
%     mhat         array containing the expected value of the
%                    model, based on the posterior distribution
%     a            array containing the vector, a
%
% Functions called
%     getfk        calculates the F vector
%     getjac       calculates the Jacobian matrix
%     plotops      creates some intermediate plots
%     rscale       row scales the Jacobian
%     golden       performs a golden section search

function [lambda,mhat,a]=...
          donewton(funfcn,jacfcn,lamiter,tollam,tolls,...
          g,upper,noise,d,ndata,beta,...
          nearzero,large);

k=0;
fignum=1;
normd=norm(d);
lambda=1*ones(ndata,1);      

% check if the initial guess satisfies the data constraint.

[fe,mhat,a]=feval(funfcn,lambda,g,beta,upper,...
          d,noise,nearzero,large);
ratio=norm(fe)/(1+normd);
if ratio > tollam 

   for i=1:lamiter

% calculate the F matrix (Equation 2.28)
     fprintf('Newton method: ITERATION %i\n',i)
     [fe,mhat,a]=feval(funfcn,lambda,g,beta,upper,d,...
          noise,nearzero,large);

% calculate the Jacobian matrix (Equation B.13)
     dfe=feval(jacfcn,lambda,g,upper,a,beta,d,noise,...
               nearzero,large);

% perform row scaling on the Jacobian matrix (page 59)
     [dfe,fe]=rscale(dfe,fe);
     fprintf('   Condition number of Jacobian: %e \n',cond(dfe))

% solve the matrix equation to obtain step direction 
% (Equations 3.13 and B.12)
     dellame=-dfe\fe;

% Use Golden section search to find the optimal step length 
% (page 59)
     alpha=golden(funfcn,lambda,dellame,tolls,g,beta,upper,...
          d,noise,nearzero,large);
     fprintf('   Step length is %f\n',alpha);

     lambda=lambda+alpha*dellame;
     [fe,mhat,a]=feval(funfcn,lambda,g,beta,upper,...
          d,noise,nearzero,large);
        fprintf('   new norm(F) is %e \n',norm(fe));

% check if convergence criteria is met (page 59)
     ratio=norm(fe)/(1+normd);
     fprintf('   Stopping criteria is %e \n',ratio);
     if ratio < tollam | isnan(ratio)
% calculate the fitted plume data
          break;
     end
   end
end
norml=norm(lambda);
fittedc=-fe+d+(sqrt(ndata)*noise(1)+normd*noise(2))*...
	lambda/norml;
if (i == lamiter & ratio > tollam)
	disp(' *** Newton method did not converge *** ')
end	
return
%%%%%%%%%%  end of function donewton.m  %%%%%%%%%%