somtoolbox2/som_probability_gmm.m
4dbef185
 function [pd,Pdm,pmd] = som_probability_gmm(D, sM, K, P)
 
 %SOM_PROBABILITY_GMM Probabilities based on a gaussian mixture model.
 %
 % [pd,Pdm,pmd] = som_probability_gmm(D, sM, K, P)
 % 
 %   [K,P] = som_estimate_gmm(sM,D);
 %   [pd,Pdm,pmd] = som_probability_gmm(D,sM,K,P);
 %   som_show(sM,'color',pmd(:,1),'color',Pdm(:,1))  
 %
 %  Input and output arguments:
 %   D    (matrix) size dlen x dim, the data for which the 
 %        (struct) data struct,     probabilities are calculated
 %   sM   (struct) map struct
 %        (matrix) size munits x dim, the kernel centers
 %   K    (matrix) size munits x dim, kernel width parameters
 %                 computed by SOM_ESTIMATE_GMM
 %   P    (matrix) size 1 x munits, a priori probabilities for each 
 %                 kernel computed by SOM_ESTIMATE_GMM
 %
 %   pd   (vector) size dlen x 1, probability of each data vector in 
 %                 terms of the whole gaussian mixture model
 %   Pdm  (matrix) size munits x dlen, probability of each vector in 
 %                 terms of each kernel
 %   pmd  (matrix) size munits x dlen, probability of each vector to 
 %                 have been generated by each kernel
 %
 % See also SOM_ESTIMATE_GMM.
 
 % Contributed to SOM Toolbox vs2, February 2nd, 2000 by Esa Alhoniemi
 % Copyright (c) by Esa Alhoniemi
 % http://www.cis.hut.fi/projects/somtoolbox/
 
 % ecco 180298 juuso 050100
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 % input arguments
 if isstruct(sM), M = sM.codebook; else M = sM; end
 [c dim] = size(M);
 
 if isstruct(D), D = D.data; end
 dlen = size(D,1);
 
 % reserve space for output variables
 pd = zeros(dlen,1); 
 if nargout>=2, Pdm = zeros(c,dlen); end
 if nargout==3, pmd = zeros(c,dlen); end
 
 % the parameters of each kernel
 cCoeff = cell(c,1);
 cCoinv = cell(c,1);
 for m=1:c, 
   co = diag(K(m,:));
   cCoinv{m} = inv(co);
   cCoeff{m} = 1 / ((2*pi)^(dim/2)*det(co)^.5);
 end
 
 % go through the vectors one by one
 for i=1:dlen, 
 
   x = D(i,:);
   
   % compute p(x|m)
   pxm = zeros(c,1); 
   for m = 1:c,
     dx     = M(m,:) - x;
     pxm(m) = cCoeff{m} * exp(-.5 * dx * cCoinv{m} * dx');
     %pxm(m) = normal(dx, zeros(1,dim), diag(K(m,:)));  
   end
   pxm(isnan(pxm(:))) = 0;
   
   % p(x|m)  
   if nargin>=2, Pdm(:,i) = pxm; end
   
   % P(x) = P(x|M) = sum( P(m) * p(x|m) )
   pd(i) = P*pxm; 
   
   % p(m|x) = p(x|m) * P(m) / P(x)
   if nargout==3, pmd(:,i) = (P' .* pxm) / pd(i); end
   
 end
 
 
 return; 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
 % subfunction normal
 %
 % computes probability of x when mean and covariance matrix
 % of a distribution are known
 
 function result = normal(x, mu, co)
 
 [l dim] = size(x);
 coinv   = inv(co);
 coeff   = 1 / ((2*pi)^(dim/2)*det(co)^.5);
 diff   = x - mu;
 result = coeff * exp(-.5 * diff * coinv * diff');