发布时间 : 星期三 文章粒子群聚类MATLAB代码更新完毕开始阅读d88c92e96294dd88d0d26b64
自己编写的粒子群(PSO)算法优化Kmeans聚类的MATLAB代码，MATLAB6.5/7.1测试通过，其它版本没测试。
粒子群算法采用惯性权重线性递减策略。
传统的Kmeans聚类依赖初始值的选择，容易陷入局部极值，也就是获得的聚类结果不一定是最优聚类。
本代码中的Kmeans聚类采用欧式距离的平方来衡量样本与类别中心的距离，目标函数是使所有样本到相应聚类中心的距离之和最小（聚类中心为每类中样本的平均值，有的还要求聚类中心之间的距离最大）。
聚类中心用粒子位置表示，每个粒子的位置包含K个聚类中心（K为聚类数目），通过调整聚类中心获得最优划分（最后的聚类中心依然是各类中样本的均值，利用粒子群算法给聚类中心加扰动，以增强跳出局部极值和寻找最优聚类的能力，如果多次扰动下聚类划分不变，则认为当前的聚类为最优聚类，也就是找到了最优聚类中心）。
输出结果解释如下：
fljg=[3 1 2 3 1 4 1]，表示1#、4#样本归为一类（类别编号为3），2#、5#、7#一类（类别编号为1），3#一类（类别编号为2），6#一类（类别编号为4），共有4类。类别编号只起类别区分作用，没有其它意义，也就是任意两个类别编号完全互换，聚类结果不变。
聚类数目K请根据需要修改，程序中为4。
每种聚类数目下的最优聚类请根据输出的适应度fg判断，适应度值越小越好，需多次运行判断。
参考了一些文献和代码，不再一一列出。
有关公式和框图如下：
公式中ij(代码中的ww)是个S×K的加权矩阵，矩阵元素非0即1，当样本i属于类j时，为1，不属于时为0。
初始种群生成代码做了改进，并用蓝色字体标出。 2009.11.30
代码如下：
%Kmeans Cluster Algorithm Based on Particle Optimization Algorithm % Author: Wang Yonglin (wylin77@126.com) clc;clear all; format long;
%------初始化------------求最小值 %数据，已经归一化
sam=[1.0000 1.0000 0.7476 0.6267 0.1696 0.0710 0.2532 0.8110
0.3188 0.3656 0.8707 0.7704 0.5559 0.5153 0.9213 0.7017
0.5548 0.7423 1.0000 0.5910 1.0000 1.0000 0.8976 1.0000
0.7800 0.7181 0.6875 1.0000 0.2115 0.0214 0.1573 0.8938
0.2680 0.3238 0.9036 0.8210 0.5874 0.3840 0.7037 0.5142
0.6928 0.6630 0.7368 0.8787 0.1818 0.0786 0.2295 0.3820
0.4256 0.4978 0.8429 0.9161 0.7133 0.3130 1.0000 0.7809]; N=50;%粒子数 c1=1.2;c2=1.2; wmax=0.9;wmin=0.4;
M=200;%代数
K=4;%类别数，根据需要修改%%%%%%%%%%%%%%%%%%%%%%%% [S D]=size(sam);%%样本数和特征维数 v=rand(N,K*D);%初始速度 %初始化分类矩阵 for i=1:N
clmat(i,:)=randperm(S);
clmat(i,clmat(i,:)>K)=ceil(rand(1,sum(clmat(i,:)>K))*K); end
%clmat(clmat>K)=fix(rand*K+1);%去掉本行 fitt=inf*ones(1,N);%初始化个体最优适应度 fg=inf;%初始化群体最优适应度 fljg=clmat(1,:);%当前最优分类 x=zeros(N,K*D);%初始化粒子群位置 y=x;%初始化个体最优解
pg=x(1,:);%初始化群体最优解 cen=zeros(K,D);%类别中心定维 fitt2=fitt;%粒子适应度定维
%------循环优化开始------------ for t=1:M for i=1:N
ww = zeros(S,K);% for ii = 1:S
ww(ii,clmat(i,ii)) = 1;%加权矩阵，元素非0即1 end
ccc=[];tmp=0; for j = 1:K
sumcs = sum(ww(:,j)*ones(1,D).*sam); countcs = sum(ww(:,j)); if countcs==0
cen(j,:) =zeros(1,D); else
cen(j,:) = sumcs/countcs; %求类别中心 end
ccc=[ccc,cen(j,:)];%串联聚类中心 aa=find(ww(:,j)==1); if length(aa)~=0
for k=1:length(aa)
tmp=tmp+(sum((sam(aa(k),:)-cen(j,:)).^2)); end end end
x(i,:)=ccc;
fitt2(i) = tmp; %Fitness value
end
%更新群体和个体最优解 for i=1:N
if fitt2(i) y(i,:)=x(i,:);%个体最优 if fitt2(i) pg=x(i,:);%群体最优 fg=fitt2(i);%群体最优适应度 fljg=clmat(i,:);%当前最优聚类 end end end bfit(t)=fg;%最优适应度记录 w = wmax - t*(wmax-wmin)/M;%更新权重 for i=1:N %更新粒子速度和位置 v(i,:)=w*v(i,:)+c1*rand(1,K*D).*(y(i,:)-x(i,:))+c2*rand(1,K*D).*(pg-x(i,:)); x(i,:)=x(i,:)+v(i,:); for k=1:K cen(k,:)=x((k-1)*D+1:k*D);%拆分粒子位置，获得K个中心 end %重新归类 for j=1:S tmp1=zeros(1,K); for k=1:K tmp1(k)=sum((sam(j,:)-cen(k,:)).^2);%每个样本关于各类的距离 end [tmp2 clmat(i,j)]=min(tmp1);%最近距离归类 end end end %------循环结束------------ fljg %最优聚类输出 fg %最优适应度输出 plot(bfit);%绘制最优适应度轨迹