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