K-means与ISODATA算法的比较研究的程序

1.K-means算法程序% k-均聚类算法clcclear;% main variablesdim = 2; % 模式样本维数k = 8; % 设有k个聚类中心load('file.txt');PM=file;% 模式样本矩阵N = size(PM,1);figure();subplot(1,2,1);for(i=1:N)plot(PM(i,1),PM(i,2), '*r'); % 绘出原始的数据点hold onendxlabel('X');ylabel('Y');title('聚类之前的数据点');CC = zeros(k,dim); % 聚类中心矩阵，CC(i,:)初始值为i号样本向量D = zeros(N,k); % D(i,j)是样本i和聚类中心j的距离C = cell(1,k); %% 聚类矩阵，对应聚类包含的样本。

初始状况下，聚类i(i<k)的样本集合为[i],聚类k的样本集合为[k,k+1,...N]for i = 1:k-1C{i} = [i];endC{k} = k:N;B = 1:N; % 上次迭代中，样本属于哪一聚类，设初值为1B(k:N) = k;for i = 1:kCC(i,:) = PM(i,:);endwhile 1change = 0;%用来标记分类结果是否变化% 对每一个样本i，计算到k个聚类中心的距离for i = 1:Nfor j = 1:k% D(i,j) = eulerDis( PM(i,:), CC(j,:) );D(i,j) = sqrt((PM(i,1) - CC(j,1))^2 + (PM(i,2) - CC(j,2))^2); endt = find( D(i,:) == min(D(i,:)) ); % i属于第t类if B(i) ~= t % 上次迭代i不属于第t类change = 1;% 将i从第B(i)类中去掉t1 = C{B(i)};t2 = find( t1==i );t1(t2) = t1(1);t1 = t1(2:length(t1));C{B(i)} = t1;C{t} = [C{t},i]; % 将i加入第t类B(i) = t;endendif change == 0 %分类结果无变化，则迭代停止break;end% 重新计算聚类中心矩阵CCfor i = 1:kCC(i,:) = 0;iclu = C{i};for j = 1:length(iclu)CC(i,:) = PM( iclu(j),: )+CC(i,:);endCC(i,:) = CC(i,:)/length(iclu);endendsubplot(1,2,2);plot(CC(:,1),CC(:,2),'o')hold onfor(i=1:N)if(B(1,i)==1)plot(PM(i,1),PM(i,2),'*b'); %作出第一类点的图形hold onelseif(B(1,i)==2)plot(PM(i,1),PM(i,2), '*r'); %作出第二类点的图形hold onelseif(B(1,i)==3)plot(PM(i,1),PM(i,2),'*g'); %作出第三类点的图形hold onelseif(B(1,i)==4)plot(PM(i,1),PM(i,2),'*m'); %作出第4类点的图形hold onelseif(B(1,i)==5)plot(PM(i,1),PM(i,2),'*c'); %作出第5类点的图形hold onelseif(B(1,i)==6)plot(PM(i,1),PM(i,2),'*b'); %作出第6类点的图形hold onelseif(B(1,i)==7)plot(PM(i,1),PM(i,2),'*y'); %作出第7类点的图形hold onelseplot(PM(i,1),PM(i,2), '*k'); %作出第四类点的图形hold onendendxlabel('X');ylabel('Y');title('聚类之后的数据点');% 打印C,CCfor i = 1:k %输出每一类的样本点标号str=['第' num2str(i) '类包含点: ' num2str(C{i})];disp(str);end;2.ISODATA算法程序Mainclc;clear all;fprintf('ISODATA分类算法\n');load('file.txt');x=file;plot(x(:,1),x(:,2),'.');%显示待分类点的坐标位置title('待分类的样本点坐标位置');[N,n]=size(x);%获得样本数目和维数，N为样本数目，n为样本维度K=4;%预期的聚类中心数目thet_N=1;%每一聚类域中最少的样本数目thet_S=1;%一个聚类域中样本距离分布的标准差thet_C=1;%两个聚类中心间的最小距离L=2;%在一次迭代运算中可以合并的聚类中心的最多对数I=200;%迭代运算的次数Nz=1;%决定初始聚类中心的个数ifNz==1z0=[1];z=[x(z0(1),:)];%初始化聚类中心fprintf('选样本点x(%d)作为初始聚类中心\n',z0(1));endifNz==3z0=[1,2,7];z=[x(z0(1),:);x(z0(2),:);x(z0(3),:)];%初始化聚类中心fprintf('选样本点x(%d)、x(%d)、x(%d)作为初始聚类中心\n',z0(1),z0(2),z0(3));endifNz==4z0=[1,2,7,8];z=[x(z0(1),:);x(z0(2),:);x(z0(3),:);x(z0(4),:)];%初始化聚类中心fprintf('选样本点x(%d)、x(%d)、x(%d)、x(%d)作为初始聚类中心\n',z0(1),z0(2),z0(3),z0(4)); endlabel=zeros(1,N);for k=1:I[z,mean_D1,mean_D,delta,Nc,Ni,label]=get_einf(x,z,label);for j=1:Ncif Ni(j)<thet_N%如果类中样本数目太少，则取消该样本子集for i=1:N;if label(i)==jfor m=i:Nx(m,:)=x(m+1,:);endendendend%[z,mean_D1,mean_D,delta,Nc,Ni,label]=get_einf(x,z,label);%重新进行聚类并计算类的基本信息if k==Ithet_C=0;[z] = merge(z,Ni,L,thet_C);continue;endif Nc<=(K/2)%如果聚类中心的数目小于或等于规定值的一半，则对已有聚类进行分裂处理[z] = split(K,z,Ni,thet_N,delta,thet_S,mean_D1,mean_D);continue;endif (Nc>=2*K)||(mod(k,2)==0)%如果Nc>=2K，或迭代运算的次数是偶数次，进行合并处理[z] = merge(z,Ni,L,thet_C);continue;else%否则进行分裂处理[z] = split(K,z,Ni,thet_N,delta,thet_S,mean_D1,mean_D);continue;endend[z,mean_D1,mean_D,delta,Nc,Ni,label]=get_einf(x,z,label);fprintf('样本被分为%d类\n',Nc);for j=1:Ncfprintf('第%d类的聚类中心坐标为:(%f,%f)\n',j,z(j,1),z(j,2));%显示聚类中心坐标endfor i=1:Nfprintf('第%d个样本的坐标:(%d,%d),属于第%d类\n',i,x(i,1),x(i,2),label(i));%显示样本坐标和属于第几类if label(i)==1plot(x(i,1),x(i,2),'.B');hold on;endif label(i)==2plot(x(i,1),x(i,2),'.R');hold on;endif label(i)==3plot(x(i,1),x(i,2),'.g');hold on;if label(i)==4plot(x(i,1),x(i,2),'.m');hold on;endif label(i)==5plot(x(i,1),x(i,2),'.c');hold on;endif label(i)==6plot(x(i,1),x(i,2),'.b');hold on;endif label(i)==7plot(x(i,1),x(i,2),'.y');hold on;endif label(i)==8plot(x(i,1),x(i,2),'.k');hold on;endtitle('分类结果图');end%==================================================================== %函数功能：基于聚类中心进行重新分类，并获得新的分类的基本信息%输入参数：x：样本% z：原聚类中心% label：原样本类别标签%输出参数：z：新的聚类中心% mean_D1：每个类中样本到聚类中心的平均距离% mean_D：全部模式样本和其对应聚类中心的总平均距离% delta：每个聚类中样本距离的标准差向量% Nc：分类的数目% Ni：每个类中样本的数目% label：新的类别标签%==================================================================== function [ z,mean_D1,mean_D,delta,Nc,Ni,label ] = get_einf( x,z,label )[N,n]=size(x);[Nc,n]=size(z);d=zeros(1,Nc);Ni=zeros(1,Nc);%用于保存每一类中样本的个数mean_D1=zeros(1,Nc);%用于保存每个类中样本到聚类中心的平均距离mean_D=0;%用于保存全部模式样本和其对应聚类中心的总平均距离delta=zeros(Nc,n);%用于保存每个聚类中样本距离的标准差向量for i=1:Nfor j=1:Ncd(j)=norm(x(i,:)-z(j,:));endd_min=d(1);for j=1:Ncif d(j)<=d_mind_min=d(j);label(i)=j;%对样本进行分类endendendfor j=1:Ncz(j,:)=zeros(1,n);for i=1:Nif label(i)==jNi(j)=Ni(j)+1;z(j,:)=z(j,:)+x(i,:);endendz(j,:)=z(j,:)/Ni(j);%计算新的聚类中心endfor j=1:Ncfor i=1:Nif label(i)==jmean_D1(j)=mean_D1(j)+norm(x(i,:)-z(j,:));for k=1:ndelta(j,k)=delta(j,k)+(x(i,k)-z(j,k))^2;endendendmean_D=mean_D+mean_D1(j);mean_D1(j)=mean_D1(j)/Ni(j);%获得每个类中样本到聚类中心的平均距离delta(j,:)=sqrt(delta(j,:)/Ni(j));%获得每个聚类中样本距离的标准差向量endmean_D=mean_D/N;%获得全部模式样本和其对应聚类中心的总平均距离end%====================================================================%函数功能：对满足条件的类进行合并%输入参数：z：原聚类中心% Ni：每个类中样本的数目% L：在一次迭代运算中可以合并的聚类中心的最多对数% thet_C：两个聚类中心间的最小距离%输出参数：z_new：新的聚类中心%====================================================================function [z_new] = merge(z,Ni,L,thet_C)[Nc,n]=size(z);D=zeros(Nc,Nc);%记录聚类中心间两两之间的距离merge=zeros(L,2);%用于保存L个最小距离的聚类中心的索引D_min=zeros(1,L);%用于保存L个最小距离D_max=[0,1,1];%保存最大距离的相关信息%第一个参数保存最大距离%第二、三个参数保存最大距离的聚类中心的索引for j=1:(Nc-1)%找到最大距离及其聚类中心的索引，用于初始化D_min和mergefor l=(j+1):NcD(j,l)=norm(z(j,:)-z(l,:));if D(j,l)>=D_maxD_max(1)=D(j,l);D_max(2)=j;D_max(3)=l;endendendfor i=1:Lmerge(i,:)=[D_max(2),D_max(3)];%初始化D_min和mergeD_min(i)=D_max(1);end%------------%以下部分的目的是找到L个最小距离及其聚类中心的索引max=[0,0];%用于保存L个最小距离中的最大距离，及其在D_min和merge中的索引for j=1:(Nc-1)for l=(j+1):Ncfor i=1:Lif max<=D_min(i)max(1)=D_min(i);max(2)=i;endendif D(j,l)<=max(1)%如果新的距离小于L个最小距离中的最大距离，用新的距离替代其中的原最大距离merge(max(2),:)=[j,l];D_min(max(2))=D(j,l);endendend%------------k=0;for i=1:Lif D_min(i)<thet_C%如果距离在最小的L个距离中，且小于thet_C，则满足合并条件，进行合并k=k+1;%以下部分将待合并的两个聚类中心进行赋值（在此处值相等）z(merge(i,1),:)=(Ni(merge(i,1))*z(merge(i,1))+Ni(merge(i,2))*z(merge(i,2)))/(Ni(merge(i,1))+Ni(me rge(i,2)));z(merge(i,2),:)=z(merge(i,1),:);endendz_new=zeros(Nc-k,n);k1=0;for j=1:Nc-1for i=j+1:Ncif z(j,:)==z(i,:)%如果出现相等的聚类中心，则去掉后面一个。