K均值聚类（Kmeans）

Sigma = [1, 0; 0, 1];
mu1 = [1, -1];
x1 = mvnrnd(mu1, Sigma, 200);
mu2 = [5.5, -4.5];
x2 = mvnrnd(mu2, Sigma, 200);
mu3 = [1, 4];
x3 = mvnrnd(mu3, Sigma, 200);
mu4 = [6, 4.5];
x4 = mvnrnd(mu4, Sigma, 200);
mu5 = [9, 0.0];
x5 = mvnrnd(mu5, Sigma, 200);
% obtain the 1000 data points to be clustered
X = [x1; x2; x3; x4; x5];
% Show the data point
plot(x1(:,1), x1(:,2), 'r.'); hold on;
plot(x2(:,1), x2(:,2), 'b.');
plot(x3(:,1), x3(:,2), 'k.');
plot(x4(:,1), x4(:,2), 'g.');
plot(x5(:,1), x5(:,2), 'm.');

save myX  %将X存储到文件中，在其他文件中load就可以了

结果如下：

 1 % 初始聚类中心
 2 mu0=[X(1,:); X(205,:);X(405,:);X(605,:);X(805,:)];
 3 mu1=zeros(5,c);
 4 lable=zeros(r,1);
 5 % first cluster
 6     dist=zeros(5,1);
 7     for k=1:r
 8         for n=1:5
 9             dist(n)=norm(X(k,:)-mu0(n,:));
10         end
11         lable(k)=find(dist==min(dist));
12     end    
13 
14 %     X1=[lable X];
15     sum=zeros(5,2);
16     count=zeros(5,1);
17     
18  % 第一次聚类
19     for n=1:5    
20         for k=1:r
21             if lable(k)==n 
22              sum(n,:)=sum(n,:)+X(k,:);
23              count(n)=count(n)+1;
24             end   
25         end   
26         mu1(n,:)=sum(n,:)/count(n);
27     end
28 
29  % square error  
30     e=zeros(5,1);
31     esum=0;
32     for n=1:5
33         for k=1:r
34             if lable(k)==n
35                  e(n,:)=e(n,:)+norm(X(k,:)-mu1(n,:));
36             end
37         end
38         esum=esum+e(n,:);
39     end
40 
41 % ***************** start Iteration **************
42 
43 countNotX=0;% 到达一定值后，停止迭代
44 s=0;
45 count1=zeros(5,1);
46 while true
47     s=s+1;
48 %random choose X to update  
49     ran=round(1+999*rand(1));%总共1000个样本
50     Xrand=X(ran,:);
51     rou=zeros(5,1);
52     rouMin=inf;
53     newClassLable=0;
54     if count(lable(ran))~=0
55         for n=1:5
56             if lable(ran)==n
57                 rou(n)=norm(Xrand-mu1(n,:))*count(n)/(count(n)-1);
58             else
59                 rou(n)=norm(Xrand-mu1(n,:))*count(n)/(count(n)+1);
60             end
61             if rou(n)<=rouMin;
62                 rouMin=rou(n);
63                 newClassLable=n;
64             end
65         end
66        if rouMin<rou(lable(ran))          
67             countNotX=0;
68 %              new mu  这里用课件中的公式，不知道哪里错了，总是不行，就用最蠢的办法了
69            for n=1:5    
70                 for k=1:r
71                     if lable(k)==n 
72                      sum(n,:)=sum(n,:)+X(k,:);
73                      count1(n)=count1(n)+1;
74                     end   
75                 end   
76                 mu1(n,:)=sum(n,:)/count1(n);
77            end
78 %             new count
79             count(lable(ran))=count(lable(ran))-1;
80             count(newClassLable)=count(newClassLable)+1;
81 %           new esum
82             e(lable(ran),:)=e(lable(ran),:)-rou(lable(ran));
83             e(newClassLable,:)=e(newClassLable,:)+rouMin;
84              esum=esum+rouMin-rou(lable(ran));  
85              lable(ran)=newClassLable ;
86          else
87 %            disp('No')
88            countNotX= countNotX+1;
89            if countNotX==600
90                %when J is not changed for 600 times,reckon OK
91                break
92            end
93        end
94     end
95 end