如何编写求K-均值聚类算法的Matlab程序?

在聚类分析中，K-均值聚类算法(k-means

algorithm)是无监督分类中的一种基本方法，其也称为C-均值算法，其基本思想是:通过迭代的方法，逐次更新各聚类中心的值，直至得到最好的聚类结果。

假设要把样本集分为c个类别，算法如下:

(1)适当选择c个类的初始中心

(2)在第k次迭代中，对任意一个样本，求其到c个中心的距离，将该样本归到距离最短的中心所在的类，

(3)利用均值等方法更新该类的中心值

(4)对于所有的c个聚类中心，如果利用(2)(3)的迭代法更新后，值保持不变，则迭代结束，否则继续迭代。

下面介绍作者编写的一个分两类的程序，可以把其作为函数调用。

%%

function

[samp1,samp2]=kmeans(samp)

作为调用函数时去掉注释符

samp=[11.1506

6.7222

2.3139

5.9018

11.0827

5.7459

13.2174

13.8243

4.8005

0.9370

12.3576]

%样本集

[l0

l]=size(samp)

%%利用均值把样本分为两类，再将每类的均值作为聚类中心

th0=mean(samp)n1=0n2=0c1=0.0c1=double(c1)c2=c1for

i=1:lif

samp(i)<th0

c1=c1+samp(i)n1=n1+1elsec2=c2+samp(i)n2=n2+1endendc1=c1/n1c2=c2/n2

%初始聚类中心t=0cl1=c1cl2=c2

c11=c1c22=c2

%聚类中心while

t==0samp1=zeros(1,l)

samp2=samp1n1=1n2=1for

i=1:lif

abs(samp(i)-c11)<abs(samp(i)-c22)

samp1(n1)=samp(i)

cl1=cl1+samp(i)n1=n1+1

c11=cl1/n1elsesamp2(n2)=samp(i)

cl2=cl2+samp(i)n2=n2+1

c22=cl2/n2endendif

c11==c1

&&

c22==c2t=1endcl1=c11cl2=c22

c1=c11c2=c22

end

%samp1,samp2为聚类的结果。

初始中心值这里采用均值的办法，也可以根据问题的性质，用经验的方法来确定，或者将样本集随机分成c类，计算每类的均值。

k-均值算法需要事先知道分类的数量，这是其不足之处。

展示如何使用MATLAB进行聚类分析

分别运用分层聚类、K均值聚类以及高斯混合模型来进行分析，然后比较三者的结果

生成随机二维分布图形，三个中心

% 使用高斯分布（正态分布）

% 随机生成3个中心以及标准差

s = rng(5,'v5normal')

mu = round((rand(3,2)-0.5)*19)+1

sigma = round(rand(3,2)*40)/10+1

X = [mvnrnd(mu(1,:),sigma(1,:),200)...

mvnrnd(mu(2,:),sigma(2,:),300)...

mvnrnd(mu(3,:),sigma(3,:),400)]

% 作图

P1 = figureclf

scatter(X(:,1),X(:,2),10,'ro')

title('研究样本散点分布图')

K均值聚类

% 距离用传统欧式距离，分成两类

[cidx2,cmeans2,sumd2,D2] = kmeans(X,2,'dist','sqEuclidean')

P2 = figureclf

[silh2,h2] = silhouette(X,cidx2,'sqeuclidean')

从轮廓图上面看，第二类结果比较好，但是第一类有部分数据表现不佳。有相当部分的点落在0.8以下。

分层聚类

eucD = pdist(X,'euclidean')

clustTreeEuc = linkage(eucD,'average')

cophenet(clustTreeEuc,eucD)

P3 = figureclf

[h,nodes] = dendrogram(clustTreeEuc,20)

set(gca,'TickDir','out','TickLength',[.002 0],'XTickLabel',[])

可以选择dendrogram显示的结点数目，这里选择20 。结果显示可能可以分成三类

重新调用K均值法

改为分成三类

[cidx3,cmeans3,sumd3,D3] = kmeans(X,3,'dist','sqEuclidean')

P4 = figureclf

[silh3,h3] = silhouette(X,cidx3,'sqeuclidean')

图上看，比前面的结果略有改善。

将分类的结果展示出来

P5 = figureclf

ptsymb = {'bo','ro','go',',mo','c+'}

MarkFace = {[0 0 1],[.8 0 0],[0 .5 0]}

hold on

for i =1:3

clust = find(cidx3 == i)

plot(X(clust,1),X(clust,2),ptsymb{i},'MarkerSize',3,'MarkerFace',MarkFace{i},'MarkerEdgeColor','black')

plot(cmeans3(i,1),cmeans3(i,2),ptsymb{i},'MarkerSize',10,'MarkerFace',MarkFace{i})

end

hold off

运用高斯混合分布模型进行聚类分析

分别用分布图、热能图和概率图展示结果 等高线

% 等高线

options = statset('Display','off')

gm = gmdistribution.fit(X,3,'Options',options)

P6 = figureclf

scatter(X(:,1),X(:,2),10,'ro')

hold on

ezcontour(@(x,y) pdf(gm,[x,y]),[-15 15],[-15 10])

hold off

P7 = figureclf

scatter(X(:,1),X(:,2),10,'ro')

hold on

ezsurf(@(x,y) pdf(gm,[x,y]),[-15 15],[-15 10])

hold off

view(33,24)

热能图

cluster1 = (cidx3 == 1)

cluster3 = (cidx3 == 2)

% 通过观察，K均值方法的第二类是gm的第三类

cluster2 = (cidx3 == 3)

% 计算分类概率

P = posterior(gm,X)

P8 = figureclf

plot3(X(cluster1,1),X(cluster1,2),P(cluster1,1),'r.')

grid onhold on

plot3(X(cluster2,1),X(cluster2,2),P(cluster2,2),'bo')

plot3(X(cluster3,1),X(cluster3,2),P(cluster3,3),'g*')

legend('第 1 类','第 2 类','第 3 类','Location','NW')

clrmap = jet(80)colormap(clrmap(9:72,:))

ylabel(colorbar,'Component 1 Posterior Probability')

view(-45,20)

% 第三类点部分概率值较低，可能需要其他数据来进行分析。

% 概率图

P9 = figureclf

[~,order] = sort(P(:,1))

plot(1:size(X,1),P(order,1),'r-',1:size(X,1),P(order,2),'b-',1:size(X,1),P(order,3),'y-')

legend({'Cluster 1 Score' 'Cluster 2 Score' 'Cluster 3 Score'},'location','NW')

ylabel('Cluster Membership Score')

xlabel('Point Ranking')

通过AIC准则寻找最优的分类数

高斯混合模型法的最大好处是给出分类好坏的标准

AIC = zeros(1,4)

NlogL = AIC

GM = cell(1,4)

for k = 1:4

GM{k} = gmdistribution.fit(X,k)

AIC(k)= GM{k}.AIC

NlogL(k) = GM{k}.NlogL

end

[minAIC,numComponents] = min(AIC)

按AIC准则给出的最优分类数为： 3 对应的AIC值为： 8647.63

后记

（1）pluskid指出K均值算法的初值对结果很重要，但是在运行时还没有发现类似的结果。也许Mathworks对该算法进行过优化。有时间会仔细研究下代码，将结果放上来。

分享：

56

喜欢

4

赠金笔

阅读(21209)┊ 评论 (4)┊ 收藏(1) ┊转载原文 ┊ 喜欢▼ ┊打印┊举报

前一篇：[转载]拉普拉斯矩阵

后一篇：[转载]用matlab做聚类分析

%%k均值聚类的示例代码：

X = [randn(100,2)+ones(100,2)...

randn(100,2)-ones(100,2)]

opts = statset('Display','final')

[idx,ctrs] = kmeans(X,2,...

'Distance','city',...

'Replicates',5,...

'Options',opts)

%5 iterations, total sum of distances = 284.671

%4 iterations, total sum of distances = 284.671

%4 iterations, total sum of distances = 284.671

%3 iterations, total sum of distances = 284.671

%3 iterations, total sum of distances = 284.671

plot(X(idx==1,1),X(idx==1,2),'r.','MarkerSize',12)

hold on

plot(X(idx==2,1),X(idx==2,2),'b.','MarkerSize',12)

plot(ctrs(:,1),ctrs(:,2),'kx',...

'MarkerSize',12,'LineWidth',2)

plot(ctrs(:,1),ctrs(:,2),'ko',...

'MarkerSize',12,'LineWidth',2)

legend('Cluster 1','Cluster 2','Centroids',...

'Location','NW')

你可以help下cluster,matlab里面还自带很多这种例子

---