em_initialize.m

EM分群,matlab程式碼,用來分群用的

function [init_par, init_weights, clust_pars, clust_flag] = EM_Initialize(data, module_type, fig_hand)

clust_flag = 1;

n_mix = length(module_type);

clust_pars = zeros(n_mix,2);

init_weights = rand(1, n_mix);
init_weights = init_weights/sum(init_weights);

init_par    = cell(n_mix,1);
n_points    = zeros(n_mix,1);

mu_init     = zeros(n_mix,1);
sig_init    = zeros(n_mix,1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
method = 'k-means';

if ~exist('kmeans')
    method = 'hand';
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(fig_hand)
hold on;


if n_mix > 1
    
    switch method
        
        case 'k-means'
            
            [IDX,mu_init,sumd]  = kmeans(data, n_mix,'display','off','emptyaction','singleton');
         
            plot(mu_init,0,'rx');              
            
            ii = 1;
            while ii<=n_mix & clust_flag>0
                n_points(ii) = sum(IDX==ii);
                if n_points(ii)>1 & sumd(ii)>0
                    sig_init(ii) = sqrt(sumd(ii)./(n_points(ii)-1));
                else
                    clust_flag = 0;
                end
                ii = ii + 1;
            end            

         
            
        case 'hand' % could be removed
            
            hold on
            
            stringa = strcat('Use mouse to indicate initial guess for n=',num2str(n_mix),' centers');
            title(stringa);        
            
            for ii=1:n_mix
                
                k = waitforbuttonpress;
                point1 = get(gca,'CurrentPoint');    % button down detected
                finalRect = rbbox;                   % return figure units
                point2 = get(gca,'CurrentPoint');    % button up detected
                point1 = point1(1,1:2);              % extract x and y
                point2 = point2(1,1:2);
                p1 = min(point1,point2);             % calculate locations
                offset = abs(point1-point2);         % and dimensions
                x = [p1(1) p1(1)+offset(1) p1(1)+offset(1) p1(1) p1(1)];
                y = [p1(2) p1(2) p1(2)+offset(2) p1(2)+offset(2) p1(2)];
                
                mu_init(ii) = 0.5*(min(x) + max(x));	
                subdata = data(abs(data-mu_init(ii))<offset(1));		
                sig_init(ii) = std(subdata);
                plot(mu_init(ii),0,'rx');
                drawnow;
                
            end
            
    end
    
else
    
    mu_init = mean(data);
    sig_init = std(data);    
    plot(mu_init,0,'rx');   
    
end

hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if clust_flag>0
       
    clust_pars = sortrows([mu_init sig_init]);
    
    mu_init  = clust_pars(:,1);
    sig_init = clust_pars(:,2);
    
    
    for ii=1:n_mix
        init_par{ii} = EM_Invert_pars(mu_init(ii), sig_init(ii), module_type{ii});
    end             

end