getjac.m

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

%%%%%%%%%% start of function getjac.m %%%%%%%%%%
%
% Function to calculate the Jacobian matrix (Equation B.14)
%
% 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: March 31, 2000 (corrected the asymptotic appoximation
%		for beta < -blarge)
%
% dfe=getjac(lambda,g,upper,a,beta,d,noise,nearzero,...
%     large);
%
% Inputs
%     lambda       array containing the Lagrange multipliers
%     g            matrix of scaled kernel functions
%     upper        array containing upper limit of prior 
%                    distributions
%     a            array containing the vector, a
%     beta         array of Lagrange multipliers
%     d      array containing sampled concentrations
%     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;
%     nearzero     value below which the asymptotic 
%                    approximation to zero is used
%     large        value above which the asymptotic 
%                    approximation to infinity is used
%
% Outputs
%     dfe          array containing the Jacobian matrix

function dfe=getjac(lambda,g,upper,a,beta,d,noise,...
       nearzero,large);

% calcualte array sizes and define new arrays
[ndata,nt]=size(g);
dfe=zeros(ndata);
dmndan=zeros(nt,1);
b=zeros(nt,1);
norml=norm(lambda);
blarge=100.;

% calculate the b vector (page 58)
for k=1:nt
     aa=0.;
     for l=1:ndata
          aa=aa+lambda(l)*g(l,k);
     end
     b(k)=aa;
end

% calculate the derivative of mhat with respect to a or b
mau=a.*upper;
for j=1:nt
     if beta(j) < -blarge

          % Equation 3.17
%          dmndan(j)=5*upper(j)/(9*beta(j))-...
%               2*(1-2*b(j)*upper(j))/(9*beta(j)^2);
	   dmndan(j)=-2/(b(j)-beta(j))^2;

     elseif abs(mau(j)) < nearzero
          au=a(j)*upper(j);

          % Equation 3.15
          dmndan(j)=-upper(j)^2*(15.-15.*au+8*au^2)/...
               (180.-180*au+105*au^2);

     elseif abs(mau(j)) > large

          % Equation 3.16
          dmndan(j)=-a(j)^(-2);

     else

          % Equations 3.14 and B.14
          dmndan(j)=((a(j)^2)*(upper(j)^2)*exp(-mau(j))...
               -1+2*exp(-mau(j))-exp(-2*mau(j)))/...
               ((a(j)^2)*(1-exp(-mau(j)))^2);

     end
end

% calculate the derivative of F with respect to lambda 

errfactor=(sqrt(ndata)*noise(1)+norm(d)*noise(2))/norml;

% Equation B.13
for j=1:ndata
     for l=1:ndata
          fl=0;
          for n=1:nt
               fl=fl-g(j,n)*dmndan(n)*g(l,n);
          end
          dfe(j,l)=fl-errfactor*lambda(j)*...
               lambda(l)/(norml^2);
     end
end
dfe=dfe+errfactor*eye(ndata);

return
%%%%%%%%%%  end of function getjac.m  %%%%%%%%%%