data=input('请输入样本数据矩阵:');m=size(data,1);n=size(data,2);counter=0;k=input('请输入聚类数目:');whilekmdisp('您输入的聚类数目过大,请输入正确的k值');k=input('请输入聚类数目:');endifk==1disp('聚类数目不能为1,请输入正确的k值');k=input('请输入聚类数目:');end%产生k个零矩阵,M用来存放聚类中心M=cell(1,m);fori=1:kM{1,i}=zeros(1,n);endMold=cell(1,m);fori=1:kMold{1,i}=zeros(1,n);end%随机选取k个样本作为初始聚类中心%第一次聚类,使用初始聚类中心p=randperm(m);%产生m个不同的随机数fori=1:kM{1,i}=data(p(i),:);endwhiletruecounter=counter+1;disp('第');disp(counter);disp('次迭代');count=zeros(1,k);%初始化聚类CC=cell(1,k);fori=1:kC{1,i}=zeros(m,n);end%聚类fori=1:mgap=zeros(1,k);ford=1:kforj=1:ngap(d)=gap(d)+(M{1,d}(j)-data(i,j))^2;endend[y,l]=min(sqrt(gap));count(l)=count(l)+1;C{1,l}(count(l),:)=data(i,:);endMold=M;disp('聚类中心为:');fori=1:kdisp(M{1,i});enddisp('聚类结果为:');fori=1:kdisp(C{1,i});endsumvar=0;fori=1:kE=0;disp('单个误差平方和为:');forj=1:count(i)forh=1:nE=E+(M{1,i}(h)-C{1,i}(j,h))^2;endenddisp(E);sumvar=sumvar+E;enddisp('总体误差平方和为:');disp(sumvar);%计算新的聚类中心,更新M,并保存旧的聚类中心fori=1:kM{1,i}=sum(C{1,i})/count(i);end%检查前后两次聚类中心是否变化,若变化则继续迭代;否则算法停止;tally=0;fori=1:kifabs(Mold{1,i}-M{1,i})1e-5*ones(1,n)tally=tally+1;continue;elsebreak;endendiftally==kbreak;endend