ica_adatap.m

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

ccp=(xip>erfclimit);ccpc=not(ccp);
xip1=ccp.*xip;
xim=-(gamma+eta)./sqrt(lambda);
ccm=(xim>erfclimit);ccmc=not(ccm);
xim1=ccm.*xim;
Dp=exp(-(xip1.^2)/2)/sqrt(2*pi);
Dm=exp(-(xim1.^2)/2)/sqrt(2*pi);
kp=Phi(xip1).*Dm;  
km=Phi(xim1).*Dp; 
f=ccp.*ccm.*(xip.*kp-xim.*km)./(sqrt(lambda).*(kp+km))+(-ccpc.*xim+ccmc.*xip)./sqrt(lambda);
df=(ccp.*ccm.*(1+xim.*xip+Dp.*Dm.*(xip+xim)./(kp+km)+sqrt(lambda).*(xip.*km-xim.*kp)./(kp+km).*f)+ccpc+ccmc)./lambda;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                       logZ0 functions                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [logZ0] = logZ0_bigauss(gamma,lambda)
sigma2=1; mu=1;
lambda=1./(1+lambda*sigma2);
logZ0terms=0.5*(log(lambda)-mu^2/sigma2+lambda.*(mu^2/sigma2+gamma.^2*sigma2))...
           -log(2)+abs(gamma*mu.*lambda)+log(1+exp(-2*abs(gamma*mu.*lambda))); % log(cosh(gamma*mu.*lambda))
logZ0=sum(sum(logZ0terms));

function [logZ0] = logZ0_binary(gamma,lambda)
logZ0terms=-0.5*lambda-log(2)+abs(gamma)+log(1+exp(-2*abs(gamma)));
logZ0=sum(sum(logZ0terms));

function [logZ0] = logZ0_binary_01(gamma,lambda)
gamma=0.5*(gamma-0.5*lambda);
logZ0terms=gamma-log(2)+abs(gamma)+log(1+exp(-2*abs(gamma)));
logZ0=sum(sum(logZ0terms));

function [logZ0] = logZ0_combi(gamma,lambda,prior)
logZ0=0;
eval(prior);

function [logZ0] = logZ0_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;
logZ0terms=cc.*(log(Phi(xi1))+0.5*log(2*pi)+0.5*xi1.^2)-(1-cc).*log(abs(xi)+cc)+log(eta)-0.5*log(lambda);
logZ0=sum(sum(logZ0terms));

function [logZ0] = logZ0_Gauss(gamma,lambda)
logZ0terms=0.5*gamma.^2./(1+lambda)-0.5*log(1+lambda);
logZ0=sum(sum(logZ0terms));

function [logZ0] = logZ0_heavy_tail(gamma,lambda)
alpha=1; % if changed change also in heavy_tail.m
logZ0terms=0.5*gamma.^2./lambda-0.5*alpha*log(1+gamma.^2./(2*lambda*alpha));
logZ0=sum(sum(logZ0terms));

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

function [logZ0] = logZ0_Laplace(gamma,lambda)
erfclimit=-25;
eta=1;
minlambda=10^-4;
lambda=lambda.*(lambda>minlambda)+minlambda.*(lambda<=minlambda);
xip=(gamma-eta)./sqrt(lambda);
ccp=(xip>erfclimit);ccpc=not(ccp);
xip1=ccp.*xip;
xim=-(gamma+eta)./sqrt(lambda);
ccm=(xim>erfclimit);ccmc=not(ccm);
xim1=ccm.*xim;
Dp=exp(-(xip1.^2)/2)/sqrt(2*pi);
Dm=exp(-(xim1.^2)/2)/sqrt(2*pi);
Phip=Phi(xip1);  
Phim=Phi(xim1); 
logZ0terms=log(0.5*eta)+0.5*log(2*pi./lambda)+... 
ccp.*ccm.*(0.5*xip1.^2+0.5*xim1.^2+log(Dm.*Phip+Dp.*Phim)-0.5*log(2*pi))+...
ccp.*ccmc.*(0.5*xip1.^2+log(Phip+Dp./(abs(xim)+ccm)))+...
ccpc.*ccm.*(0.5*xim1.^2+log(Phim+Dm./(abs(xip)+ccp)))+...
ccpc.*ccmc.*(log(1./(abs(xip)+ccp)+1./(abs(xim)+ccm))-0.5*log(2*pi));
logZ0=sum(sum(logZ0terms));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                        Phi function                                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Phi function

function y=Phi(x)
%
z=abs(x/sqrt(2));
t=1.0./(1.0+0.5*z);
y=0.5*t.*exp(-z.*z-1.26551223+t.*(1.00002368+...
    t.*(0.37409196+t.*(0.09678418+t.*(-0.18628806+t.*(0.27886807+...
    t.*(-1.13520398+t.*(1.48851587+t.*(-0.82215223+t.*0.17087277)))))))));
y=(x<=0.0).*y+(x>0).*(1.0-y);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                   loglikelihood function                                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [loglikelihood]=loglikelihood_f(X,S,dfdg,h,hcav,V_0,V,J,G,likelihood_cmd,likelihood_arg,Sigma_prior,Sigma,InvSigma);
% calculates loglikelihood per S variable
%
[M N]=size(S);
D=size(X,1);
eval(likelihood_cmd);
free=free+sum(sum(hcav.*S))+0.5*sum(sum(V.*S.*S))+0.5*sum(sum(log(1+dfdg.*V)));
for i=1:N
  free=free-0.5*S(:,i)'*J*S(:,i)-0.5*log(det(squeeze(G(:,:,i))));
end
switch size(Sigma,1)*size(Sigma,2)
  case 1   % 'isotropic'
    logdetSigma=D*log(Sigma);  
    XInvSigmaX=InvSigma*sum(sum(X.^2));
  case D   % 'diagonal'
    logdetSigma=sum(log(Sigma)); 
    XInvSigmaX=InvSigma'*sum(X.^2,2);
  case D^2 % 'free'
    [L,U]=lu(Sigma);
    logdetSigma=sum(log(abs(diag(U)))); % assuming positive determinant
    XInvSigmaX=sum(sum(X.*(InvSigma*X)));
end
free=free+0.5*N*(D*log(2*pi)+logdetSigma)+0.5*XInvSigmaX;
loglikelihood=-free/N;

% initialization of method, A, S and Sigma

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                       Initialize method                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [prior,par]=method_init(prior,par)
try 
  switch prior.method
    case 'constant'
      prior.A='constant'; prior.Sigma='constant';
    case 'fa'
      prior.A='free'; prior.S='Gauss'; prior.Sigma='diagonal';
    case 'neg_kurtosis'
      prior.S='bigauss';
    case 'pos_kurtosis'
      prior.S='heavy_tail'; 
    case 'positive'
      prior.S='exponential'; prior.A='positive';
    case 'ppca' % probabilistic PCA
      prior.A='free'; prior.S='Gauss'; prior.Sigma='isotropic';
  end
end     

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                       Initialize A                                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [A,A_cmd,A_prior,M]=A_init(prior,par,D);
try A_prior=prior.A; catch A_prior='free'; end                % default mixing matrix
switch A_prior
  case 'constant'
    A_cmd='';  
  case 'MacKay'
    A_cmd=sprintf('A=A_free(traceSS-tracechi,X(:,I)*(S(:,I))''-A*tracechi,A);');
  otherwise
    A_cmd=sprintf('A=A_%s(traceSS,X(:,I)*(S(:,I))'',A);',A_prior);
end
try
  A=par.A_init;
  M=size(A,2);
catch
  try M=par.sources; catch M=D; end                           % quadratic mixing matrix is default 
  aMD=abs(D-M); scale=1; A=scale*randn(D,M)/10;                                        
  if (D>M)                                                    % default is asymmetric toeplitz plus Gaussian noise 
    a=(M+1:-1:M-D+2)/(M+1); A=A+scale*toeplitz(a.*(a>0),(1:M)<=1);
  else 
    a=(D+1:-1:D-M+2)/(D+1); A=A+scale*toeplitz((1:D)<=1,a.*(a>0));
  end
  if strcmp('positive',A_prior) A=A.*sign(A);
  else A=sign(rand(D,M)-0.5).*A; end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                       Initialize S                                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [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);
try S_prior=prior.S; catch S_prior='heavy_tail'; end          % default prior
switch S_prior                                                % set up combi prior
  case 'combi'
    cS_prior={};
    for i=1:50 % maximum different priors
      try S_arg=sprintf('cS_prior{i}=prior.S%i; cM(i)=prior.M%i;',i,i); eval(S_arg); catch break; end 
    end
    if i==1 
      S_prior='heavy_tail';  
      S_mean_cmd_seq=sprintf('[f,dfdg] = %s(h(cindx,I)+hcav,V_0(cindx,I)-V(cindx,I));',S_prior);
      S_mean_cmd_all=sprintf('[f,dfdg] = %s(h(:,I)+hcav,V_0(:,I)-V(:,I));',S_prior); 
      likelihood_cmd=sprintf('free = -logZ0_%s(h+hcav,V_0-V);',S_prior);
    end
    bindx=1; S_arg=''; S_arg_seq={}; likelihood_arg='';
    for j=1:i-1
     eindx=min(M,bindx+cM(j)-1);
     if bindx>M break; end
     if (j==i-1 | bindx+cM(j)>M) eindx=M; end
     S_arg=sprintf('%s[f(%i:%i,:),df(%i:%i,:)] = %s(gamma(%i:%i,:),lambda(%i:%i,:));\n',...
                   S_arg,bindx,eindx,bindx,eindx,cS_prior{j},bindx,eindx,bindx,eindx);
     likelihood_arg=sprintf('%slogZ0 = logZ0 + logZ0_%s(gamma(%i:%i,:),lambda(%i:%i,:));\n',...
                            likelihood_arg,cS_prior{j},bindx,eindx,bindx,eindx);
     for k=bindx:eindx S_arg_seq{k}=sprintf('[f,df] = %s(gamma,lambda);',cS_prior{j}); end
     bindx=bindx+cM(j); 
    end  
    S_mean_cmd_seq=sprintf('[f,dfdg] = combi(h(cindx,I)+hcav,V_0(cindx,I)-V(cindx,I),S_arg_seq{cindx});');
    S_mean_cmd_all=sprintf('[f,dfdg] = combi(h(:,I)+hcav,V_0(:,I)-V(:,I),S_arg);'); 
    likelihood_cmd=sprintf('free = -logZ0_combi(h+hcav,V_0-V,likelihood_arg);');
  otherwise  
    likelihood_arg=''; cS_prior=''; cM=0; S_arg_seq='';
    S_mean_cmd_seq=sprintf('[f,dfdg] = %s(h(cindx,I)+hcav,V_0(cindx,I)-V(cindx,I));',S_prior);
    S_mean_cmd_all=sprintf('[f,dfdg] = %s(h(:,I)+hcav,V_0(:,I)-V(:,I));',S_prior); 
    likelihood_cmd=sprintf('free = -logZ0_%s(h+hcav,V_0-V);',S_prior);
end
try
  S=par.S_init;
catch
  S=zeros(M,N);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                       Initialize Sigma                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Sigma,Sigma_cmd,Sigma_eta,Sigma_prior,Sigma_rel]=Sigma_init(prior,par,A_prior,S_prior,X,A,S)  
try Sigma_eta=par.Sigma_eta; catch Sigma_eta=1; end           % learning rate for Sigma
try 
  Sigma_prior=prior.Sigma;
catch
  Sigma_prior='isotropic';                                      % default prior
end
switch Sigma_prior
  case 'constant' 
    Sigma_cmd='';
  case 'isotropic'
    Sigma_cmd=sprintf('Sigma = (1-Sigma_eta)*Sigma+Sigma_eta*(sum(sum((X(:,I)-A*S(:,I)).^2))+sum(sum((A*tracechi).*A)))/(NI*D);');    
  case 'diagonal'
    Sigma_cmd=sprintf('Sigma = (1-Sigma_eta)*Sigma+Sigma_eta*(sum((X(:,I)-A*S(:,I)).^2,2)+sum((A*tracechi).*A,2))/NI;');    
  case 'free'
    Sigma_cmd=sprintf('Sigma = (1-Sigma_eta)*Sigma+Sigma_eta*((X(:,I)-A*S(:,I))*(X(:,I)-A*S(:,I))''+A*tracechi*A'')/NI;');      
end
if strcmp('free',Sigma_prior) & strcmp('positive',A_prior)
  fprintf('   Prior analyse: uses isotropic noise approximation in A-update.\n   Full update is computationally expensive. Not implemented yet.');   
end
try 
  Sigma_rel=par.Sigma_rel;
catch
  if strcmp(S_prior,'heavy_tail') % | strcmp(S_prior,'heavy_tail_plus_delta')
    Sigma_rel=1;                                            % default relative noise level 
  else
    Sigma_rel=1;
  end
end 
[D,N]=size(X);
try
  Sigma=par.Sigma_init;  
catch                               
  switch Sigma_prior
    case {'isotropic','constant'}
      Sigma=Sigma_rel*sum(sum((X-A*S).^2))/(D*N);   
    case 'diagonal'
      Sigma=Sigma_rel*sum((X-A*S).^2,2)/N; 
    case 'free'
      Sigma=Sigma_rel*(X-A*S)*(X-A*S)'/N; 
  end
end