about summary refs log tree commit diff
path: root/sourcecodes/bnt-master/BNT/examples/static/Zoubin/mfa_cl.m
blob: b90bab187f34bcd21b162693bf9d05d140cf2dbd (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
% function [lik, likv]=mfa_cl(X,Lh,Ph,Mu,Pi);
% 
% Calculates log likelihoods of a data set under a mixture of factor
% analysis model.
%
% X - data matrix
% Lh - factor loadings 
% Ph - diagonal uniquenesses matrix
% Mu - mean vectors
% Pi - priors
%
% lik - log likelihood of X 
% likv - vector of log likelihoods
% 
% If 0 or 1 output arguments requested, lik is returned. If 2 output
% arguments requested, [lik likv] is returned.

function [lik, likv]=mfa_cl(X,Lh,Ph,Mu,Pi);

N=length(X(:,1));
D=length(X(1,:));
K=length(Lh(1,:));
M=length(Pi);

if (abs(sum(Pi)-1) > 1e-6) 
  disp('ERROR: Pi should sum to 1');
  return;
elseif ((size(Lh) ~= [D*M K]) | (size(Ph) ~= [D 1]) | (size(Mu) ~= [M D]) ...
  | (size(Pi) ~= [M 1] & size(Pi) ~= [1 M]))   
  disp('ERROR in input matrix sizes');
  return;
end;  

tiny=exp(-744);
const=(2*pi)^(-D/2);

I=eye(K);
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*sum((XM.*Xk)'))'; 	
end;

Hsum=rsum(H); 				

likv=log(Hsum+(Hsum==0)*tiny);
lik=sum(likv);