Mercurial > hg > camir-aes2014

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 % function [Lh,Ph,Mu,Pi,LL]=mfa(X,M,K,cyc,tol);
wolffd@0 2 %
wolffd@0 3 % Maximum Likelihood Mixture of Factor Analysis using EM
wolffd@0 4 %
wolffd@0 5 % X - data matrix
wolffd@0 6 % M - number of mixtures (default 1)
wolffd@0 7 % K - number of factors in each mixture (default 2)
wolffd@0 8 % cyc - maximum number of cycles of EM (default 100)
wolffd@0 9 % tol - termination tolerance (prop change in likelihood) (default 0.0001)
wolffd@0 10 %
wolffd@0 11 % Lh - factor loadings
wolffd@0 12 % Ph - diagonal uniquenesses matrix
wolffd@0 13 % Mu - mean vectors
wolffd@0 14 % Pi - priors
wolffd@0 15 % LL - log likelihood curve
wolffd@0 16 %
wolffd@0 17 % Iterates until a proportional change < tol in the log likelihood
wolffd@0 18 % or cyc steps of EM
wolffd@0 19
wolffd@0 20 function [Lh, Ph, Mu, Pi, LL] = mfa(X,M,K,cyc,tol)
wolffd@0 21
wolffd@0 22 if nargin<5 tol=0.0001; end;
wolffd@0 23 if nargin<4 cyc=100; end;
wolffd@0 24 if nargin<3 K=2; end;
wolffd@0 25 if nargin<2 M=1; end;
wolffd@0 26
wolffd@0 27 N=length(X(:,1));
wolffd@0 28 D=length(X(1,:));
wolffd@0 29 tiny=exp(-700);
wolffd@0 30
wolffd@0 31 %rand('state',0);
wolffd@0 32
wolffd@0 33 fprintf('\n');
wolffd@0 34
wolffd@0 35 if (M==1)
wolffd@0 36 [Lh,Ph,LL]=ffa(X,K,cyc,tol);
wolffd@0 37 Mu=mean(X);
wolffd@0 38 Pi=1;
wolffd@0 39 else
wolffd@0 40 if N==1
wolffd@0 41 mX = X;
wolffd@0 42 else
wolffd@0 43 mX=mean(X);
wolffd@0 44 end
wolffd@0 45 cX=cov(X);
wolffd@0 46 scale=det(cX)^(1/D);
wolffd@0 47 randn('state',0);
wolffd@0 48 Lh=randn(D*M,K)*sqrt(scale/K);
wolffd@0 49 Ph=diag(cX)+tiny;
wolffd@0 50 Pi=ones(M,1)/M;
wolffd@0 51 %randn('state',0);
wolffd@0 52 Mu=randn(M,D)*sqrtm(cX)+ones(M,1)*mX;
wolffd@0 53 oldMu=Mu;
wolffd@0 54 I=eye(K);
wolffd@0 55
wolffd@0 56 lik=0;
wolffd@0 57 LL=[];
wolffd@0 58
wolffd@0 59 H=zeros(N,M); % E(w|x)
wolffd@0 60 EZ=zeros(N*M,K);
wolffd@0 61 EZZ=zeros(K*M,K);
wolffd@0 62 XX=zeros(D*M,D);
wolffd@0 63 s=zeros(M,1);
wolffd@0 64 const=(2*pi)^(-D/2);
wolffd@0 65 %%%%%%%%%%%%%%%%%%%%
wolffd@0 66 for i=1:cyc;
wolffd@0 67
wolffd@0 68 %%%% E Step %%%%
wolffd@0 69
wolffd@0 70 Phi=1./Ph;
wolffd@0 71 Phid=diag(Phi);
wolffd@0 72 for k=1:M
wolffd@0 73 Lht=Lh((k-1)*D+1:k*D,:);
wolffd@0 74 LP=Phid*Lht;
wolffd@0 75 MM=Phid-LP*inv(I+Lht'*LP)*LP';
wolffd@0 76 dM=sqrt(det(MM));
wolffd@0 77 Xk=(X-ones(N,1)*Mu(k,:));
wolffd@0 78 XM=Xk*MM;
wolffd@0 79 H(:,k)=const*Pi(k)*dM*exp(-0.5*rsum(XM.*Xk));
wolffd@0 80 EZ((k-1)*N+1:k*N,:)=XM*Lht;
wolffd@0 81 end;
wolffd@0 82
wolffd@0 83 Hsum=rsum(H);
wolffd@0 84 oldlik=lik;
wolffd@0 85 lik=sum(log(Hsum+(Hsum==0)*exp(-744)));
wolffd@0 86
wolffd@0 87 Hzero=(Hsum==0); Nz=sum(Hzero);
wolffd@0 88 H(Hzero,:)=tiny*ones(Nz,M)/M;
wolffd@0 89 Hsum(Hzero)=tiny*ones(Nz,1);
wolffd@0 90
wolffd@0 91 H=rdiv(H,Hsum);
wolffd@0 92 s=csum(H);
wolffd@0 93 s=s+(s==0)*tiny;
wolffd@0 94 s2=sum(s)+tiny;
wolffd@0 95
wolffd@0 96 for k=1:M
wolffd@0 97 kD=(k-1)*D+1:k*D;
wolffd@0 98 Lht=Lh(kD,:);
wolffd@0 99 LP=Phid*Lht;
wolffd@0 100 MM=Phid-LP*inv(I+Lht'*LP)*LP';
wolffd@0 101 Xk=(X-ones(N,1)*Mu(k,:));
wolffd@0 102 XX(kD,:)=rprod(Xk,H(:,k))'*Xk/s(k);
wolffd@0 103 beta=Lht'*MM;
wolffd@0 104 EZZ((k-1)*K+1:k*K,:)=I-beta*Lht +beta*XX(kD,:)*beta';
wolffd@0 105 end;
wolffd@0 106
wolffd@0 107 %%%% log likelihood %%%%
wolffd@0 108
wolffd@0 109 LL=[LL lik];
wolffd@0 110 fprintf('cycle %g \tlog likelihood %g ',i,lik);
wolffd@0 111
wolffd@0 112 if (i<=2)
wolffd@0 113 likbase=lik;
wolffd@0 114 elseif (lik<oldlik)
wolffd@0 115 fprintf(' violation');
wolffd@0 116 elseif ((lik-likbase)<(1 + tol)*(oldlik-likbase)||~isfinite(lik))
wolffd@0 117 break;
wolffd@0 118 end;
wolffd@0 119
wolffd@0 120 fprintf('\n');
wolffd@0 121
wolffd@0 122 %%%% M Step %%%%
wolffd@0 123
wolffd@0 124 % means and covariance structure
wolffd@0 125
wolffd@0 126 Ph=zeros(D,1);
wolffd@0 127 for k=1:M
wolffd@0 128 kD=(k-1)*D+1:k*D;
wolffd@0 129 kK=(k-1)*K+1:k*K;
wolffd@0 130 kN=(k-1)*N+1:k*N;
wolffd@0 131
wolffd@0 132 T0=rprod(X,H(:,k));
wolffd@0 133 T1=T0'*[EZ(kN,:) ones(N,1)];
wolffd@0 134 XH=EZ(kN,:)'*H(:,k);
wolffd@0 135 T2=inv([s(k)*EZZ(kK,:) XH; XH' s(k)]);
wolffd@0 136 T3=T1*T2;
wolffd@0 137 Lh(kD,:)=T3(:,1:K);
wolffd@0 138 Mu(k,:)=T3(:,K+1)';
wolffd@0 139 T4=diag(T0'*X-T3*T1')/s2;
wolffd@0 140 Ph=Ph+T4.*(T4>0);
wolffd@0 141 end;
wolffd@0 142
wolffd@0 143 Phmin=exp(-700);
wolffd@0 144 Ph=Ph.*(Ph>Phmin)+(Ph<=Phmin)*Phmin; % to avoid zero variances
wolffd@0 145
wolffd@0 146 % priors
wolffd@0 147 Pi=s'/s2;
wolffd@0 148
wolffd@0 149 end;
wolffd@0 150 fprintf('\n');
wolffd@0 151 end;
wolffd@0 152
wolffd@0 153