实验3Matlab聚类分析

Matlab提供了两种方法进行聚类分析。一种是利用clusterdata函数对样本数据进行一次聚类，其缺点为可供用户选择的面较窄，不能更改距离的计算方法；另一种是分步聚类：（1）找到数据集合中变量两两之间的相似性和非相似性，用pdist函数计算变量之间的距离；（2）用linkage函数定义变量之间的连接；（3）用cophenetic函数评价聚类信息；（4）用cluster函数创建聚类。1．Matlab中相关函数介绍1.1pdist函数调用格式：Y=pdist(X,’metric’)说明：用‘metric’指定的方法计算X数据矩阵中对象之间的距离。’X：一个m×n的矩阵，它是由m个对象组成的数据集，每个对象的大小为n。metric’取值如下：‘euclidean’：欧氏距离（默认）；‘seuclidean’：标准化欧氏距离；‘mahalanobis’：马氏距离；‘cityblock’：布洛克距离；‘minkowski’：明可夫斯基距离；‘cosine’：‘correlation’：‘hamming’：‘jaccard’：‘chebychev’：Chebychev距离。1.2squareform函数调用格式：Z=squareform(Y,..)说明：强制将距离矩阵从上三角形式转化为方阵形式，或从方阵形式转化为上三角形式。1.3linkage函数调用格式：Z=linkage(Y,’method’)说明：用‘method’参数指定的算法计算系统聚类树。Y：pdist函数返回的距离向量；method：可取值如下：‘single’：最短距离法（默认）；‘complete’：最长距离法；‘average’：未加权平均距离法；‘weighted’：加权平均法；‘centroid’：质心距离法；‘median’：加权质心距离法；‘ward’：内平方距离法（最小方差算法）返回：Z为一个包含聚类树信息的（m-1）×3的矩阵。1.4dendrogram函数调用格式：[H，T，…]=dendrogram(Z,p，…)说明：生成只有顶部p个节点的冰柱图（谱系图）。1.5cophenet函数调用格式：c=cophenetic(Z,Y)说明：利用pdist函数生成的Y和linkage函数生成的Z计算cophenet相关系数。1.6cluster函数调用格式：T=cluster(Z,…)说明：根据linkage函数的输出Z创建分类。1.7clusterdata函数调用格式：T=clusterdata(X,…)说明：根据数据创建分类。T=clusterdata(X,cutoff)与下面的一组命令等价：Y=pdist(X,’euclid’);Z=linkage(Y,’single’);T=cluster(Z,cutoff);2.Matlab程序2.1一次聚类法X=[1197812.593.531908;…;5750067.6238.015900];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.94698Step4创建聚类，并作出谱系图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',[.0020],'XTickLabel',[]);可以选择dendrogram显示的结点数目，这里选择20。结果显示可能可以分成三类重新调用K均值法改为分成三类[cidx3,cmeans3,sumd3,D3]=kmeans(X,3,'dist','sqEuclidean');P4=figure;clf;[silh3,h3]=silhouette(X,cidx3,'sqeuclidean');图上看，比前面的结果略有改善。将分类的结果展示出来P5=figure;clfptsymb={'bo','ro','go',',mo','c+'};MarkFace={[001],[.800],[0.50]};holdonfori=1:3clust=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});endholdoff运用高斯混合分布模型进行聚类分析分别用分布图、热能图和概率图展示结果等高线%等高线options=statset('Display','off');gm=gmdistribution.fit(X,3,'Options',options);P6=figure;clfscatter(X(:,1),X(:,2),10,'ro');holdonezcontour(@(x,y)pdf(gm,[x,y]),[-1515],[-1510]);holdoffP7=figure;clfscatter(X(:,1),X(:,2),10,'ro');holdonezsurf(@(x,y)pdf(gm,[x,y]),[-1515],[-1510]);holdoffview(33,24)热能图cluster1=(cidx3==1);cluster3=(cidx3==2);%通过观察，K均值方法的第二类是gm的第三类cluster2=(cidx3==3);%计算分类概率P=posterior(gm,X);P8=figure;clfplot3(X(cluster1,1),X(cluster1,2),P(cluster1,1),'r.')gridon;holdonplot3(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,'Component1PosteriorProbability')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({'Cluster1Score''Cluster2Score''Cluster3Score'},'location','NW');ylabel('ClusterMembershipScore');xlabel('PointRanking');通过AIC准则寻找最优的分类数高斯混合模型法的最大好处是给出分类好坏的标准AIC=zeros(1,4);NlogL=AIC;GM=cell(1,4);fork=1:4GM{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