Daniel@0: % function [Lh,Ph,Mu,Pi,LL]=mfa(X,M,K,cyc,tol);
Daniel@0: % 
Daniel@0: % Maximum Likelihood Mixture of Factor Analysis using EM
Daniel@0: %
Daniel@0: % X - data matrix
Daniel@0: % M - number of mixtures (default 1)
Daniel@0: % K - number of factors in each mixture (default 2)
Daniel@0: % cyc - maximum number of cycles of EM (default 100)
Daniel@0: % tol - termination tolerance (prop change in likelihood) (default 0.0001)
Daniel@0: %
Daniel@0: % Lh - factor loadings 
Daniel@0: % Ph - diagonal uniquenesses matrix
Daniel@0: % Mu - mean vectors
Daniel@0: % Pi - priors
Daniel@0: % LL - log likelihood curve
Daniel@0: %
Daniel@0: % Iterates until a proportional change < tol in the log likelihood 
Daniel@0: % or cyc steps of EM 
Daniel@0: 
Daniel@0: function [Lh, Ph,  Mu, Pi, LL] = mfa(X,M,K,cyc,tol)
Daniel@0: 
Daniel@0: if nargin<5   tol=0.0001; end;
Daniel@0: if nargin<4   cyc=100; end;
Daniel@0: if nargin<3   K=2; end;
Daniel@0: if nargin<2   M=1; end;
Daniel@0: 
Daniel@0: N=length(X(:,1));
Daniel@0: D=length(X(1,:));
Daniel@0: tiny=exp(-700);
Daniel@0: 
Daniel@0: %rand('state',0);
Daniel@0: 
Daniel@0: fprintf('\n');
Daniel@0: 
Daniel@0: if (M==1)
Daniel@0:   [Lh,Ph,LL]=ffa(X,K,cyc,tol);
Daniel@0:   Mu=mean(X);
Daniel@0:   Pi=1;
Daniel@0: else
Daniel@0:   if N==1
Daniel@0:     mX = X;
Daniel@0:   else
Daniel@0:     mX=mean(X);
Daniel@0:   end
Daniel@0:   cX=cov(X);
Daniel@0:   scale=det(cX)^(1/D);
Daniel@0:   randn('state',0); 
Daniel@0:   Lh=randn(D*M,K)*sqrt(scale/K);
Daniel@0:   Ph=diag(cX)+tiny;
Daniel@0:   Pi=ones(M,1)/M;
Daniel@0:   %randn('state',0); 
Daniel@0:   Mu=randn(M,D)*sqrtm(cX)+ones(M,1)*mX;
Daniel@0:   oldMu=Mu;
Daniel@0:   I=eye(K);
Daniel@0: 
Daniel@0:   lik=0;
Daniel@0:   LL=[];
Daniel@0: 
Daniel@0:   H=zeros(N,M); 	% E(w|x) 
Daniel@0:   EZ=zeros(N*M,K);
Daniel@0:   EZZ=zeros(K*M,K);
Daniel@0:   XX=zeros(D*M,D);
Daniel@0:   s=zeros(M,1);
Daniel@0:   const=(2*pi)^(-D/2);
Daniel@0:   %%%%%%%%%%%%%%%%%%%%
Daniel@0:   for i=1:cyc;
Daniel@0: 
Daniel@0:     %%%% E Step %%%%
Daniel@0: 
Daniel@0:     Phi=1./Ph;
Daniel@0:     Phid=diag(Phi);
Daniel@0:     for k=1:M
Daniel@0:       Lht=Lh((k-1)*D+1:k*D,:);
Daniel@0:       LP=Phid*Lht;
Daniel@0:       MM=Phid-LP*inv(I+Lht'*LP)*LP';
Daniel@0:       dM=sqrt(det(MM));      	
Daniel@0:       Xk=(X-ones(N,1)*Mu(k,:)); 
Daniel@0:       XM=Xk*MM;
Daniel@0:       H(:,k)=const*Pi(k)*dM*exp(-0.5*rsum(XM.*Xk)); 	
Daniel@0:       EZ((k-1)*N+1:k*N,:)=XM*Lht;
Daniel@0:     end;
Daniel@0:     
Daniel@0:     Hsum=rsum(H);
Daniel@0:     oldlik=lik;
Daniel@0:     lik=sum(log(Hsum+(Hsum==0)*exp(-744)));
Daniel@0: 
Daniel@0:     Hzero=(Hsum==0); Nz=sum(Hzero); 
Daniel@0:     H(Hzero,:)=tiny*ones(Nz,M)/M; 
Daniel@0:     Hsum(Hzero)=tiny*ones(Nz,1);
Daniel@0:     
Daniel@0:     H=rdiv(H,Hsum); 				
Daniel@0:     s=csum(H);
Daniel@0:     s=s+(s==0)*tiny;
Daniel@0:     s2=sum(s)+tiny;
Daniel@0:     
Daniel@0:     for k=1:M  
Daniel@0:       kD=(k-1)*D+1:k*D;
Daniel@0:       Lht=Lh(kD,:);
Daniel@0:       LP=Phid*Lht;
Daniel@0:       MM=Phid-LP*inv(I+Lht'*LP)*LP';
Daniel@0:       Xk=(X-ones(N,1)*Mu(k,:)); 
Daniel@0:       XX(kD,:)=rprod(Xk,H(:,k))'*Xk/s(k); 
Daniel@0:       beta=Lht'*MM;
Daniel@0:       EZZ((k-1)*K+1:k*K,:)=I-beta*Lht +beta*XX(kD,:)*beta'; 
Daniel@0:     end;
Daniel@0: 
Daniel@0:     %%%% log likelihood %%%%
Daniel@0: 
Daniel@0:     LL=[LL lik];
Daniel@0:     fprintf('cycle %g   \tlog likelihood %g ',i,lik);
Daniel@0:     
Daniel@0:     if (i<=2)
Daniel@0:       likbase=lik;
Daniel@0:     elseif (lik<oldlik) 
Daniel@0:       fprintf(' violation');
Daniel@0:     elseif ((lik-likbase)<(1 + tol)*(oldlik-likbase)||~isfinite(lik)) 
Daniel@0:       break;
Daniel@0:     end;
Daniel@0: 
Daniel@0:     fprintf('\n');
Daniel@0:     
Daniel@0:     %%%% M Step %%%%
Daniel@0:     
Daniel@0:     % means and covariance structure
Daniel@0:     
Daniel@0:     Ph=zeros(D,1);
Daniel@0:     for k=1:M
Daniel@0:       kD=(k-1)*D+1:k*D;
Daniel@0:       kK=(k-1)*K+1:k*K;
Daniel@0:       kN=(k-1)*N+1:k*N;
Daniel@0: 
Daniel@0:       T0=rprod(X,H(:,k));
Daniel@0:       T1=T0'*[EZ(kN,:) ones(N,1)];
Daniel@0:       XH=EZ(kN,:)'*H(:,k);
Daniel@0:       T2=inv([s(k)*EZZ(kK,:) XH; XH' s(k)]);
Daniel@0:       T3=T1*T2;
Daniel@0:       Lh(kD,:)=T3(:,1:K);
Daniel@0:       Mu(k,:)=T3(:,K+1)';
Daniel@0:       T4=diag(T0'*X-T3*T1')/s2;
Daniel@0:       Ph=Ph+T4.*(T4>0); 
Daniel@0:     end;
Daniel@0: 
Daniel@0:     Phmin=exp(-700);
Daniel@0:     Ph=Ph.*(Ph>Phmin)+(Ph<=Phmin)*Phmin; % to avoid zero variances
Daniel@0: 
Daniel@0:     % priors
Daniel@0:     Pi=s'/s2;
Daniel@0:     
Daniel@0:   end;
Daniel@0:   fprintf('\n');
Daniel@0: end;
Daniel@0: 
Daniel@0: