view toolboxes/FullBNT-1.0.7/bnt/examples/static/Zoubin/mfa.m @ 0:e9a9cd732c1e tip

first hg version after svn

% function [Lh,Ph,Mu,Pi,LL]=mfa(X,M,K,cyc,tol);
% 
% Maximum Likelihood Mixture of Factor Analysis using EM
%
% X - data matrix
% M - number of mixtures (default 1)
% K - number of factors in each mixture (default 2)
% cyc - maximum number of cycles of EM (default 100)
% tol - termination tolerance (prop change in likelihood) (default 0.0001)
%
% Lh - factor loadings 
% Ph - diagonal uniquenesses matrix
% Mu - mean vectors
% Pi - priors
% LL - log likelihood curve
%
% Iterates until a proportional change < tol in the log likelihood 
% or cyc steps of EM 

function [Lh, Ph,  Mu, Pi, LL] = mfa(X,M,K,cyc,tol)

if nargin<5   tol=0.0001; end;
if nargin<4   cyc=100; end;
if nargin<3   K=2; end;
if nargin<2   M=1; end;

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

%rand('state',0);

fprintf('\n');

if (M==1)
  [Lh,Ph,LL]=ffa(X,K,cyc,tol);
  Mu=mean(X);
  Pi=1;
else
  if N==1
    mX = X;
  else
    mX=mean(X);
  end
  cX=cov(X);
  scale=det(cX)^(1/D);
  randn('state',0); 
  Lh=randn(D*M,K)*sqrt(scale/K);
  Ph=diag(cX)+tiny;
  Pi=ones(M,1)/M;
  %randn('state',0); 
  Mu=randn(M,D)*sqrtm(cX)+ones(M,1)*mX;
  oldMu=Mu;
  I=eye(K);

  lik=0;
  LL=[];

  H=zeros(N,M); 	% E(w|x) 
  EZ=zeros(N*M,K);
  EZZ=zeros(K*M,K);
  XX=zeros(D*M,D);
  s=zeros(M,1);
  const=(2*pi)^(-D/2);
  %%%%%%%%%%%%%%%%%%%%
  for i=1:cyc;

    %%%% E Step %%%%

    Phi=1./Ph;
    Phid=diag(Phi);
    for k=1:M
      Lht=Lh((k-1)*D+1:k*D,:);
      LP=Phid*Lht;
      MM=Phid-LP*inv(I+Lht'*LP)*LP';
      dM=sqrt(det(MM));      	
      Xk=(X-ones(N,1)*Mu(k,:)); 
      XM=Xk*MM;
      H(:,k)=const*Pi(k)*dM*exp(-0.5*rsum(XM.*Xk)); 	
      EZ((k-1)*N+1:k*N,:)=XM*Lht;
    end;
    
    Hsum=rsum(H);
    oldlik=lik;
    lik=sum(log(Hsum+(Hsum==0)*exp(-744)));

    Hzero=(Hsum==0); Nz=sum(Hzero); 
    H(Hzero,:)=tiny*ones(Nz,M)/M; 
    Hsum(Hzero)=tiny*ones(Nz,1);
    
    H=rdiv(H,Hsum); 				
    s=csum(H);
    s=s+(s==0)*tiny;
    s2=sum(s)+tiny;
    
    for k=1:M  
      kD=(k-1)*D+1:k*D;
      Lht=Lh(kD,:);
      LP=Phid*Lht;
      MM=Phid-LP*inv(I+Lht'*LP)*LP';
      Xk=(X-ones(N,1)*Mu(k,:)); 
      XX(kD,:)=rprod(Xk,H(:,k))'*Xk/s(k); 
      beta=Lht'*MM;
      EZZ((k-1)*K+1:k*K,:)=I-beta*Lht +beta*XX(kD,:)*beta'; 
    end;

    %%%% log likelihood %%%%

    LL=[LL lik];
    fprintf('cycle %g   \tlog likelihood %g ',i,lik);
    
    if (i<=2)
      likbase=lik;
    elseif (lik<oldlik) 
      fprintf(' violation');
    elseif ((lik-likbase)<(1 + tol)*(oldlik-likbase)||~isfinite(lik)) 
      break;
    end;

    fprintf('\n');
    
    %%%% M Step %%%%
    
    % means and covariance structure
    
    Ph=zeros(D,1);
    for k=1:M
      kD=(k-1)*D+1:k*D;
      kK=(k-1)*K+1:k*K;
      kN=(k-1)*N+1:k*N;

      T0=rprod(X,H(:,k));
      T1=T0'*[EZ(kN,:) ones(N,1)];
      XH=EZ(kN,:)'*H(:,k);
      T2=inv([s(k)*EZZ(kK,:) XH; XH' s(k)]);
      T3=T1*T2;
      Lh(kD,:)=T3(:,1:K);
      Mu(k,:)=T3(:,K+1)';
      T4=diag(T0'*X-T3*T1')/s2;
      Ph=Ph+T4.*(T4>0); 
    end;

    Phmin=exp(-700);
    Ph=Ph.*(Ph>Phmin)+(Ph<=Phmin)*Phmin; % to avoid zero variances

    % priors
    Pi=s'/s2;
    
  end;
  fprintf('\n');
end;