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