ffa.m

来自「贝叶斯matlab程序设计算法」· M 代码 · 共 76 行

% function [L,Ph,LL]=ffa(X,K,cyc,tol);
% 
% Fast Maximum Likelihood Factor Analysis using EM
%
% X - data matrix
% K - number of factors
% cyc - maximum number of cycles of EM (default 100)
% tol - termination tolerance (prop change in likelihood) (default 0.0001)
%
% L - factor loadings 
% Ph - diagonal uniquenesses matrix
% LL - log likelihood curve
%
% Iterates until a proportional change < tol in the log likelihood 
% or cyc steps of EM 
%

function [L,Ph,LL]=ffa(X,K,cyc,tol);

if nargin<4  tol=0.0001; end;
if nargin<3  cyc=100; end;

N=length(X(:,1));
D=length(X(1,:));
tiny=exp(-700);

X=X-ones(N,1)*mean(X);
XX=X'*X/N;
diagXX=diag(XX);

randn('seed', 0);
cX=cov(X);
scale=det(cX)^(1/D);
L=randn(D,K)*sqrt(scale/K);
Ph=diag(cX);

I=eye(K);

lik=0; LL=[];

const=-D/2*log(2*pi);


for i=1:cyc;

  %%%% E Step %%%%
  Phd=diag(1./Ph);
  LP=Phd*L;
  MM=Phd-LP*inv(I+L'*LP)*LP';
  dM=sqrt(det(MM));
  beta=L'*MM;
  XXbeta=XX*beta';
  EZZ=I-beta*L +beta*XXbeta;

  %%%% Compute log likelihood %%%%
  
  oldlik=lik;
  lik=N*const+N*log(dM)-0.5*N*sum(diag(MM*XX));
  fprintf('cycle %i lik %g \n',i,lik);
  LL=[LL lik];
  
  %%%% M Step %%%%

  L=XXbeta*inv(EZZ);
  Ph=diagXX-diag(L*XXbeta');

  if (i<=2)    
    likbase=lik;
  elseif (lik<oldlik)     
    disp('VIOLATION');
  elseif ((lik-likbase)<(1+tol)*(oldlik-likbase)|~finite(lik))  
    break;
  end;

end