%  模糊聚類分析
%  模糊ISODATA聚類分析法
%  驗證程序正確性(樣本：P364例)

clc;
clear;
x=[0.80 17 0.22 0.67 0.09; 9.42 18 0.06 1.00 0.14;
   5.91 11 0.36 1.00 0.21; 1.12 17 0.52 0.67 0.12;
   2.96 17 0.57 1.00 0.02; 2.42 11 0.54 1.00 0.01;
   0.99 13 0.23 0.63 0.06; 1.00 13 0.49 0.60 0.02;
   1.26 13 0.55 0.69 0.15; 1.05 16 0.30 0.71 0.11;
   1.06 12 0.43 0.67 0.02; 1.45 15 0.25 0.92 0.08;
   1.21 12 0.24 0.97 0.04; 2.28 15 0.16 1.00 0.01;
   2.25 12 0.18 1.00 0.05; 2.58 15 0.19 1.00 0.08;
   3.02 13 0.16 1.00 0.05; 3.55 15 0.31 1.00 0.27;
   3.79 13 0.31 0.98 0.11; 1.05 13 0.29 0.80 0.02];

%  數(shù)據(jù)規(guī)格化
[n,m]=size(x);
xmax=max(x);xmin=min(x);
for i=1:n
    for j=1:m
        y(i,j)=(x(i,j)-xmin(j))/(xmax(j)-xmin(j));
    end
end

%  分類數(shù)c=3；初始模糊分類矩陣R(c,n)，q=2；精確度e=0.001
c=3;q=2;e=0.000001;p=12;
R=[0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.4 0.4 0.5 0.5 0.5 0.5 0.6 0.6;
   0.2 0.1 0.3 0.4 0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 0.1 0.2 0.1 0.2 0.3 0.4 0.1 0.2;
   0.7 0.8 0.6 0.5 0.7 0.6 0.5 0.4 0.6 0.5 0.4 0.3 0.5 0.4 0.4 0.3 0.2 0.1 0.3 0.2];
V=ones(c,m);
Rl=R;b=1;run_time=0;
%while run_time<10000
while b>e
    run_time=run_time+1
    R=Rl;
    %  用R求V
    Vp=R.^q;
    for i=1:c
        vh=zeros(1,m);vl=0;
        for k=1:m
            vh=vh+Vp(i,k)*y(k,:);
            vl=vl+Vp(i,k);
        end
        V(i,:)=vh/vl;
    end
    %  用V求Rl
    %  參考過渡量：XV(k,i)=||xk-Vi||^(1/(q-2))
    for k=1:n
        for i=1:c
            xvp=0;
            for j=1:m
                xvp=xvp+(abs(y(k,j)-V(i,j)))^p;
            end
            XV(k,i)=xvp^(1/(p*(q-1)));
        end
    end
    for i=1:c
        for k=1:n
            Rlp=0;
            for j=1:c
                Rlp=Rlp+XV(k,i)/XV(k,j);
            end
            Rl(i,k)=1/Rlp;
        end
    end                
    %  門限計算
    A=Rl-R;
    B=(abs(A).^p).^(1/p);
    C=max(B);
    b=max(C);
end

%   求聚類指標Fc與Hc
Fcp=0;Hcp=0;
for i=1:c
    for j=1:n
        Fcp=Fcp+Rl(i,j)^2;
        Hcp=Hcp+Rl(i,j)*logm(Rl(i,j));
    end
end
Fc=Fcp/n;
Hc=Hcp/(-n);

%   判斷聚類R(c,m)
No1=0;No2=0;No3=0;No4=0;No5=0;
for i=1:n
    num(i)=find(R(:,i)==max(R(:,i)));
    switch num(i)
        case 1
            No1=[No1,i];
        case 2
            No2=[No2,i];
        case 3
            No3=[No3,i];
        case 4
            No4=[No4,i];
        case 5
            No5=[No5,i];
        otherwise
            disp('Error!');
    end
end
No1
No2
No3
No4
No5
V
Fc
Hc