domre.m

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

%%%%%%%%%% start of function domre.m %%%%%%%%%%
%
% Function to solves an inverse problem using Minimum Relative 
%     Entropy Inversion
%
% 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: July 12, 1999 (fixed the criteria for switching to
%		the asymptotic approximation for p5 and p95)
%		
% [x,lambda,beta,a,p5,p95] = ...
%                            domre(g,d,upper,expvalue,...
%                             noise,leftbegin,rightbegin,tolbeta,...
%                             maxiter,tollam,lamiter,tolls,nearzero,large);
%
% Inputs
%     d      array containing sampled concentrations
%     upper        array containing upper limit of prior 
%                    distributions
%     expvalue     array containing prior expected value
%     params       array of all transport parameters values
%     noise        noise level.
%     leftbegin    lower limit of range for bisection method 
%     rightbegin   upper limit of range for bisection method 
%     tolbeta      tolerance for beta-fitting
%     maxiter      maximum number of iterations in beta-fitting
%     tollam       tolerance for lambda-fitting
%     lamiter      maximum number of iterations in 
%                    lambda-fitting
%     tolls        tolerance for golden section search
%     nearzero     value below which the asymptotic 
%                    approximation to zero is used
%     large        value above which the asymptotic 
%                    approximation to infinity is used
%
% Outputs
%     x      array containing reconstructed solution of 
%                    source history 
%     lambda       array containing the Lagrange multipliers
%     beta         array containing the Lagrange multipliers
%     p5           array containing the 5th percentile 
%                    probability level
%     p95          array containing the 95th percentile 
%                    probability level
%
% Functions called
%     bisection    calculates the values of the Lagrange 
%                    multipliers, beta, using bisection method
%     getrep       calculates the values of the kernel matrix
%     donewton     calculates the values of the Lagrange 
%                    multipliers, lambda, using Newton's 
%                    method, and, in the process, calculates 
%                    the expected value of the posterior 
%                    distribution

function [x,lambda,beta,a,p5,p95] = ...
     domre(g,d,upper,expvalue,...
          noise,leftbegin,rightbegin,tolbeta,...
          maxiter,tollam,lamiter,tolls,nearzero,large);

% define size variables
     nt=size(g,2);
     ndata=size(d,1);

% calculate betas using bisection method  

     beta=bisection(upper,expvalue,leftbegin,rightbegin,...
          maxiter,tolbeta,nt,nearzero,large);

% use Newton-Raphson method to calculate values of Lagrange 
% multipliers, lambda, and, in the process, the expected value
% of the distribution

     [lambda,x,a]=...
          donewton('getfk','getjac',lamiter,tollam,tolls,...
          g,upper,noise,d,ndata,beta,...
          nearzero,large);

%
% calculate the probability levels of the distribution
% (page 36)
% NEED TO MODIFIY TO ACCOUNT FOR ASYMPTOTIC APPROX OF A
     p5=-log(0.05 .*(exp(-a.*upper)-1)+1)./a;
     p95=-log(0.95 .*(exp(-a.*upper)-1)+1)./a;
     alist=find(abs(a) < nearzero);
     p5(alist)=0.05*upper(alist);
     p95(alist)=0.95*upper(alist);

%%%%%%%%%%  end of function domre.m  %%%%%%%%%%