isodata.m

模式识别下的ISODATA算法 例子为IRIS数据库数据

clear
clc
kl
s='y';ii=0;ccc=0;
while s~='n'&ccc<10;
ii=ii+1;ccc=ccc+1

clear si sil si2
data1=xlsread('data1.xls'); 
data2=xlsread('data2.xls');
data3=xlsread('data3.xls');%从workspace读入3类50组数据以备分析.
num=randperm(150);
data=[data1;data2;data3];%xlsread('Y.xls');%%取所有数据为训练数据.
n=length(data);
X=data(num,:);%随机排列训练数据.
si=(num<=50)+((num<=100)&(num>50))*2+(num>100)*3;%标记待训练数据中各数据的类别状态，x∈datai,标为i.
%－－－－－－－－－－－－－－准备训练数据并标记－－－－－－－－－－－－－－－－－

K=3;%input('Please input the K:');        %期望的聚类数.
Min_N=13;%input('Please input the min_N:');%聚类中的最少样本数.
Min_S=.56;%input('Please input the min_S:');%标准偏差参数.
Min_C=1.2;%input('Please input the min_C:');%两聚类中心的最小距离.
L=1;%input('Please input the L:');        %每次迭代允许合并的最大聚类对数.
I0=100;%input('Please input the I:');        %允许迭代的次数.
c=5;%input('Please input the c:');        %初始聚类数.
a0=.5;%input('Please inout the a:');        %分裂距离因子，0<a<1.
m=X(1:c,:);m0=[m';si(1:c)]';            %初始聚类中心.
%－－－－－－－－－－－－－－设置初始数据及控制参数(1)－－－－－－－－－－－－－

count=0;%迭代次数计数器.
%－－－－－－－－－－－－－－开始ISODATA迭代算法－－－－－－－－－－－－－－－－
N=1;N0=0;m1=zeros(c,4);min_size=c;
tic
while ((N(1:min(length(N),length(N0)))~=N0(1:min(length(N),length(N0))))|sum(sum(m(1:min_size,:)~=m1(1:min_size,:))))&count<I0
    N0=N;m1=m;
    clear d N av_d dd;
    cc=(1:c);
    for i=1:n
        for j=1:c
            d(j)=norm(X(i,:)-m(j,:));%计算欧式距离.
        end
        f=find(d==min(d));
        sil(i)=f(1);%分类后的数据类标识，作用同si.
    end
%%－－－－－－－－－－以最小欧式距离将所有样本数据分类(2)－－－－－－－－－－－－    
    for i=1:c
        N(i)=length(find(sil==i));
    end
    %N
    N_th=find(N<Min_N);
    if isempty(N_th)~=1
        for i=1:length(N_th)
            TH=find(sil==N_th(i));
            sil(TH)=0;   %更新后的分类标识.
        end
    end
    c=c-length(N_th);  %更新后的聚类数.
    cc(N_th)=[];
    clear m;
%%－－－－－－－－－－－－取消基数小于Min_N的类(3)－－－－－－－－－－－－－－－ 
    for i=1:c
        ff=find(sil==cc(i));
        m(i,:)=mean(X(ff,:));%计算新的聚类中心.
        for j=1:length(ff)
            d1(j)=norm(X(ff(j),:)-m(i,:));
        end
        av_d(i)=mean(d1);%计算各类中样本至中心的平均距离.
    end
    av_d;
%%－－－－－－－－计算更新的聚类中心和各类样本的平均距离av_d(4),(5)－－－－－－  
    to_av_d=(sum(av_d.*N(cc)))/(sum(N(cc)));
%%－－－－－－－－计算所有样本距其相应中心的平均距离to_av_d,(6)－－－－－－－－
    clear S
    if (c<=K/2)| ((mod(count+1,2)==1)&(c<=2*K))
        for i=1:c
            Y=X((find(sil==cc(i))),:);
            clear mm;
            for j=1:length(Y)
                mm(j,:)=m(i,:);
            end
            S(i,:)=sqrt(mean((Y-mm).^2));%每一类各分量的标准偏差S.
        end
        S_max=max(S'); %每一类中偏差最大的分量S_max
%%－－－－计算每一类各分量的标准偏差s及每一类中偏差最大的分量S_max,(8),(9)－－－
        f3=find(S_max>Min_S);
        nu=c;
        if isempty(f3)==0
            for i=1:length(f3)
                if (nu<=K/2)|(N(f3(i))>2*(Min_N+1))&(av_d(f3(i))>to_av_d)
                    m(c+1,:)=m(f3(i),:)+a0*S_max(f3(i));
                    m(f3(i),:)=m(f3(i),:)-a0*S_max(f3(i));
                    N(c+1)=fix(N(f3(i))/2);
                    N(f3(i))=fix(N(f3(i))/2)+rem(N(f3(i)),2);
                    c=c+1; %类分裂后新的聚类个数.
                end
            end
            m; %类分裂后新的聚类中心.
        end
%%－－－－－－－－－－－－－－－分裂过程(10)－－－－－－－－－－－－－－－－－－ 
    end
%%－－－－－－－－－－－－－－－判断是否需要考虑分裂(7)－－－－－－－－－－－－－
    for i=1:c-1
        for j=i+1:c
            dd(i,j)=norm(m(i,:)-m(j,:)); %两中心间距.
        end
    end
%%－－－－－－－－－－－－对所有中心，求两两间距离dd.(11)－－－－－－－－－－－－－
    clear f4;
    if isempty(find(dd>0&dd<Min_C))==0
        [f4(:,1),f4(:,2)]=find(dd>0&dd<Min_C);
        d2=find((dd<Min_C)&(dd>0));
        f5=sortrows([(dd(d2))';(f4(:,1))';(f4(:,2))']');
        l=min(L,length(f4(:,1)));
        for i=1:l
            m(f5(i,2),:)=(m(f5(i,2),:)*N(f5(i,2))+m(f5(i,3),:)*N(f5(i,3)))/(N(f5(i,2))+N(f5(i,3)));
            N(f5(i,2))=N(f5(i,2))+N(f5(i,3));
        end
        i=1:l;
        m(f5(i,3),:)=[];%合并后新的中心向量.
        N(f5(i,3))=[] %合并后新的类元素个数.
        siz_m=size(m);
        c=siz_m(1);%合并后新的类个数.
    end
%%－－－－－－－－比较dd和Min_C，合并后求新中心和聚类个数.(12)(13)－－－－－－－    
    count=count+1; 
    m_size=size(m);
    m1_size=size(m1);
    min_size=min(m_size(1),m1_size(1));
end
time(ii)=toc;
I(ii)=count;
%－－－－－－－－－－－－－－ISODATA迭代算法结束－－－－－－－－－－－－－－－－ 
sil
clear a d6 b N0 m1
for i=1:c
    f10=find(sil==i);
    f11=si(f10);
    for j=1:c
        n1(j)=length(find(f11==j));
    end
    n2=find(n1==max(n1));
    a(i,n2(1))=1;
end
a; %分类标识转换矩阵.
[d6(:,1),d6(:,2)]=find(a==1);
for i=1:length(d6(:,1))
    b=find(sil==d6(i,1));
    sil(b)=d6(i,2)+10;
    N0(d6(i,2))=N(d6(i,1));
    m1(d6(i,2),:)=m(d6(i,1),:);
end
a_size=size(a);
if a_size(1)<length(N)
    N0(a_size(1)+1:length(N))=N(a_size(1)+1:length(N));
end
m=m1;
zero_N=N0==0;
if zero_N<3
    N0(zero_N)=[];
    m(zero_N,:)=[];
end

%clear N
if length(N0)<3
    Num(ii,:)=[N0,zeros(1,3-length(N0))]; %最终各类中元素个数.
else
    Num(ii,:)=N0(1:3);
end   
sil=sil-10;
%－－－－－－－－－－－－－－将最终分类结果sil对照si做标识－－－－－－－－－－－
error_rate(ii)=1-((sum(si==sil))/n); %计算分类错误率.
%M=[mean(data1);mean(data2);mean(data3)]; %理想的各类中心.
M=[mean(data(1:50,:));mean(data(51:100,:));mean(data(101:150,:))]; %理想的各类中心.
for i=1:n
    clear d f
    for j=1:3
        d(j)=norm(X(i,:)-M(j,:));%以理想中心计算欧式距离.
    end
    f=find(d==min(d));
    si2(i)=f(1);%以理想中心分类后的数据类标识.
end
s_error_rate(ii)=1-(sum(si==si2))/n;%理想中心分类错误率.

%m_size=size(m);
%for i=1:m_size(1)
 %   clear d;
 %   F=find(sil==i);
 %   XX=X(F,:);
 %   for j=1:length(F)
 %       d(i,j)=(norm(XX(j,:)-m(i,:)))^2;
 %   end
%end
%J(ii)=sum(sum(d));    %计算分类误差平方和.
%s_J(ii)=sum(sum(dd));%理想中心分类的误差平方和.
%－－－－－－－－－－－－－－基于最小欧式距离的聚类分析结束－－－－－－－－－－－

s=a;%input('Anykey to continue,Space to exit.','s')
end

m
Num
I
time
error_rate
%J

%plot(si,'or');title('实际数据分类分布');xlabel('数据次序'),ylabel('类别');
%grid %描绘实际数据分类分布.
%figure(2)
%plot(sil,'or');title('ISODATA算法分类结果');xlabel('数据次序'),ylabel('类别');
%grid %描绘ISODATA算法分类结果.
%figure(3)
%plot3(data1(:,1),data1(:,3),data1(:,4),'or',data2(:,1),data2(:,3),data2(:,4),'ok',data3(:,1),data3(:,3),data3(:,4),'og');
%title('实际数据的三维分量空间分布');xlabel('第一维'),ylabel('第三维'),zlabel('第四维')
%grid %描绘实际数据的三维分量空间分布.
%figure(4)
%f_N=find(Num(:,3)~=0);
%KK=.02:.02:ccc;
K=1:ccc;%(f_N');
%plot(K,Num);title('各类中数据个数');xlabel('a'),ylabel('个数'); 
%grid %描绘K取不同值时所分类各类中数据个数.
%figure(5)
%plot(K,I);title('迭代次数');xlabel('a'),ylabel('次数'); 
%grid %描绘K取不同值时所需的迭代次数.
%figure(6)
s_error=0.0667+zeros(1,length(K));
%s_J=s_error-0.067+43.4850;
%plot(K,error_rate,'b',K,s_error,'r');title('错误率');xlabel('a'),ylabel('error_rate'); 
grid %描绘K取不同值时的分类错误率.
%figure(7)
%plot(K,J,'b',K,s_J,'r');title('误差平方和');xlabel('a'),ylabel('J'); 
%grid %描绘K取不同值时的误差平方和.
%figure(8)
plot(K,time);title('消耗时间');xlabel('a'),ylabel('time/second'); 
%grid %描绘K不同时ISODATA算法分类的时间分布.
%figure(9);title('分类错误和正确的个数');xlabel('a'),ylabel('个数'); 
%hist(J)
%grid %描绘在100次实验中由于聚类中心初始随机性而分类错误和正确的个数.