wolffd@0: % function [L,Ph,LL]=ffa(X,K,cyc,tol);
wolffd@0: % 
wolffd@0: % Fast Maximum Likelihood Factor Analysis using EM
wolffd@0: %
wolffd@0: % X - data matrix
wolffd@0: % K - number of factors
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: % L - factor loadings 
wolffd@0: % Ph - diagonal uniquenesses matrix
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: 
wolffd@0: function [L,Ph,LL]=ffa(X,K,cyc,tol);
wolffd@0: 
wolffd@0: if nargin<4  tol=0.0001; end;
wolffd@0: if nargin<3  cyc=100; end;
wolffd@0: 
wolffd@0: N=length(X(:,1));
wolffd@0: D=length(X(1,:));
wolffd@0: tiny=exp(-700);
wolffd@0: 
wolffd@0: X=X-ones(N,1)*mean(X);
wolffd@0: XX=X'*X/N;
wolffd@0: diagXX=diag(XX);
wolffd@0: 
wolffd@0: randn('seed', 0);
wolffd@0: cX=cov(X);
wolffd@0: scale=det(cX)^(1/D);
wolffd@0: L=randn(D,K)*sqrt(scale/K);
wolffd@0: Ph=diag(cX);
wolffd@0: 
wolffd@0: I=eye(K);
wolffd@0: 
wolffd@0: lik=0; LL=[];
wolffd@0: 
wolffd@0: const=-D/2*log(2*pi);
wolffd@0: 
wolffd@0: 
wolffd@0: for i=1:cyc;
wolffd@0: 
wolffd@0:   %%%% E Step %%%%
wolffd@0:   Phd=diag(1./Ph);
wolffd@0:   LP=Phd*L;
wolffd@0:   MM=Phd-LP*inv(I+L'*LP)*LP';
wolffd@0:   dM=sqrt(det(MM));
wolffd@0:   beta=L'*MM;
wolffd@0:   XXbeta=XX*beta';
wolffd@0:   EZZ=I-beta*L +beta*XXbeta;
wolffd@0: 
wolffd@0:   %%%% Compute log likelihood %%%%
wolffd@0:   
wolffd@0:   oldlik=lik;
wolffd@0:   lik=N*const+N*log(dM)-0.5*N*sum(diag(MM*XX));
wolffd@0:   fprintf('cycle %i lik %g \n',i,lik);
wolffd@0:   LL=[LL lik];
wolffd@0:   
wolffd@0:   %%%% M Step %%%%
wolffd@0: 
wolffd@0:   L=XXbeta*inv(EZZ);
wolffd@0:   Ph=diagXX-diag(L*XXbeta');
wolffd@0: 
wolffd@0:   if (i<=2)    
wolffd@0:     likbase=lik;
wolffd@0:   elseif (lik<oldlik)     
wolffd@0:     disp('VIOLATION');
wolffd@0:   elseif ((lik-likbase)<(1+tol)*(oldlik-likbase)||~isfinite(lik))  
wolffd@0:     break;
wolffd@0:   end;
wolffd@0: 
wolffd@0: end