EM算法代码

function EMfc
%This is an implementation of the EM algorithm for Clustering via Gaussian
%mixture models using graphical on-line representation.

%Panagiotis Braimakis （s6990029）

%load simulated gaussian data.
clear;clc
%for p=1:1000
load data
%load M5    %uncomment only if you want to take results from k-means
            %algorithm

%if i load the M5 matrix then i get as starting values tha ones given via
%the kmeans algorithm which as he saw detects only spherical clusters
%（which does not always hold）
%So lets give initial random id‘s to the observations.

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%OPTIONAL SUB-SAMPLING%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5%%%%%%
 %sample from X
 t=10000;
%  Y=randperm（t）;
%  X=X（Y:）;
 
cidx=unidrnd（5t1）;
%in order to give a SEED for the Z matrix which classifies 
%each row of data to it‘s cluster we use the results of another clustering method  
tic;
Z=zeros（t5）;   %Setting up the storing matrix.
for i=1:t
for j=1:5
if （cidx（i）==j）     %cidx was the （t*1） id-matrix of each record （from 1 to 5）
Z（ij）=1;           %The new （t*5） Z matrix has only ones where needed
end %if
end %for
end %for

%Since we have the seed we can now implement the M-step of the algorithm &
%get the loop started between the “Maximization“ & “Expectation“ steps.

    w=0;    %for the graphs.
    c=2;    %starting counter.
    ll（1）=-88888888;  %starting values just to set the first difference between the log-likelihood‘s
                %in the first step （remember l（-1） & l（0） does not exist!!!）
    ll（2）=-77777777;
while ll（c）-ll（c-1） > 10^（-10）;
     
    c=c+1;  
    g=c-2   %iteration number.
    %all these are the estimates for the parameters of mixture models with
    %normal components.（Geoffrey J. Mclachlan &Kaye E. Basford）
    
    %nk row matrix is an estimate of the number of “points“ on each
    %cluster.
    nk=sum（Z1）;
    
    %just for the plots
    
    for i=1:t
max（i）=Z（i1）;
k=1;
for j=2:5
if  （Z（ij）>max（i））
max（i）=Z（ij）;
P（ik）=0;
k=j;
else
P（ij）=0;
end
end
P（ik）=1;
end
    
    n=sum（P1）;
        
        %Now the tk row-matrix estimates each time the probabilities of belonging
    %to the jth cluster （j=1:5）
    tk=nk./t;
    %%Now comes the mean‘s matrix if we suppose that i element of X belongs to the k-th
    %%cluster （that is