ica_adatap.m

ICA方法用于脑电图分析的程序

function [S,A,loglikelihood,Sigma,chi,exitflag]=ica_adatap(X,prior,par,draw);

% ICA_ADATAP  Mean field independent component analysis (ICA)
%    [S,A,LL,SIGMA,CHI,EXITFLAG]=ICA_ADATAP(X) performs linear 
%    instanteneous mixing ICA by solving the adaptive TAP or mean field 
%    or variational mean field equations [1-4].
%   
%    Input and output arguments: 
%      
%    X        : Mixed signal
%    A        : Estimated mixing matrix
%    S        : Estimated source mean values                                   
%    SIGMA    : Estimated noise covariance
%    LL       : Log likelihood, log P(X|A,SIGMA), divided by  
%               number of samples - mean field estimate
%    CHI      : Estimated source covariances
%    EXITFLAG : 0, no convergence; 1,  maximum number of iterations  
%               reached and 2, convergence criteria met.  
%
%    Outputs are sorted in descending order according to 'energy'
%    sum(A.^2,1))'.*sum(S.^2,2).
%
%    [S,A,LL,SIGMA,CHI]=ICA_ADATAP(X,PRIOR) allows for specification 
%    of priors on A, S, SIGMA and the method to be used. PRIOR.method 
%    overrides other choices of priors. PRIOR.A, PRIOR.S, Prior.Sigma 
%    and PRIOR.method are character strings that can take the following 
%    values
%    
%    PRIOR.A :            constant      constant mixing matrix 
%                         free          free likelihood optimization
%                         positive      constrained likelihood 
%                                       optimization (uses QUADPROG)
%                         MacKay        corresponding to MacKay's 
%                                       hyperparameter update
%
%    PRIOR.S :            bigauss       sum of two Gaussians, default
%                                       variance 1 centered at +/-1.  
%                         binary        binary +/-1 
%                         binary_01     binary 0/1             
%                         combi         combinations of other priors
%                                       E.g. M1 binary_01's and M2
%                                       exponential's are specified by
%                                       PRIOR.S1='binary_01';
%                                       PRIOR.M1=M1; 
%                                       PRIOR.S2='exponential'; 
%                                       PRIOR.M2=M2;
%                         exponential   exponential (positive)
%                         heavy_tail    heavy_tailed (non-analytic 
%  				        power law tail)
%                         Laplace       Laplace (double exponential)
%                         Gauss         Gaussian for factor analysis
%                                       PRIOR.Sigma='diagonal' and
%                                       probabilistic PCA 
%                                       PRIOR.Sigma='isotropic'
%
%    PRIOR.Sigma :        constant      constant noise covariance
%                         free          free likelihood optimization 
%                         isotropic     isotropic noise covariance
%                                       (scalar noise variance)
%                         diagonal      different noise variance for 
%                                       each sensor.
%
%    PRIOR.method :       constant      constant A and Sigma 
%                                       for test sets. A and Sigma
%                                       should initialized, see 
%                                       below. Sets prior.A='constant'
%                                       and prior.Sigma='constant'.
%                         fa            factor analysis (FA) sets
%                                       prior.A='free', 
%                                       prior.S='Gauss' and
%                                       prior.Sigma='diagonal'.
%                         neg_kurtosis  negative kurtosis set
%                                       prior.S='bigauss'.
%                         pos_kurtosis  positive kurtosis sets
%                                       prior.S='heavy_tail'. 
%                         positive      positive ICA sets
%                                       prior.S='exponential' and
%                                       prior.A='positive'.
%                         ppca          probabilistic PCA (pPCA)
%                                       sets prior.A='free', 
%                                       prior.S='Gauss' and 
%                                       prior.Sigma='isotropic'.
%
%    Default values are PRIOR.A='free', PRIOR.S='Laplace',
%    PRIOR.Sigma='isotropic' and PRIOR.method not set.
%
%    [S,A,LL,SIGMA,CHI]=ICA_ADATAP(X,PRIOR,PAR) is used for giving 
%    additional arguments. PAR is a structure array with the 
%    following fields 
%   
%       sources           Number of sources (default is quadratic
%                         size(X,1)). This parameter is overridden 
%                         by size(PAR.A_init,2) if PAR.A_init is
%                         defined.
%       solver            method for solving mean field equations
%                         can take the following character string
%                         values, default is sequential: 
%
%                         beliefprop2   Belief propagation -      
%                                       Sequential second order 
%                                       method by T. Minka. 
%                         sequential    normal sequential iterative
%                                       update (first order method) 
%                         
%       A_init            Initial mixing matrix (default is Toeplitz
%                         type with small randmom values).        
%       S_init            Initial source values (default is zero).
%       Sigma_init        Initial noise covariance (scalar for 
%                         isotropic). Default is Sigma_rel times
%                         empirical covariance.
%       Sigma_rel         See above, default 1 for PRIOR.A=
%                         'A_positive' and 100 otherwise. 
%       max_ite           Maximum of A and SIGMA updates (default
%                         is 20)
%       S_max_ite         Maximum number of S update steps in each
%                         E-step (default is 100)
%       tol               Termination tolerance for program in 
%                         terms of relative change in 1-norm of 
%                         A and Sigma (default is 10^-3).
%       S_tol             Termination tolerance for S updates in 
%                         terms of the squared error of the S
%                         mean field equations (default is 10^-10).
%
%    [S,A,LL,SIGMA,CHI]=ICA_ADATAP(X,PRIOR,PAR,DRAW) with DRAW 
%    different from zero will output runtime information. DRAW
%    equal to 2 will output the runtime value of the log 
%    likelihood. Default value is 1.
%
%    The memory requirements are O(M^2*N+D*N), where D is the 
%    number of sensors, N the number of samples (X is D*N) and M is 
%    the number of sources. The computational complexity is 
%    O(M^3*N) plus for positive mixing matrix, D times the 
%    computational complexity of a M-dimensional quadratic 
%    programming problem.   
%
%    It is recommended as a preprocessing step to normalize the 
%    data, e.g. by dividing by the largest element X=X/max(max(X)).

%
% References [1-4]:
%
%    Tractable Approximations for Probabilistic Models: The Adaptive 
%    Thouless-Anderson-Palmer Mean Field Approach
%    M. Opper and O. Winther
%    Phys. Rev. Lett. 86, 3695-3699 (2001).
%
%    Adaptive and Self-averaging Thouless-Anderson-Palmer Mean Field 
%    Theory for Probabilistic Modeling
%    M. Opper and O. Winther
%    Physical Review E 64, 056131 (2001). 
% 
%    Mean Field Approaches to Independent Component Analysis
%    Pedro A.d.F.R. H鴍en-S鴕ensen, Ole Winther and Lars Kai Hansen
%    Neural Computation 14, 889-918 (2002).
%  
%    TAP Gibbs Free Energy, Belief Propagation and Sparsity
%    L. Csato, M. Opper and O. Winther
%    In Advances in Neural Information Processing Systems 14 (NIPS'2001), 
%    MIT Press (2002).
 
% - by Ole Winther 2002 - IMM, Technical University of Denmark
% - http://isp.imm.dtu.dk/staff/winther/
% - version 2.2

% Uses adaptive TAP convention refs. [1,2] besides in the definition of the mean 
% function where it uses the ICA-convention of ref. [3]. 

[D,N]=size(X);                                                % number of sensor, samples
try prior=prior; catch prior=[]; end                          % make sure prior is defined
%try par=par; catch par=[]; end                               % make sure par is defined
[prior]=method_init(prior);                                   % initialize method
[A,A_cmd,A_prior,M]=A_init(prior,par,D);                      % initialize A
A_norm_old=norm(abs(A),1);                                    % quantify chance in A
[S,S_mean_cmd_all,S_mean_cmd_seq,likelihood_cmd,likelihood_arg,...
 S_arg_seq,cS_prior,cM,S_prior]=S_init(prior,par,M,N);        % initialize S
[Sigma,Sigma_cmd,Sigma_eta,Sigma_prior,Sigma_rel]=...
Sigma_init(prior,par,A_prior,S_prior,X,A,S);                  % initialize Sigma
dim_Sigma=size(Sigma,1)*size(Sigma,2);                        % number of elements in Sigma 
Sigma_norm_old=norm(abs(Sigma),1);                            % quantify chance in Sigma.
try max_ite=par.max_ite; catch max_ite=20; end                % default maximum number of EM-steps 
try S_max_ite=par.S_max_ite; catch S_max_ite=100; end         % default maximum number of S updates in each E-step 
try S_min_ite=par.S_min_ite; catch S_min_ite=1; end           % default minimum number of S updates in each E-step 
try tol=par.tol; catch tol=10^-3; end                         % termination tolerance relative chance in norm(Sigma)+norm(A) 
try S_tol=par.S_tol; catch S_tol=10^-10; end; S_tolN=S_tol/N; % termination tolerance chance in dS.*dS 
try                                                           % default solver is sequential
  switch par.solver
    case {'beliefprop2','convergent'}
      solver=par.solver;    
    otherwise
      solver='sequential';
  end
catch solver='sequential'; end
try draw=draw; catch draw=1; end                              % set draw 

% here starts the algorithm! 

if strcmp(solver,'beliefprop2')
  fprintf(' ICA Adaptive TAP - %s solver\n',solver); 
else  
  fprintf(' ICA Linear Response - %s solver\n',solver); 
end
fprintf(' %s source prior\n',S_prior); %s A and %s noise optimization\n',S_prior,A_prior,Sigma_prior); 
if strcmp(S_prior,'combi')
  for i=1:length(cM) 
    if cM(i)==1 fprintf('   1 %s source prior\n',cS_prior{i});  
    else fprintf('   %i %s source priors\n',cM(i),cS_prior{i}); end
  end
end
fprintf(' %s A and %s noise optimization\n',A_prior,Sigma_prior); 
fprintf(' sensors: D = %i\n sources: M = %i\n samples: N = %i\n\n',D,M,N); 

dS=zeros(M,N); S_old=zeros(M,N);
chi=zeros(M,M,N); diagchi=zeros(M,N);
G=zeros(M,M,N);
diagG=zeros(M,N);
hpred=zeros(M,1); V=zeros(M,N); dV=zeros(M,N); 
traceSS=zeros(M,M); tracechi=zeros(M,M);

dA_rel=Inf; dSigma_rel=Inf;

if draw 
  if dim_Sigma==D
    fprintf(' noise variance = %g \n\n',sum(Sigma)/D);
  else       
    fprintf(' noise variance = %g \n\n',trace(Sigma)/size(Sigma,1));
  end  
end

EM_step=0;

while EM_step<max_ite & (dA_rel>tol | dSigma_rel>tol) 

  EM_step=EM_step+1;

  if dim_Sigma==D
    InvSigma=1./Sigma;  
    for i=1:M  % set up matrices for M-step.
      for j=i:M J(i,j)=-(A(:,i))'*(InvSigma.*A(:,j)); J(j,i)=J(i,j); end
      for k=1:N h(i,k)=(A(:,i))'*(InvSigma.*X(:,k)); end
    end      
  else
    InvSigma=inv(Sigma); 
    J=-A'*InvSigma*A;
    h=A'*InvSigma*X;
  end

  if EM_step>1 % heuristic for initializing V so that V<V_0
    Vrel=repmat(diag(J)./diagJ,1,N); V=Vrel.*V;
  else 
    V=zeros(M,N);
  end
  diagJ=diag(J);  

  V_0=-repmat(diag(J),1,N);
  J=J-diag(diag(J)); 

  S_ite=0;  
  dSdS=Inf;
  dSdSN=Inf*ones(N,1);
  I=1:N;

  switch solver 

    case 'beliefprop2'

      alpha=S;     
      Omega=zeros(M,N);
      for i=1:N
        G(:,:,i)=J;
      end