杜丽萍的回答:本文重点是展示如何用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值为: 8647.63 后记 (1)pluskid指出K均值算法的初值对结果很重要,但是在运行时还没有发现类似的结果。也许Mathworks对该算法进行过优化。有时间会仔细研究下代码,将结果放上来。 转载仅供参考,版权属于原作者 |