efcm.m

在ECM的基础上加入了模糊关系

%%%%注意，对于聚类不许调整Dthr出现聚类数结果为1的情况，
%当出现聚类数为1时其目标函数值是假的，不是真正的样本聚类结果所具有的目标函数值
%注意：用prestd进行归一化一定要将数据按行数小于列数的方式排列才是正确的
%也就是说，行代表的是变量个数，列代表样本数那才是正确的

clc
clear
load data1
data=data1;
[m n]=size(data);
data(:,4:6)=[];
iitr=[1:4:m 2:4:m 3:4:m];
iitst=[4:4:m];
dddd=[];
dddd=data(iitr,:);
data=[];
data=dddd;
[m n]=size(data);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%以上部分是为了得到样本类标签而改动过得
% data(:,4)=[];
% [data,meanp,stdp] = prestd(data');
% data=data';
data=guiyihua(data);
symbol=1;%初始化标志位
mm=3;%隶属度上面的次数

Cc=[];%聚类中心
Ru=[];%聚类中心的半径
result=[];%记录每个样本归属哪个类
%利用第一个样本初始化第一个聚类中心点
Cc(1,:)=data(1,:);%n维数据
Ru(1)=0;%一维数据
result(1)=1;%第一个数据属于第一个类
u=zeros(m,10);%初始化隶属度矩阵,最大化为聚类个数为样本个数
cluster=[];%初始化类别矩阵，用来存放各类的数据
Dthr=0.5;
% Dthr=0.045; %对应310数据，误差0.0401
% Dthr=0.5;  %对应703数据，误差0.0074
% Dthr=0.5;    %对应花的数据，误差0.0246
% Dthr=0.24;    %对应酒的数据，误差0.0743
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
% if symbol~=2        
      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(:).^(-4/(mm-1)))'*ones(t,1);%由于d已经取了根号，所以这里的-2/（m-1）,改为-4/（m-1）
           for s1=1:t
             u(i,s1)=(d(s1)^(-4/(mm-1)))/sss;
           end
      end
% else
%         for m1=1:i
%             d=[];
%             [t,tl]=size(Cc);
%             for s1=1:t
%               su=0;
%               su=(data(m1,:)-Cc(s1,:))*(data(m1,:)-Cc(s1,:))';
%               d(s1)=sqrt(su);
%             end
%             %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             for s1=1:t
%               if data(m1,:)==Cc(s1,:);
%                   u(m1,:)=0;
%                   u(m1,s1)=1;
%               end
%             end
%             if  1==isempty(find(1==u(m1,:)))
%                sss=(d(:).^(-4/(mm-1)))'*ones(t,1);%由于d已经取了根号，所以这里的-2/（m-1）,改为-4/（m-1）
%                  for s1=1:t
%                    u(m1,s1)=(d(s1)^(-4/(mm-1)))/sss;
%                  end
%             end
%         end
%         symbol=1;
%         u
% 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;%新建一个类，这个样本为新类的中心
        %i=i-1;%只要新建立了类就将i值减1，这是为了让它能够更新隶属度矩阵的值
        symbol=2;%意思创建了新类
    else
        Ru(nu)=S/2;
        if Ru(nu)<Dthr
        ls=d(nu)/Ru(nu);
       %修改中心坐标
        for f1=1:n
        Cc(nu,f1)=data(i,f1)-(data(i,f1)-Cc(nu,f1))/ls;
        end
        result(i)=nu;
       % i=i-1;
        symbol=2;%代表类中心和类半径修改了，要对隶属度矩阵进行更新
        else    %如果Ru(nu)>=Dthr，则创建新类
              Cc(t+1,:)=data(i,:);
              Ru(t+1)=0;
              result(i)=t+1;%这个样本为新建类的第一个样本
         %     i=i-1;%只要新建立了类就将i值减1，这是为了让它能够更新隶属度矩阵的值
              symbol=2;%意思创建了新类
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%

   if symbol==2
       u=zeros(m,10);
         for m1=1:i
              d=[];
            [t,tl]=size(Cc);
            for s1=1:t
              su=0;
              su=(data(m1,:)-Cc(s1,:))*(data(m1,:)-Cc(s1,:))';
              d(s1)=sqrt(su);
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%
            for s1=1:t
              if data(m1,:)==Cc(s1,:);
                  u(m1,:)=0;
                  u(m1,s1)=1;
              end
            end
            if  1==isempty(find(1==u(m1,:)))
               sss=(d(:).^(-4/(mm-1)))'*ones(t,1);%由于d已经取了根号，所以这里的-2/（m-1）,改为-4/（m-1）
                 for s1=1:t
                   u(m1,s1)=(d(s1)^(-4/(mm-1)))/sss;
                 end
            end
        end
        symbol=1;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else
    result(i)=sub;
    symbol=1;%没有类中心被修改，也没有增加新的类
end
end
U=u
% 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,:))';
        z=z+u(i,j)*(ts);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%用于测量类内样本到中心的距离
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=dis(:);
p=length(dis);
chaoguo=0;
for i=1:p
    if dis(i)>Dthr
        chaoguo=100000;
    end
end


% Cc
% u
 Ru
 z


S=measure(data,U,Cc)
chaoguo


figure(45)
zhongbiao=Cc;
%%%%%%计算各样本的类别标签：Cc中类中心坐标是按行存储的，
[tt1,tt2]=size(data);
[dd1,dd2]=size(Cc);
leiad=zeros(tt1,1);
for i=1:tt1
    juli=zeros(1,dd1);
    for j=1:dd1
        juli(j)=sum((data(i,:)-Cc(j,:)).^2);
    end
    [xxxxxx,leiad(i)]=min(juli);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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');