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