Mercurial > hg > camir-aes2014
annotate toolboxes/FullBNT-1.0.7/bnt/examples/static/Zoubin/ffa.m @ 0:e9a9cd732c1e tip
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 [L,Ph,LL]=ffa(X,K,cyc,tol); |
| wolffd@0 | 2 % |
| wolffd@0 | 3 % Fast Maximum Likelihood Factor Analysis using EM |
| wolffd@0 | 4 % |
| wolffd@0 | 5 % X - data matrix |
| wolffd@0 | 6 % K - number of factors |
| wolffd@0 | 7 % cyc - maximum number of cycles of EM (default 100) |
| wolffd@0 | 8 % tol - termination tolerance (prop change in likelihood) (default 0.0001) |
| wolffd@0 | 9 % |
| wolffd@0 | 10 % L - factor loadings |
| wolffd@0 | 11 % Ph - diagonal uniquenesses matrix |
| wolffd@0 | 12 % LL - log likelihood curve |
| wolffd@0 | 13 % |
| wolffd@0 | 14 % Iterates until a proportional change < tol in the log likelihood |
| wolffd@0 | 15 % or cyc steps of EM |
| wolffd@0 | 16 % |
| wolffd@0 | 17 |
| wolffd@0 | 18 function [L,Ph,LL]=ffa(X,K,cyc,tol); |
| wolffd@0 | 19 |
| wolffd@0 | 20 if nargin<4 tol=0.0001; end; |
| wolffd@0 | 21 if nargin<3 cyc=100; end; |
| wolffd@0 | 22 |
| wolffd@0 | 23 N=length(X(:,1)); |
| wolffd@0 | 24 D=length(X(1,:)); |
| wolffd@0 | 25 tiny=exp(-700); |
| wolffd@0 | 26 |
| wolffd@0 | 27 X=X-ones(N,1)*mean(X); |
| wolffd@0 | 28 XX=X'*X/N; |
| wolffd@0 | 29 diagXX=diag(XX); |
| wolffd@0 | 30 |
| wolffd@0 | 31 randn('seed', 0); |
| wolffd@0 | 32 cX=cov(X); |
| wolffd@0 | 33 scale=det(cX)^(1/D); |
| wolffd@0 | 34 L=randn(D,K)*sqrt(scale/K); |
| wolffd@0 | 35 Ph=diag(cX); |
| wolffd@0 | 36 |
| wolffd@0 | 37 I=eye(K); |
| wolffd@0 | 38 |
| wolffd@0 | 39 lik=0; LL=[]; |
| wolffd@0 | 40 |
| wolffd@0 | 41 const=-D/2*log(2*pi); |
| wolffd@0 | 42 |
| wolffd@0 | 43 |
| wolffd@0 | 44 for i=1:cyc; |
| wolffd@0 | 45 |
| wolffd@0 | 46 %%%% E Step %%%% |
| wolffd@0 | 47 Phd=diag(1./Ph); |
| wolffd@0 | 48 LP=Phd*L; |
| wolffd@0 | 49 MM=Phd-LP*inv(I+L'*LP)*LP'; |
| wolffd@0 | 50 dM=sqrt(det(MM)); |
| wolffd@0 | 51 beta=L'*MM; |
| wolffd@0 | 52 XXbeta=XX*beta'; |
| wolffd@0 | 53 EZZ=I-beta*L +beta*XXbeta; |
| wolffd@0 | 54 |
| wolffd@0 | 55 %%%% Compute log likelihood %%%% |
| wolffd@0 | 56 |
| wolffd@0 | 57 oldlik=lik; |
| wolffd@0 | 58 lik=N*const+N*log(dM)-0.5*N*sum(diag(MM*XX)); |
| wolffd@0 | 59 fprintf('cycle %i lik %g \n',i,lik); |
| wolffd@0 | 60 LL=[LL lik]; |
| wolffd@0 | 61 |
| wolffd@0 | 62 %%%% M Step %%%% |
| wolffd@0 | 63 |
| wolffd@0 | 64 L=XXbeta*inv(EZZ); |
| wolffd@0 | 65 Ph=diagXX-diag(L*XXbeta'); |
| wolffd@0 | 66 |
| wolffd@0 | 67 if (i<=2) |
| wolffd@0 | 68 likbase=lik; |
| wolffd@0 | 69 elseif (lik<oldlik) |
| wolffd@0 | 70 disp('VIOLATION'); |
| wolffd@0 | 71 elseif ((lik-likbase)<(1+tol)*(oldlik-likbase)||~isfinite(lik)) |
| wolffd@0 | 72 break; |
| wolffd@0 | 73 end; |
| wolffd@0 | 74 |
| wolffd@0 | 75 end |