www.pudn.com > flexica.zip > flexica.m

% 
% 
%   flexica.m 
% 
%		seungjin@postech.ac.kr 
% 
%		Last updated: 
%		February 6, 2003 
% 
%		function for batch version of flexible ICA 
% 
%		function [W]=flexica(x,n,nb,ga_W,maxits) 
% 
%		W: demixing matrix 
%		x: multivariate input data matrix 
%			m by N where m is the number of sensors and 
%			N is the number of data points 
%		n: number of sources (default=m) 
%     nb: number of blocks (default=5) 
%		ga_W: learning rate (default=.07) 
%		maxits: number of maximal iterations (default=100) 
% 
%		You can find independent components by y=Wx. 
% 
%     Tips: 1. Better to choose the number of blocks such that 
%              the number of data points in the block is about 3000. 
%           2. Sometimes you might want to increase maxits. 
%   
% 
 
 
function [W]=flexica(x,n,nb,ga_W,maxits); 
 
if nargin==0;  
   fprintf('function [W]=flexica(x,n,nb,ga_W,maxits) \n'); 
   break;  
end; 
if nargin<5; maxits=100; end; 
if nargin<4; ga_W=.07; end; 
if nargin<3; nb=5; end; 
if nargin<2; n=length(x(:,1)); end; 
 
 
% figure out the size of data matrix x 
m=length(x(:,1)); 
N=length(x(1,:)); 
 
 
% compute the block size 
bsize=floor((2*N)/(nb+1)); 
if mod(bsize,2) ~= 0 
   bsize=bsize-1; 
end 
 
 
% preprocessing: data sphering 
for i=1:m 
   x(i,:)=x(i,:)-mean(x(i,:)); 
end 
Rx=cov(x'); 
[u d v]=svd(Rx); 
Q=sqrt(inv(d(1:n,1:n)))*u(:,1:n)'; 
tempx=Q*x; 
clear x; 
x=tempx; 
clear tempx; 
 
 
% initial conditions 
W=eye(n,n); 
I=eye(n,n); 
nu=2*ones(n,1); 
loop=1; 
dellik=1; 
tol=0.0005; 
 
% initial likelihood 
yy=W*x; 
templik=0; 
for i=1:n 
   templik=templik+log(nu(i))-log(2*gamma(1/nu(i)))-mean(abs(yy(i,:)).^nu(i)); 
end 
lik=log(det(W))+templik; 
oldlik=lik; 
oldoldlik=lik;  
 
 
% main algorithm 
% big loop (sweep) 
while (dellik > tol) & (loop <= maxits), 
 
   oldoldlik=oldlik; 
   oldlik=lik; 
 
   % small loop depending on the number of blocks 
   for k=1:nb 
     
      itsnum=(loop-1)*(nb-1)+k; 
      y=W*x(:,(k-1)*bsize/2+1:(k+1)*bsize/2); 
 
      % nonlinear function 
      for i=1:n 
         fy(i,:)=sign(y(i,:)).*abs(y(i,:)).^(nu(i)-1); 
      end 
 
      % update W 
      scalingW=1/(n*bsize)*abs(trace(fy*y')); 
      change=1/N*(y*y'+1/(1+ga_W*scalingW)*(fy*y'-y*fy')); 
      W=W+ga_W*( I-change )*W; 
 
   end 
 
   % update nonlinear functions 
   yy=W*x; 
   % estimate moments 
   m2=mean(yy.^2'); 
   m4=mean(yy.^4'); 
   kurt=m4./(m2.^2)-3*ones(1,n); 
   % choose Gaussian exponent 
   for i=1:n 
      if kurt(i) > 20 
	      nu(i)=.8; 
	   elseif kurt(i) > 5 
         nu(i)=1; 
	   elseif kurt(i) > 0 
         nu(i)=1.3; 
	   else 
	      nu(i)=4; 
      end 
   end 
   % likelihood calculation 
   templik=0; 
   for i=1:n 
      templik=templik+log(nu(i))-log(2*gamma(1/nu(i)))-mean(abs(yy(i,:)).^nu(i)); 
   end 
   loop=loop+1; 
   lik=log(det(W))+templik; 
   dellik=(abs(lik-oldlik)+abs(oldlik-oldoldlik))/2; 
 
   % anealing learning rate 
   if dellik<.01 
      ga_W=.05; 
   elseif dellik<.005 
      ga_W=.03 
   else 
      ga_W=.07; 
   end 
 
   % display 
   fprintf('sweep= %d, difference of likelihood= %f\n', loop-1, dellik); 
 
end % end of big loop 
 
 
% demixing matrix 
W=W*Q;