www.pudn.com > gmm_utilities.zip > gmm_em_auto.m

function [g, w] = gmm_em_auto(s, K, N, type) 
%function [g, w] = gmm_em_auto(s, K, N) 
% 
% INPUTS: 
%   s - samples 
%   K - number of components in resultant gmm 
%   N - number of iterations of EM 
% 
% OUTPUT: 
%   g - resultant gmm 
%   w - negative log-likelihood of fit 
% 
% Tim Bailey 2005. 
if nargin == 3, type = 1; end 
 
NA = 10; % number of total search iterations 
NK = 5; % number of kmeans iterations 
NE = 5; % number of trial EM iterations 
 
% Auto search for initial gmm components 
switch type 
case 1 % kmeans generates entire gmm 
    Ps = repmat(eye(size(s,1))*eps, [1,1,K]); 
    wb = 1e99; 
    for i=1:NA 
        g.x = s(:, ceil(stratified_random(K)*N)); 
        [g.x,g.P,g.w] = kmeans(s, g.x, NK); 
        g.P = g.P + Ps; % protect against rank deficiency 
        [gem, w] = gmm_em(s, g, NE); 
        if w < wb 
            gbest = gem; 
            wb = w; 
        end 
    end 
     
case 2 % kmeans generates means only, P and w are computed from ensemble statistics 
    g.w = ones(1,K)/K; 
    [xe,Pe] = sample_mean(s); 
    g.P = repmat(Pe/K, [1,1,K]); 
     
    wb = 1e99; 
    for i=1:NA 
        g.x = s(:, ceil(stratified_random(K)*N)); 
        [g.x,dmy,dmy] = kmeans(s, g.x, NK); 
        [gem, w] = gmm_em(s, g, NE); 
        if w < wb 
            gbest = gem; 
            wb = w; 
        end 
    end     
     
otherwise 
    error('Invalid initialisation type') 
end 
 
% Optimise best gmm 
[g, w] = gmm_em(s, gbest, N);