实验3Matlab聚类分析

一种是利用 clusterdata函数对样本数据进行一次聚类，其缺点为可供用户选择的面较窄，不能更改距离的计算方法；

另一种是分步聚类：（1）找到数据集合中变量两两之间的相似性和非相似性，用pdist函数计算变量之间的距离；（2）用 linkage函数定义变量之间的连接；（3）用 cophenetic函数评价聚类信息；（4）用cluster函数创建聚类。

1．Matlab中相关函数介绍

1.1 pdist函数

调用格式：Y=pdist(X,’metric’)

说明：用 ‘metric’指定的方法计算 X 数据矩阵中对象之间的距离。’

X：一个m×n的矩阵，它是由m个对象组成的数据集，每个对象的大小为n。

metric’取值如下：

‘euclidean’：欧氏距离（默认）；‘seuclidean’：标准化欧氏距离；

‘mahalanobis’：马氏距离；‘cityblock’：布洛克距离；

‘minkowski’：明可夫斯基距离；‘cosine’：

‘correlation’： ‘hamming’：

‘jaccard’： ‘chebychev’：Chebychev距离。

1.2 squareform函数

调用格式：Z=squareform(Y,..)

说明： 强制将距离矩阵从上三角形式转化为方阵形式，或从方阵形式转化为上三角形式。

1.3 linkage函数

调用格式：Z=linkage(Y,’method’)

说 明：用‘method’参数指定的算法计算系统聚类树。

Y：pdist函数返回的距离向量；

method：可取值如下：

‘single’：最短距离法（默认）； ‘complete’：最长距离法；

‘average’：未加权平均距离法； ‘weighted’： 加权平均法；

‘centroid’：质心距离法； ‘median’：加权质心距离法；

‘ward’：内平方距离法（最小方差算法）

返回：Z为一个包含聚类树信息的（m-1）×3的矩阵。

1.4 dendrogram函数

调用格式：[H，T，…]=dendrogram(Z,p，…)

说明：生成只有顶部p个节点的冰柱图（谱系图）。

1.5 cophenet函数

调用格式：c=cophenetic(Z,Y)

说明：利用pdist函数生成的Y和linkage函数生成的Z计算cophenet相关系数。

1.6 cluster 函数

调用格式：T=cluster(Z,…)

说明：根据linkage函数的输出Z 创建分类。

1.7 clusterdata函数

调用格式：T=clusterdata(X,…)

说明：根据数据创建分类。

T=clusterdata(X,cutoff)与下面的一组命令等价：

Y=pdist(X,’euclid’);

Z=linkage(Y,’single’);

T=cluster(Z,cutoff);

2. Matlab程序

2.1 一次聚类法

X=[11978 12.5 93.5 31908;…;57500 67.6 238.0 15900];

T=clusterdata(X,0.9)

2.2 分步聚类

Step1 寻找变量之间的相似性

用pdist函数计算相似矩阵，有多种方法可以计算距离，进行计算之前最好先将数据用zscore函数进行标准化。

X2=zscore(X); %标准化数据

Y2=pdist(X2); %计算距离

Step2 定义变量之间的连接

Z2=linkage(Y2);

Step3 评价聚类信息

C2=cophenet(Z2,Y2); //0.94698

Step4 创建聚类，并作出谱系图

T=cluster(Z2,6);

H=dendrogram(Z2);

分类结果：{加拿大}，{中国，美国，澳大利亚}，{日本，印尼}，{巴西}，{前苏联}

用MATLAB做聚类分析

内容

• 展示如何使用MATLAB进行聚类分析 • 生成随机二维分布图形，三个中心 • K均值聚类 • 分层聚类

• 重新调用K均值法 • 将分类的结果展示出来

• 运用高斯混合分布模型进行聚类分析 • 通过AIC准则寻找最优的分类数

展示如何使用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 = figure;clf;

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

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

K均值聚类

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

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

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

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

分层聚类

eucD = pdist(X,'euclidean');

clustTreeEuc = linkage(eucD,'average'); cophenet(clustTreeEuc,eucD);

P3 = figure;clf;

[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 = figure;clf;

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

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

将分类的结果展示出来

P5 = figure;clf

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 = figure;clf

scatter(X(:,1),X(:,2),10,'ro'); hold on

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

P7 = figure;clf

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 = figure;clf

plot3(X(cluster1,1),X(cluster1,2),P(cluster1,1),'r.') grid on;hold 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 = figure;clf

[~,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值为： 87.63