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