Mercurial > hg > camir-aes2014
annotate toolboxes/FullBNT-1.0.7/bnt/examples/static/Zoubin/mfa_cl.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 [lik, likv]=mfa_cl(X,Lh,Ph,Mu,Pi); |
| wolffd@0 | 2 % |
| wolffd@0 | 3 % Calculates log likelihoods of a data set under a mixture of factor |
| wolffd@0 | 4 % analysis model. |
| wolffd@0 | 5 % |
| wolffd@0 | 6 % X - data matrix |
| wolffd@0 | 7 % Lh - factor loadings |
| wolffd@0 | 8 % Ph - diagonal uniquenesses matrix |
| wolffd@0 | 9 % Mu - mean vectors |
| wolffd@0 | 10 % Pi - priors |
| wolffd@0 | 11 % |
| wolffd@0 | 12 % lik - log likelihood of X |
| wolffd@0 | 13 % likv - vector of log likelihoods |
| wolffd@0 | 14 % |
| wolffd@0 | 15 % If 0 or 1 output arguments requested, lik is returned. If 2 output |
| wolffd@0 | 16 % arguments requested, [lik likv] is returned. |
| wolffd@0 | 17 |
| wolffd@0 | 18 function [lik, likv]=mfa_cl(X,Lh,Ph,Mu,Pi); |
| wolffd@0 | 19 |
| wolffd@0 | 20 N=length(X(:,1)); |
| wolffd@0 | 21 D=length(X(1,:)); |
| wolffd@0 | 22 K=length(Lh(1,:)); |
| wolffd@0 | 23 M=length(Pi); |
| wolffd@0 | 24 |
| wolffd@0 | 25 if (abs(sum(Pi)-1) > 1e-6) |
| wolffd@0 | 26 disp('ERROR: Pi should sum to 1'); |
| wolffd@0 | 27 return; |
| wolffd@0 | 28 elseif ((size(Lh) ~= [D*M K]) | (size(Ph) ~= [D 1]) | (size(Mu) ~= [M D]) ... |
| wolffd@0 | 29 | (size(Pi) ~= [M 1] & size(Pi) ~= [1 M])) |
| wolffd@0 | 30 disp('ERROR in input matrix sizes'); |
| wolffd@0 | 31 return; |
| wolffd@0 | 32 end; |
| wolffd@0 | 33 |
| wolffd@0 | 34 tiny=exp(-744); |
| wolffd@0 | 35 const=(2*pi)^(-D/2); |
| wolffd@0 | 36 |
| wolffd@0 | 37 I=eye(K); |
| wolffd@0 | 38 Phi=1./Ph; |
| wolffd@0 | 39 Phid=diag(Phi); |
| wolffd@0 | 40 for k=1:M |
| wolffd@0 | 41 Lht=Lh((k-1)*D+1:k*D,:); |
| wolffd@0 | 42 LP=Phid*Lht; |
| wolffd@0 | 43 MM=Phid-LP*inv(I+Lht'*LP)*LP'; |
| wolffd@0 | 44 dM=sqrt(det(MM)); |
| wolffd@0 | 45 Xk=(X-ones(N,1)*Mu(k,:)); |
| wolffd@0 | 46 XM=Xk*MM; |
| wolffd@0 | 47 H(:,k)=const*Pi(k)*dM*exp(-0.5*sum((XM.*Xk)'))'; |
| wolffd@0 | 48 end; |
| wolffd@0 | 49 |
| wolffd@0 | 50 Hsum=rsum(H); |
| wolffd@0 | 51 |
| wolffd@0 | 52 likv=log(Hsum+(Hsum==0)*tiny); |
| wolffd@0 | 53 lik=sum(likv); |
| wolffd@0 | 54 |