ecm.m

一种在线聚类程序

%注意：用prestd进行归一化一定要将数据按行数小于列数的方式排列才是正确的
%也就是说，行代表的是变量个数，列代表样本数那才是正确的
clc
clear
load data
% [data,meanp,stdp] = prestd(data');
% data=data';
data=guiyihua(data);
mm=2;%隶属度上面的次数
[m n]=size(data);
Cc=[];%聚类中心
Ru=[];%聚类中心的半径
result=[];%记录每个样本归属哪个类
%利用第一个样本初始化第一个聚类中心点
Cc(1,:)=data(1,:);%n维数据
Ru(1)=0;%一维数据
result(1)=1;%第一个数据属于第一个类
%u=zeros(m,m);%初始化隶属度矩阵,最大化为聚类个数为样本个数
cluster=[];%初始化类别矩阵，用来存放各类的数据
Dthr=0.5;
%u(1,1)=1;%初始化第一个数据的隶属度
for i=2:m%总的样本循环
d=[];%存放每个样本到每个聚类中心的距离
[t,tl]=size(Cc);
%%%%%%%%%%%%%%%%%%%%%%%%%%Step1
for s1=1:t
    su=0;
    su=((data(i,:)-Cc(s1,:))*(data(i,:)-Cc(s1,:))');
     d(s1)=sqrt(su);
end
%%%%%%%%%%%计算uij
% for s1=1:t
%     if data(i,:)==Cc(s1,:)
%        u(i,:)=0;
%        u(i,s1)=1;
%     end 
% end
% if  1==isempty(find(1==u(i,:)))
%    sss=(d(i,:).^(-2/(mm-1)))*ones( t,1);
% for s1=1:t
%     u(i,s1)=(d(i,s1)^(-2/(mm-1)))/sss;
% end
% end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%step2
[dmin,sub]=min(d);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%step3,4,5
if dmin>=Ru(sub)            %%%%%%%%%%%%%%%
%取min(Sij)
    inter=[];
    for s1=1:t
   inter(s1)=(d(s1))+Ru(s1);
    end
    [S,nu]=min(inter);
    %增加新的聚类中心点
    if S>2*Dthr
        Cc(t+1,:)=data(i,:);
        Ru(t+1)=0;
        result(i)=t+1;%这个样本位新建类的第一个样本
    else
        Ru(nu)=S/2;
        if Ru(nu)<Dthr
        ls=(d(nu))/Ru(nu);
       %修改中心坐标
        %for f1=1:n
        Cc(nu,:)=data(i,:)-(data(i,:)-Cc(nu,:))/ls;
        %end
        result(i)=nu;
        else    %如果Ru(nu)>=Dthr，则创建新类
              Cc(t+1,:)=data(i,:);
              Ru(t+1)=0;
              result(i)=t+1;%这个样本为新建类的第一个样本
        end
    end
else
    result(i)=sub;
end
end
Cc
Ru
result

% %计算目标函数值
% z=0;
% [t,tl]=size(Cc);
% for i=1:m
%     for j=1:t
%         ts=0;
%         ts=(data(i,:)-Cc(j,:))*(data(i,:)-Cc(j,:))'/sqrt(n);
%         z=z+sqrt(ts);
%     end
% end
% z


%计算目标函数值
% z=0;
% [t,tl]=size(Cc);
% for i=1:m
%     for j=1:t
%         ts=0;
%         ts=(data(i,:)-Cc(j,:))*(data(i,:)-Cc(j,:))';
%         z=z+sqrt(ts);
%     end
% end


[k,n]=size(Cc);
for i=1:k
    ind=[]
    ind = find(result==i);%ind中记录的是所有属于第i类的样本的下标
    tt=length(ind);
    dist = sum((data(ind,:)-repmat(Cc(i,:),tt,1)).^2,2);
    %试验过了，当属于某类的元素个数为0时，执行上面一条指令的结果为0！！！
    sums(i)=sum(dist);
end
z=sum(sums);





dis=zeros(t,m);
for i=1:t
    record=[];
    record=find(i==result);
    tp=length(record);
    for j=1:tp
        dis(i,j)=sqrt((data(record(j),:)-Cc(i,:))*(data(record(j),:)-Cc(i,:))');
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%
dis=dis(:);
p=length(dis);
chaoguo=0;
for i=1:p
    if dis(i)>Dthr
        chaoguo=100000;
    end
end
chaoguo
Ru
z





figure(53)
zhongbiao=Cc;
for i=1:t
    x=zhongbiao(i,1);
    y=zhongbiao(i,2);
    plot(x,y,'ro');
    hold on;
end


maxres=result;
[p1,p2]=size(data);
wy=zeros(1,p2);%纯粹为了是画图语句方便写
data(p1+1,:)=wy;
x=[];
y=[];
for i=1:t
    ind=[];
    indt=repmat(p1+1,1,p1+1);
    ind=find(maxres==i);
    z=length(ind);
    indt(1:z)=ind;
    x(:,i)=data(indt,1);
    y(:,i)=data(indt,2);
end
x=x';
y=y';
hold on
plot(x(1,:),y(1,:),'g+');
plot(x(2,:),y(2,:),'bd');
plot(x(3,:),y(3,:),'ks');
plot(x(4,:),y(4,:),'mo');
plot(x(5,:),y(5,:),'bp');
plot(x(6,:),y(6,:),'gd');
plot(x(7,:),y(7,:),'m+');
plot(x(8,:),y(8,:),'b+');
plot(x(9,:),y(9,:),'kp');
plot(x(10,:),y(10,:),'mh');