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(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
% 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');