ica_adatap.m

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

      while (~isempty(I) & dSdS>S_tol & S_ite<S_max_ite)

        S_ite=S_ite+1;
  
        indx=randperm(M);

        for sindx=1:M

          cindx=indx(sindx); % current index 

          % update hcav and V
          if M==1
            hpred=(squeeze(G(1,1,I)))'.*alpha(1,I);
          else
            hpred=sum(shiftdim(G(cindx,:,I)).*alpha(:,I),1);
          end
          Gc=squeeze(G(cindx,cindx,I))';
          dV(cindx,I)=Gc./(1+Gc.*Omega(cindx,I))-V(cindx,I);
          if S_ite>1
            Vfrac=0.7; Vind=(abs(dV(cindx,I))<=Vfrac*(V_0(cindx,I)-V(cindx,I)));
            dV(cindx,I)=Vind.*dV(cindx,I)+Vfrac*(1-Vind).*sign(dV(cindx,I)).*(V_0(cindx,I)-V(cindx,I)); 
            V(cindx,I)=V(cindx,I)+dV(cindx,I);
          end 
          hcav=hpred-V(cindx,I).*(Omega(cindx,I).*hpred+alpha(cindx,I)); 
       
          % update S and derivative
          S_old=S(cindx,I);
          eval(S_mean_cmd_seq); S(cindx,I)=f;  % e.g [f,dfdg] = exponential(h(cindx,I)+hcav,V_0(cindx,I)-V(cindx,I));
          dS(cindx,I)=S(cindx,I)-S_old;

          % update alpha and Omega
          alpha_old=alpha(cindx,I); alpha(cindx,I)=(S(cindx,I)-dfdg.*hcav)./(1+dfdg.*V(cindx,I));
          Omega_old=Omega(cindx,I); Omega(cindx,I)=dfdg./(1+dfdg.*V(cindx,I));                   
          
          % update G 
          for j=1:length(I)
            i=I(j); dOmega=Omega(cindx,i)-Omega_old(j); sG=squeeze(G(:,:,i));
            if any((diag(sG)+dOmega/(1-dOmega*sG(cindx,cindx)-eps)*sG(:,cindx).^2)<zeros(M,1)) 
              dOmega=0; Omega(cindx,i)=Omega_old(j)+dOmega;
              alpha(cindx,i)= 0.5*(alpha_old(j)+alpha(cindx,i));
            end
            if (dOmega~=0) % use Sherman-Woodbury
              b = sG(:,cindx); G(:,:,i)=sG+dOmega/(1-dOmega*sG(cindx,cindx))*b*b'; 
            end
          end
 
        end % over variables - sindx

        dSdSN=sum(dS.*dS,1);
        I=find(dSdSN > S_tolN);
        dSdS=sum(dSdSN);

      end; % epoch loop - S_ite 
 
    case 'sequential'
     
      while (~isempty(I) & dSdS>S_tol & S_ite<S_max_ite)     

        S_ite=S_ite+1;  

        indx=randperm(M);

        for sindx=1:M

          cindx=indx(sindx); % current index 

          % update hcav           
          hcav=J(cindx,:)*S(:,I)-V(cindx,I).*S(cindx,I);
   
          % update S and derivative
          S_old=S(cindx,I);
          eval(S_mean_cmd_seq);  S(cindx,I)=f;
          dS(cindx,I)=S(cindx,I)-S_old;        
 
        end   
 
        dSdSN=sum(dS.*dS,1);
        I=find(dSdSN > S_tolN);         
        dSdS=sum(dSdSN);

      end      

  end % switch solver

  % calculate G
  switch solver 
    case 'beliefprop2' 
      for i=1:N
        diagOmega=diag(Omega(:,i)); 
        G(:,:,i)=eye(M)+diagOmega*squeeze(G(:,:,i));
      end
    otherwise
      I=1:N; hcav=J*S-V.*S; eval(S_mean_cmd_all); Omega=dfdg./(1+dfdg.*V);
      for i=1:N
        diagOmega=diag(Omega(:,i)); 
        Gi=inv(eye(M)-diagOmega*J); diagG(:,i)=diag(Gi); G(:,:,i)=Gi;  
      end
  end  

  % calculate chi
  for i=1:N
    diagOmega=diag(Omega(:,i));
    chi(:,:,i)=squeeze(G(:,:,i))*diagOmega;
  end 

  if draw
    hcav=J*S-V.*S; I=1:N; eval(S_mean_cmd_all);     
    for i=1:N
      diagchi(:,i)=diag(squeeze(chi(:,:,i)));
    end    
    zerochi=(dfdg<eps & diagchi<eps);
    dVdV=sum(sum(((1-zerochi).*(1./(dfdg+zerochi)-1./(diagchi+zerochi))).^2));
    fprintf(' Iteration %d\n E-step: S_ite = %d   dSdS = %g   dVdV = %g\n',EM_step,S_ite,dSdS,dVdV);    
    if draw==2
      loglikelihood=loglikelihood_f(X,f,dfdg,h,hcav,V_0,V,J,G,likelihood_cmd,...
                                    likelihood_arg,Sigma_prior,Sigma,InvSigma);
      fprintf(' loglikelihood = %6.3f\n',loglikelihood);
    end
  end

  find(dSdSN < S_tol);
  if ~isempty(I) 
    NI=length(I);
    for i=1:M  % set up matrices for M-step.
      for j=i:M
        tracechi(i,j)=sum(squeeze(chi(i,j,I)));
        tracechi(j,i)=tracechi(i,j);
        traceSS(i,j)=sum(squeeze(chi(i,j,I))+S(i,I)'.*S(j,I)');
        traceSS(j,i)=traceSS(i,j);
      end;
    end;
    eval(A_cmd);     % e.g. A=A_positive(traceSS,X(:,I)*(S(:,I))',A);   this is correct EM as long as A is independent of Sigma.
    eval(Sigma_cmd);
  else  
    disp('No samples converged - exiting'); exitflag=0; break
  end 

  % termination criteria 
  A_norm=norm(abs(A),1); dA_rel=abs(A_norm-A_norm_old)/A_norm; A_norm_old=A_norm;
  Sigma_norm=norm(abs(Sigma),1); dSigma_rel=abs(Sigma_norm-Sigma_norm_old)/Sigma_norm; Sigma_norm_old=Sigma_norm;

  if draw
    fprintf(' M-step: d|A|/|A| = %g   d|Sigma|/|Sigma| = %g\n',dA_rel,dSigma_rel);
    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

  % update parameters of the prior 

end % EM_step

% calculate loglikelihood
hcav=J*S-V.*S; I=1:N; eval(S_mean_cmd_all);
loglikelihood=loglikelihood_f(X,f,dfdg,h,hcav,V_0,V,J,G,likelihood_cmd,...
                              likelihood_arg,Sigma_prior,Sigma,InvSigma); 
if dim_Sigma==D 
  fprintf(' %d EM steps, loglikelihood = %g, noise variance = %g\n',...
          EM_step,loglikelihood,sum(Sigma)/D);   
else       
  fprintf(' %d EM steps, loglikelihood = %g, noise variance = %g\n',...
          EM_step,loglikelihood,trace(Sigma)/size(Sigma,1)); 
end   
loglikelihood=real(loglikelihood); % small imaginary contributions might occur

% return variables sorted according to 'energy'  
[Y,I]=sort((sum(A.^2,1))'.*sum(S.^2,2));
I=flipud(I);
S=S(I,:); A=A(:,I); 
for i=1:N
  chi(:,:,i)=squeeze(chi(I,I,i));
end 
if EM_step==max_ite exitflag=1;
elseif (dA_rel<tol & dSigma_rel<tol) exitflag=2; end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                        end of main program                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                          help functions                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                    mixing matrix (A) functions                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [A] = A_free(traceSS,XSt,A) 
A=XSt*inv(traceSS);

function [A] = A_positive(traceSS,XSt,A) 

A=A+10^-3*(A<eps);
amp=XSt./(A*traceSS);
Aneg=any(any(amp<eps));
KT_max_ite=200; KT_ite=0;
while ~Aneg & KT_ite<KT_max_ite 
  KT_ite=KT_ite+1;
  A=A.*amp;
  amp=XSt./(A*traceSS);
  Aneg=any(any(amp<eps));
end
if Aneg % use quadratic programming instead
  M=size(traceSS,1);
  D=size(XSt,1);
  options=optimset('Display','off','TolX',10^5*eps,'TolFun',10^5*eps);
  for i=1:D
    B=quadprog(traceSS,-XSt(i,:)',[],[],[],[],zeros(M,1),[],A(i,:)',options);
    A(i,:)=B';
  end 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                       mean (S) functions                                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [f,df] = bigauss(gamma,lambda)
sigma2=1; mu=1;
lambda=1./(1+lambda*sigma2);
f=tanh(mu*gamma.*lambda);
df=lambda.*(sigma2+mu^2.*lambda.*(1-f.^2));
f=lambda.*(sigma2*gamma+mu*f);

function [f,df] = binary(gamma,lambda)
f=tanh(gamma);
df=1-f.^2;

function [f,df] = binary_01(gamma,lambda)
tmp=exp(0.5*lambda+gamma);
f=1./(1+exp(0.5*lambda-gamma));
df=f.*(1-f);

function [f,df] = combi(gamma,lambda,prior)
eval(prior);

function [f,df] = exponential(gamma,lambda)
eta=1;
erfclimit=-35;
minlambda=10^-4;
lambda=lambda.*(lambda>minlambda)+minlambda.*(lambda<=minlambda);
xi=(gamma-eta)./sqrt(lambda);
cc=(xi>erfclimit);
xi1=xi.*cc;
epfrac=exp(-(xi1.^2)/2)./(Phi(xi1)*sqrt(2*pi));
f=cc.*(xi1+epfrac)./sqrt(lambda);            % need to go to higher order to get asymptotics right
df=(1./lambda+f.*(xi1./sqrt(lambda)-f)); % need to go to higher order to get asymptotics right -fix at some point!!!

function [f,df] = Gauss(gamma,lambda)
f=gamma./(1+lambda);
df=1./(1+lambda);

function [f,df] = heavy_tail(gamma,lambda)
alpha=1; % if changed change also in logZ0_heavy_tail
f=gamma./lambda-alpha*gamma./(2*alpha*lambda+gamma.^2);
df=1./lambda+alpha*(gamma.^2-2*alpha*lambda)./(2*alpha*lambda+gamma.^2).^2;

function [f,df] = heavy_tail_plus_delta(gamma,lambda)
alpha=1; % if changed change also in logZ0_heavy_tail_plus_delta
beta=0.3; % proporation delta if changed change also in logZ0_heavy_tail_plus_delta
Z0ht=exp(0.5*gamma.^2./lambda).*(1+gamma.^2./(2*alpha*lambda)).^(-0.5*alpha);
f=(1-beta)*(gamma./lambda-alpha*gamma./(2*alpha*lambda+gamma.^2))./...
  (beta./Z0ht+(1-beta));
df=(1-beta)./(beta./Z0ht+(1-beta)).*...
   (1./lambda+alpha*(gamma.^2-2*alpha*lambda)./(2*alpha*lambda+gamma.^2).^2+...
   (gamma./lambda-alpha*gamma./(2*alpha*lambda+gamma.^2)).^2)-f.^2;

function [f,df] = Laplace(gamma,lambda)
erfclimit=-25;
eta=1;
minlambda=10^-4;
lambda=lambda.*(lambda>minlambda)+minlambda.*(lambda<=minlambda);
xip=(gamma-eta)./sqrt(lambda);