wolffd@0: % function [lik, likv]=mfa_cl(X,Lh,Ph,Mu,Pi);
wolffd@0: % 
wolffd@0: % Calculates log likelihoods of a data set under a mixture of factor
wolffd@0: % analysis model.
wolffd@0: %
wolffd@0: % X - data matrix
wolffd@0: % Lh - factor loadings 
wolffd@0: % Ph - diagonal uniquenesses matrix
wolffd@0: % Mu - mean vectors
wolffd@0: % Pi - priors
wolffd@0: %
wolffd@0: % lik - log likelihood of X 
wolffd@0: % likv - vector of log likelihoods
wolffd@0: % 
wolffd@0: % If 0 or 1 output arguments requested, lik is returned. If 2 output
wolffd@0: % arguments requested, [lik likv] is returned.
wolffd@0: 
wolffd@0: function [lik, likv]=mfa_cl(X,Lh,Ph,Mu,Pi);
wolffd@0: 
wolffd@0: N=length(X(:,1));
wolffd@0: D=length(X(1,:));
wolffd@0: K=length(Lh(1,:));
wolffd@0: M=length(Pi);
wolffd@0: 
wolffd@0: if (abs(sum(Pi)-1) > 1e-6) 
wolffd@0:   disp('ERROR: Pi should sum to 1');
wolffd@0:   return;
wolffd@0: elseif ((size(Lh) ~= [D*M K]) | (size(Ph) ~= [D 1]) | (size(Mu) ~= [M D]) ...
wolffd@0:   | (size(Pi) ~= [M 1] & size(Pi) ~= [1 M]))   
wolffd@0:   disp('ERROR in input matrix sizes');
wolffd@0:   return;
wolffd@0: end;  
wolffd@0: 
wolffd@0: tiny=exp(-744);
wolffd@0: const=(2*pi)^(-D/2);
wolffd@0: 
wolffd@0: I=eye(K);
wolffd@0: Phi=1./Ph;
wolffd@0: Phid=diag(Phi);
wolffd@0: for k=1:M  
wolffd@0:   Lht=Lh((k-1)*D+1:k*D,:);
wolffd@0:   LP=Phid*Lht;
wolffd@0:   MM=Phid-LP*inv(I+Lht'*LP)*LP';
wolffd@0:   dM=sqrt(det(MM));      	
wolffd@0:   Xk=(X-ones(N,1)*Mu(k,:)); 
wolffd@0:   XM=Xk*MM; 
wolffd@0:   H(:,k)=const*Pi(k)*dM*exp(-0.5*sum((XM.*Xk)'))'; 	
wolffd@0: end;
wolffd@0: 
wolffd@0: Hsum=rsum(H); 				
wolffd@0: 
wolffd@0: likv=log(Hsum+(Hsum==0)*tiny);
wolffd@0: lik=sum(likv);
wolffd@0: