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