ISOData（迭代自组织分析算法）

代码实现

clear all;close all;clc;%   一、载入待分类的样本模式集load mydata.mat%设置控制参数C = 3 ;%预期类数Nc = 1 ;%初始聚类中心个数Tn =15;%每一类的最少模式数目Ts = 5;%类内距离标准差上界Td = 24;%类间最小距离下界L = 2;%在每次迭代中可以合并的类的最多对数I = 10;%允许的最大迭代次数Ip = 1;%初始迭代次数%获取样本的长度N = size(X);N = N(1);%选取初始聚类中心R=randperm(N);Sample(:,1) = X;Sample(:,2) = Y;Z = zeros(Nc,2);%聚类中心for i = 1:Nc    Z(i,:) = Sample(R(i),:);end%绘制样本% grid on;axis([-10 50 -10 50]);hold on% plot(X,Y,'mo');%开始迭代while(Ip <= I)    t = 1;    while(t)        %   二、按最短距离归类        di = zeros(N,Nc);        for i = 1:N            for j = 1:Nc                di(i,j) = norm(Sample(i,:) - Z(j,:));            end        end        dn = zeros(N,1);        for i = 1:N            [m ,n] = min(di(i,:));%m为最短距离，n为所在位置            dn(i,1) = n;%归类        end        %   三、根据Tn，即每类中允许的最少模式        ni = zeros(Nc,1);%各个聚类中心含有的模式数目        for i = 1:Nc            ni(i,1) = length(find(dn == i));%个数            if ni(i,1) < Tn                Nc = Nc-1;                Z(i,:) = 0;            end        end        v = find(Z(:,1)~=0);        K = zeros(Nc,2);        for i = 1:Nc            K(i,:) = Z(v(i),:);        end        Z = K;%合并后的聚类中心        if Nc >= j  %相等时，说明每类样本数满足参数要求            break        end    end    %设置颜色    colour = zeros(Nc,3);    for i = 1:Nc        colour(i,:) = rand(1,3);    end    %  四、计算分类后的参数：更新聚类中心，类内平均距离，总体平均距离    %更新聚类中心    dj = zeros(Nc,1);%类内平均距离    for i = 1:Nc        v = find(dn == i);%位置        u = ni(i,1);%个数        M = zeros(u,2);%用来存储所属同一类的点        D = zeros(u,1);%用来存储类内距离        for j = 1:u            M(j,:) = Sample(v(j),:);            D(j,1) = norm(M(j,:) - Z(i,:));        end        P = plot(M(:,1),M(:,2),'d');        set(P,'color',colour(i,:));        K(i,:) = mean(M);%新的聚类中心        %计算各类中模式到类心的平均距离dj        dj(i,1) = mean(D);%类内平均距离     end    Z1 = K;    %总体平均距离d    d = sum(ni.*dj)/N;    %  五、依据Ip、Nc判断停止，分裂，合并    if Ip >= I        Td = 0;        jump  = 2;%转至9    elseif Nc <= C/2        jump  = 1; %转至6    elseif Nc >= 2*C        jump  = 2;%转至9    elseif Nc > C/2 && Nc < 2*C        if mod(Ip,2) == 1  %如果Ip为奇数            jump  = 1; %转至6，分裂        else            jump  = 2;%转至9 ,合并        end    end    switch  (jump)        case 1    %%%%%%%%%分裂            %  六、计算类内距离标准差矢量si, 矢量维数:n = 2;类编号:j = Nc;            si = zeros(Nc,2);            for i = 1:Nc                v = find(dn == i);%位置                u = ni(i,1);%个数                d1 = zeros(u,1);%X与Z(i,1)做差                d2 = zeros(u,1);%X与Z(i,2)做差                for j = 1:u                    d1(j,1) = (X(v(j),:)-Z1(i,1))^2;                    d2(j,1) = (Y(v(j),:)-Z1(i,2))^2;                end                si(i,1) = mean(d1)^(1/2);                si(i,2) = mean(d2)^(1/2);            end            %  七、对于每一聚类，求出类内距离标准差矢量si中的最大分类simax            simax = zeros(Nc,1);            for i = 1:Nc                [m,n] = max(si(i,2));                simax(i,1) = m;            end            %  八、simax与类内允许的最大类内标准差Ts进行比较            for i = 1:Nc                if simax(i) > Ts                    if (dj(i)>d  &&  ni(i)>2*(Tn+1)) || Nc <= C/2                        Nc = Nc+1;                        K(i,1) = Z1(i,1)+ 0.5*simax(i);                        K(Nc,1) = Z1(i,1)- 0.5*simax(i);                        K(Nc,2) = Z1(i,2);                    end                end            end            Z1 = K;%分裂后的            if i < Nc  %判断是否进行了分裂                Ip = Ip+1;                Z = Z1;            else                %  九、计算各对聚类中心间的距离（类间距离）Dij                Dij = zeros(Nc-1,Nc);                for i = 1:Nc-1                    for j = i+1:Nc                        Dij(i,j) = norm(Z1(i,:) - Z1(j,:));                    end                end                %  十、依据Td判断合并 Dij<Td 合并                k = 1;                dij = zeros(k,1);                for i = 1:Nc-1                    for j = 1:Nc                        if Dij(i,j) > 0 && Dij(i,j) < Td;                            dij(k,1) = Dij(i,j);                            k = k+1;                        end                    end                end                dij = sort(dij);%排序                %根据可以合并的最多对数L进行合并                g = size(dij);                g = g(1);                Zi = zeros(L,1);Zj = zeros(L,1);                for i = 1:L                    if dij(i) == 0;                        break                    end                    [Zi(i,1),Zj(i,1)] = find(Dij == dij(i));                    if  i == 1 || Zi(i,1) ~= Zi(i-1,1)                        Z1(Zi(i),:) = (ni(Zi(i),1)*Z1(i,:)+ni(Zj(i),1)*Z1(j,:))/(ni(Zi(i),1)+ni(Zj(i),1));                        Z1(Zj(i),:) = 0;                        Nc = Nc-1;                    end                    g = g-1;                    if g == 0                        break                    end                end                v = find(Z1(:,1)~=0);                K = zeros(Nc,2);                for i = 1:Nc                    K(i,:) = Z1(v(i),:);                end                Z1 = K;%合并后的聚类中心                %  十一、判断结束、或调整参数、或继续迭代                if Ip >= I || (isequal(size(Z1),size(Z)) && abs(sum(sum(Z1-Z)))<0.0001 )                    Z = Z1;                    break                else                    Ip = Ip+1;                end            end        case 2   %%%%%%%%%合并            %  九、计算各对聚类中心间的距离（类间距离）Dij            Dij = zeros(Nc-1,Nc);            for i = 1:Nc-1                for j = i+1:Nc                    Dij(i,j) = norm(Z1(i,:) - Z1(j,:));                end            end            %  十、依据Td判断合并 Dij<Td 合并            k = 1;            dij = zeros(k,1);            for i = 1:Nc-1                for j = 1:Nc                    if Dij(i,j) > 0 && Dij(i,j) < Td;                        dij(k,1) = Dij(i,j);                        k = k+1;                    end                end            end            dij = sort(dij);%排序            %根据可以合并的最多对数L进行合并            g = size(dij);            g = g(1);            Zi = zeros(L,1);Zj = zeros(L,1);            for i = 1:L                if dij == 0                    break                end                [Zi(i,1),Zj(i,1)] = find(Dij == dij(i));                if  i == 1 || Zi(i,1) ~= Zi(i-1,1)                    Z1(Zi(i),:) = (ni(Zi(i),1)*Z1(i,:)+ni(Zj(i),1)*Z1(j,:))/(ni(Zi(i),1)+ni(Zj(i),1));                    Z1(Zj(i),:) = 0;                    Nc = Nc-1;                end                g = g-1;                if g == 0                    break                end            end            v = find(Z1(:,1)~=0);            K = zeros(Nc,2);            for i = 1:Nc                K(i,:) = Z1(v(i),:);            end            Z1 = K;%合并后的聚类中心            %  十一、判断结束、或调整参数、或继续迭代            if Ip >= I || ( isequal(size(Z1),size(Z)) && abs(sum(sum(Z1-Z)))<0.0001 )                Z = Z1;                break            else                Ip = Ip+1;            end    endend%显示距离中心for i = 1:Nc    plot(Z(i,1),Z(i,2),'kp ');end

效果图