wolffd@0: function [K,P] = som_estimate_gmm(sM, sD)
wolffd@0: 
wolffd@0: %SOM_ESTIMATE_GMM Estimate a gaussian mixture model based on map.
wolffd@0: %
wolffd@0: % [K,P] = som_estimate_gmm(sM, sD)
wolffd@0: %
wolffd@0: %  Input and output arguments:
wolffd@0: %   sM    (struct) map struct
wolffd@0: %   sD    (struct) data struct
wolffd@0: %         (matrix) size dlen x dim, the data to use when estimating
wolffd@0: %                  the gaussian kernels
wolffd@0: %
wolffd@0: %   K     (matrix) size munits x dim, kernel width parametes for 
wolffd@0: %                  each map unit
wolffd@0: %   P     (vector) size 1 x munits, a priori probability of each map unit
wolffd@0: %
wolffd@0: % See also SOM_PROBABILITY_GMM.
wolffd@0: 
wolffd@0: % Reference: Alhoniemi, E., Himberg, J., Vesanto, J.,
wolffd@0: %   "Probabilistic measures for responses of Self-Organizing Maps", 
wolffd@0: %   Proceedings of Computational Intelligence Methods and
wolffd@0: %   Applications (CIMA), 1999, Rochester, N.Y., USA, pp. 286-289.
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 250400
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: 
wolffd@0: [c, dim] = size(sM.codebook);
wolffd@0: M = sM.codebook;
wolffd@0: 
wolffd@0: if isstruct(sD), D = sD.data; else D = sD; end
wolffd@0: dlen = length(D(:,1));
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%
wolffd@0: % compute hits & bmus
wolffd@0: 
wolffd@0: [bmus, qerrs] = som_bmus(sM, D);
wolffd@0: hits = zeros(1,c);
wolffd@0: for i = 1:c, hits(i) = sum(bmus == i); end
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%
wolffd@0: % a priori  
wolffd@0: 
wolffd@0: % neighborhood kernel
wolffd@0: r  = sM.trainhist(end).radius_fin; % neighborhood radius
wolffd@0: if isempty(r) | isnan(r), r=1; end
wolffd@0: Ud = som_unit_dists(sM);
wolffd@0: Ud = Ud.^2; 
wolffd@0: r = r^2; 
wolffd@0: if r==0, r=eps; end % to get rid of div-by-zero errors
wolffd@0: switch sM.neigh, 
wolffd@0:  case 'bubble',   H = (Ud<=r); 
wolffd@0:  case 'gaussian', H = exp(-Ud/(2*r)); 
wolffd@0:  case 'cutgauss', H = exp(-Ud/(2*r)) .* (Ud<=r);
wolffd@0:  case 'ep',       H = (1-Ud/r) .* (Ud<=r);
wolffd@0: end  
wolffd@0: 
wolffd@0: % a priori prob. = hit histogram weighted by the neighborhood kernel
wolffd@0: P = hits*H;
wolffd@0: P = P/sum(P);              
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%
wolffd@0: % kernel widths (& centers)
wolffd@0: 
wolffd@0: K = ones(c, dim) * NaN; % kernel widths
wolffd@0: for m = 1:c,
wolffd@0:   w = H(bmus,m);
wolffd@0:   w = w/sum(w);
wolffd@0:   for i = 1:dim,
wolffd@0:     d = (D(:,i) - M(m,i)).^2;   % compute variance of ith
wolffd@0:     K(m,i) = w'*d;              % variable of centroid m
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%