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 argumentsif isstruct(sM), M = sM.codebook; else M = sM; end[c dim] = size(M);if isstruct(D), D = D.data; enddlen = size(D,1);% reserve space for output variablespd = zeros(dlen,1); if nargout>=2, Pdm = zeros(c,dlen); endif nargout==3, pmd = zeros(c,dlen); end% the parameters of each kernelcCoeff = 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 onefor 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  endreturn; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subfunction normal%% computes probability of x when mean and covariance matrix% of a distribution are knownfunction 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');