EMfc.zip_EM聚类_K.

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(5,t,1); %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(t,5); %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(i,j)=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(Z,1); %just for the plots for i=1:t max(i)=Z(i,1); k=1; for j=2:5 if (Z(i,j)>max(i)) max(i)=Z(i,j); P(i,k)=0; k=j; else P(i,j)=0; end end P(i,k)=1; end n=sum(P,1); %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 done via the Z matrix each time). sumclusterX=X'*Z; diagnk=diag(nk); diag1nk=inv(diagnk); mean=sumclusterX*diag1nk; %Next one is the estimated intra-cluster VAriance-COvariance MAtrix of %the data. covma=cell(5,1); for j=1:5 Q1=X-ones(t,1)*mean(:,j)'; %x's-means Q2=Q1'; A=cell(t,1); for i=1:t A(i)={Q1(i,:)}; end B=cell(t,1); for i=1:t B(i)={Q2(:,i)}; end T=[Z(:,j) Z(:,j)]; C=cell(1,1); for e=1:t C{e}=(T(e,:)'.*B{e})*A{e}; end nom=[0 0;0 0]; for k = 1:length(C) nom=C{k}+nom; end %for covma{j}=nom./nk(j); end %for %PLOTS the 95% confidence ellipses & the (X,Y) pairs of each cluster in each iteration. % w=w+1; % subplot(2,2,w); % hold on, % % plot(X(:,1),X(:,2),'m.'); confreg(mean(:,1),covma{1},n(1),0.05); confreg(mean(:,2),covma{2},n(2),0.05); confreg(mean(:,3),covma{3},n(3),0.05); confreg(mean(:,4),covma{4},n(4),0.05); confreg(mean(:,5),covma{5},n(5),0.05); axis equal; zoom out; xlabel(g); title('Progress Graph'); ylabel('Data vs Clusters'); grid on; hold off pause(0.1); if w==4 %tests' whether your initial graph counter is 4 in order to reset him and subplot to w=0; %the correct position. else w=w; end %if %f(xi|theta-hat) matrix (n*k) -which is for every xi and every %mean,var-cov matrix. for k=1:5 for i=1:t Q11=X(i,:)-mean(:,k)'; %xi-mean Q22=Q11'; f(i,k)=(((2*pi)^-1)*(det(covma{k}))^(-1/2))*exp(-(1/2)*(Q11*(covma{k}^-1))*Q22); %The l(i,k) matrix is only computed in these two for's in order to avoid extra aggravation %of our program. l(i,k)=tk(k)*f(i,k); end %for end %for %the log-likelihood of our data given the iteration c is equal to: ll(c)=sum(log(sum(l,2))); % E-step (posterior expected value) for i=1:t for k=1:5 Z(i,k)=(tk(k)*f(i,k))/sum(tk'.*f(i,:)'); end end save sampleresults1 Z mean covma n P ll ; end %while % disp('Time to Find Clusters: '); % Q(p)=toc % LL{p}=ll; % MEANS{p}=mean; % COVS{p}=covma; % save MEANS&COVS100 MEANS COVS Q LL; % p %end %for p disp('Have a nice morning')