MATLAB实现模糊聚类与模糊综合评价(三个程序,其中许峰教授做的)

function Fuzzy_method_all

 

function Fuzzy_method_all
% 本程序通过一个较为复杂的例子演示了如何调用相关子程序实现模糊聚类与模糊综合评价。
clear,clc,format short g
% 原始数据矩阵.
%   人口    GDP      工业   财政收入  财政支出 固定资产 居民储蓄 环境治理
A0=[%花色苷	酒石酸(g/L)	苹果酸(g/L)	柠檬酸(g/L)	多酚氧化酶活力	褐变度	DPPH自由基
408.028 	2.060 	18.210 	1.830 	33.753 	1119.853 	0.4301 
224.367 	9.930 	4.750 	0.770 	30.904 	762.525 	0.4644 
157.939 	8.080 	2.960 	1.050 	19.303 	266.640 	0.4090 
79.685 	    3.770 	5.230 	0.550 	15.534 	72.905 	0.2655 
120.606 	9.490 	3.770 	1.440 	31.536 	143.513 	0.3961 
46.186 	    2.830 	2.210 	0	36.774 	115.943 	0.2750 
60.767 	    5.820 	7.740 	0.540 	25.591 	433.751 	0.1756 
241.397 	5.710 	13.550 	2.510 	50.434 	1305.595 	0.4148 
240.843 	13.230 	4.120 	1.100 	16.869 	424.108 	0.6658 
44.203 	    2.450 	2.300 	0.240 	10.427 	459.569 	0.3255 
7.787 	    9.290 	8.610 	1.900 	14.260 	91.468 	0.2790 
32.343 	    6.080 	5.330 	1.130 	21.080 	132.216 	0.1973 
65.324 	    4.300 	0.830 	1.150 	28.076 	99.881 	0.4406 
140.257 	5.730 	4.120 	1.630 	41.577 	991.046 	0.3597 
52.792 	    6.230 	3.630 	2.060 	25.743 	157.997 	0.2189 
60.660 	    9.030 	7.280 	2.380 	13.648 	529.969 	0.2367 
59.424 	    5.880 	5.110 	0.880 	17.174 	129.581 	0.3585 
40.228 	    3.60 	5.590 	0.520 	27.077 	158.870 	0.2256 
115.704 	5.560 	4.270 	0.130 	30.408 	202.962 	0.3796 
23.523 	3.510 	0.920 	0.440 	12.439 	89.770 	0.2819 
89.282 	15.510 	2.930 	2.380 	18.123 	194.262 	0.3793 
74.027 	6.490 	7.730 	0.770 	21.824 	417.665 	0.2837 
172.626 	4.080 	5.200 	0.390 	16.406 	427.028 	0.5725 
144.881 	8.360 	4.600 	1.700 	15.066 	144.729 	0.2830 
49.643 	2.870 	2.480 	0.160 	14.280 	140.946 	0.3509 
58.469 	7.150 	1.400 	0.820 	32.026 	82.359 	0.3172 
34.190 	6.230 	1.390 	1.260 	23.035 	592.199 	0.2649 ]




%[l,k]=find(A0==0);
%ratio=mean(A0([1:l-1,l+1:end],k)./A0([1:l-1,l+1:end],2));
%A0(l,k)=A0(l,2)*ratio;
A=[]; % A记录七个人均指标.
for i=2:size(A0,2)
    A=[A A0(:,i)./A0(:,1)];
end
Dynamic_clustering(A,7);
B=round(A); % A取整为B,B将写入正文.
[m n]=size(A);
% 标准化原始数据矩阵 A,标准化矩阵仍记为 A.
A=Standard(A);
% 将标准化矩阵 A 写入 excel,文件名为 result.xls.
you=fopen('result.xls','w');
fprintf(you,'       \n');
fprintf(you,'标准化矩阵\n');
geshi1=[]; % 写入excel的书写格式.
for i=1:m-1
    geshi1=[geshi1 '%f\t'];
end
geshi1=[geshi1 '%f\n'];
for i = 1:n
    fprintf(you,geshi1,A(i,:));
end;
% 第一个问题:模糊聚类
R=Fuzzy_similarity_matrix(A,7);
% 再应用平方法计算R的传递闭包,仍记闭包为R.
R=Transtive_closure(R);
% 将传递闭包矩阵R写入resultforbook.xls.
fprintf(you,'       \n');
fprintf(you,'模糊等价矩阵的传递闭包\n');
geshi2=[]; % 写入excel的书写格式.
for i=1:m-1
    geshi2=[geshi2 '%f\t'];
end
geshi2=[geshi2 '%f\n'];
for i = 1:m
    fprintf(you,geshi2,R(i,:));
end;
% 求R的lamda截矩阵,这需要分类水平lamda:
lamda=Classification_level(R);
% 将分类水平 lamda 写入 result_forbook.xls.
fprintf(you,'       \n');
fprintf(you,'分类水平\n');
fprintf(you,geshi2,lamda);
for p=1:length(lamda)
    %显示 lamda(p) 水平上分类情况.
    M(:,:,p)=R>=lamda(p);
    [cls nmb]=Computing_cls(M(:,:,p));
    disp(' '); %显示一空行是为了方便阅读显示内容.
    disp(['在分类水平 ',num2str(lamda(p)),'上分为 ',num2str(nmb),'类: ']);
    cities{1}='1 '; cities{2}='2 '; cities{3}='3 '; cities{4}='4 '; cities{5}='5 '; 
    cities{6}='6 '; cities{7}='7 '; cities{8}='8 '; cities{9}='9 '; cities{10}='10 '; 
    cities{11}='11 '; cities{12}='12'; cities{13}='13 '; cities{14}='14 '; cities{15}='15 '; 
    cities{16}='16 '; cities{17}='17 '; cities{18}='18 '; cities{19}='19 '; cities{20}='20 '; 
    cities{21}='21 '; cities{22}='22 '; cities{23}='23 '; cities{24}='24 '; cities{25}='25 ';
    cities{26}='26 '; cities{27}='27 ';
    % cities{28}='成都 '; cities{29}='贵阳 '; cities{30}='昆明 '; cities{31}='西安 '; cities{32}='兰州 '; cities{33}='西宁 '; cities{34}='银川 '; cities{35}='乌鲁木齐 '; 
    for i=1:nmb
        if length(cls{i})==1
            cities{cls{i}}(end)=[];
            disp(['"',cities{cls{i}},'"','  自成一类']);
        else
            ct=[];
            for t=1:length(cls{i})
                ct=[ct cities{cls{i}(t)}];
            end
            ct(end)=[];
            disp([' "',ct,'"',' 归为一类']);
        end
    end
    %求水平lamda(p)上的 F 值--F(p).
    F(p)=F_statistic(M(:,:,p),A);
end
disp('       ')
[v,ind]=max(F);
disp(['在所有分类中分为 ' num2str(ind) ' 类是最合理的.'])
% 第二个问题:模糊模式识别.
% 计算七个类聚类中心.
I=[2 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 26 28 29 30 32 34 35];
II=[25 27];
alpha1=mean(A(I,:)); alpha2=mean(A(II,:));
W=[A(1,:);alpha1;A(3,:);A(24,:);alpha2;A(31,:);A(33,:)];
% C-待识别的甲乙丙三城市的原始数据.
C=[89962 43704 13134 12756 19356 43662 658
   37959 11662 8197  8251  12993 26005 857
   16642 9037  11071 12112 28451 39037 1302];
C=Standard(C);
D=dist(C,W'); % C到聚类中心W的欧氏距离.
disp('    ');
jyb='甲乙丙';
for i=1:size(D,1)
    [val,ind]=min(D(i,:));
    disp(['城市' jyb(i) '属于' num2str(ind) '类.'])
end
% 第三个问题:模糊综合评价.
O=[1 1 1 1 1 1 0]; %虚拟的最优对象O.
E=dist(A,O'); % 35个城市到最优对象聚O的欧氏距离.
[hh ss]=sort(E); % ss记录了35个城市按距离的排序结果.
disp('       ');
disp('距离法给出35个城市排序为: ');
disp( num2str(ss'));
    
function S=Standard(x)
% 本子程序功能是通过标定和压缩对原始数据进行标准化。
n=size(x,2); mu=mean(x); sig=std(x);
for j=1:n
    R(:,j)=(x(:,j)-mu(j))/sig(j);
    S(:,j)=(R(:,j)-min(R(:,j)))/(max(R(:,j))-min(R(:,j)));
end

function r=Fuzzy_similarity_matrix(x,method)
% 本子程序功能是分别用7种方法建立模糊相似矩阵。
[m,n]=size(x);
for i=1:m
    for j=1:m
        ed(i,j)=sqrt(sum((x(i,:)-x(j,:)).^2));
        hd(i,j)=sum(abs(x(i,:)-x(j,:)));
    end
end
for i=1:m
    for j=1:m
        K(i,j)=sum(x(i,:).*x(j,:));
    end
end
for i=1:m
    K(i,i)=0;
end
E=max(max(ed)); H=max(max(hd)); M=max(max(K));
disp('  ');
if method==1
    disp('夹角余弦法建立模糊相似矩阵');
elseif method==2
    disp('相似系数法建立模糊相似矩阵');
elseif method==3
    disp('欧氏距离法建立模糊相似矩阵');
elseif method==4
    disp('海明距离法建立模糊相似矩阵');
elseif method==5
    disp('指数距离法建立模糊相似矩阵');
elseif method==6
    disp('数量积法建立模糊相似矩阵');
elseif method==7
    disp('算术平均最小法建立模糊相似矩阵');
end
for i=1:m
    for j=1:m
        if method==1
            r(i,j)=sum(x(i,:).*x(j,:))/sqrt(sum(x(i,:).^2)*sum(x(j,:).^2));
        elseif method==2
            r(i,j)=sum(abs((x(i,:)-mean(x(i,:))).*(x(j,:)-mean(x(j,:)))))/sqrt(sum((x(i,:)-mean(x(i,:))).^2*sum((x(j,:)-mean(x(j,:))).^2)));
        elseif method==3
            r(i,j)=1-ed(i,j)/E;
        elseif method==4
            r(i,j)=1-hd(i,j)/H;
        elseif method==5
            r(i,j)=exp(-sum(abs(x(i,:)-x(j,:))));
        elseif method==6
            if i==1
                r(i,j)=1;
            else
                r(i,j)=sum(x(i,:).*x(j,:))/M;
            end
        elseif method==7
            mink=[];
            for k=1:n
                mink=[mink min(x(i,k),x(j,k))];
            end
            r(i,j)=2*sum(mink)/sum(x(i,:)+x(j,:));
        end
    end
end

function b=FSM_square(x)
% 本子程序功能是求模糊相似矩阵的平方;输入参数x必须是一个模糊方阵,输出为x的平方。
[m,n]=size(x);
if m~=n
    disp('输入矩阵不是方阵,请重新输入一个方阵');
    return
end
for i=1:n
    for j=1:m
        for k=1:n
            b(k)=min(x(i,k),x(k,j));
        end
        xsquare(i,j)=max(b);
    end
end
b=xsquare;

function bibao=Transtive_closure(x)
% 本子程序功能是求模糊相似矩阵x的传递闭包,以用于模糊聚类。
[m,n]=size(x);
if m~=n
    disp('输入矩阵不是方阵,请重新输入一个方阵');
    return
end
xsquare=FSM_square(x);
while sum(sum(x==xsquare))~=n^2;
    x=xsquare;
    xsquare=FSM_square(x);
end
bibao=xsquare;

function A=Classification_level(x)
% 本子程序功能是求动态模糊聚类的分类水平。
[m,n]=size(x);
B=sort(reshape(x,1,m*n));
M=B(end);
i=1;
while B(i)<M
    L=find(B==B(i));
    B(L(2:end))=[];
    i=i+1;
end
L=find(B==M);
B(L(2:end))=[];
A=B;

function [cls,nmb]=Computing_cls(x)
% 本子程序功能是求存储分类的数组cls及相应分类下包含对象的数量。
m=size(x,1);
for i=1:m
    J=[];
    for j=1:m
        if x(i,j)==1
            J=[J,j];
        end
    end
    cl{i}=J;
end
for i=1:m
    for j=i+1:m
        if length(cl{j})==length(cl{i})
            if sum(cl{j}==cl{i})==length(cl{i})
                cl{i}=[];
            end
        end
    end
end
cls=[];
for i=1:m
    if ~isempty(cl{i})
        cls=[cls cl(i)];
    end
end
nmb=length(cls);
        
function F=F_statistic(r,x)
% 本子程序功能是求F统计量,以确定最佳分类数目。
[cls,nmb]=Computing_cls(r);
m=size(r,1);
ubar=mean(x);
for k=1:nmb
    s=cls{k};
    if length(s)==1
        ukbar(k,:)=x(s,:);
    else
        ukbar(k,:)=mean(x(s,:));
    end
    fz(k)=length(s)*sum((ukbar(k,:)-ubar).^2)/(nmb-1);
    for l=1:length(s)
        fm0(k,l)=sum((x(s(l))-ukbar(k,:)).^2);
    end
    fm(k)=sum(fm0(k,:)/(m-nmb));
end
F=sum(fz)/sum(fm);

function Dynamic_clustering(x,method)
% 本子程序功能是动态聚类图,返回的是原始数据的动态聚类图。
close all;
A=Standard(x);
B=Fuzzy_similarity_matrix(A,method);
R=Transtive_closure(B);
lamda=Classification_level(R);
m=length(lamda);
for i=1:m
    M(:,:,i)=R>=lamda(i);
    cls=Computing_cls(M(:,:,i));
    allcls{i}=cls;
end
rr=[];
for i=1:length(allcls)-1
    ind=[];
    for j=1:i
        ind0=[];
        for k=1:i+1
            ind0=[ind0 ismember(allcls{i+1}{k}(1),allcls{i}{j})];
        end
        ind=[ind;ind0];
    end
    tmp=allcls{i+1};
    for r=1:i
        if sum(ind(r,:))==2
            tp=ind(r,:);
            s=find(tp==1);
            allcls{i+1}{r}=tmp{s(1)};
            allcls{i+1}{r+1}=tmp{s(2)};
            t1=[1:r-1,r+1:i];
            t2=[1:r-1,r+2:i+1];
            for t=1:i-1
                allcls{i+1}{t2(t)}=allcls{i}{t1(t)};
            end
            rr=[rr r];
        end
    end
end
seq0=allcls{end};
seq=[];
for i=1:m
    seq=[seq allcls{end}{i}];
end
figure,hold on,title('动态聚类图');
wide=30; 
high=25;
x0=wide*(m-1); 
y0=high*m;
par=(wide+high)/2;
axis([-2.5*par,x0+1.5*par,-par,y0+par]);
axis('off');
lbx{1}=wide*[0:m-1];
for i=2:m
    lbx{i}(rr(m+1-i))=(lbx{i-1}(rr(m+1-i))+lbx{i-1}(rr(m+1-i)+1))/2;
    for u=1:rr(m+1-i)-1
        lbx{i}(u)=lbx{i-1}(u);
    end
    for v=rr(m+1-i)+1:m+1-i
        lbx{i}(v)=lbx{i-1}(v+1);
    end
end
for i=1:m
    text((i-1)*wide-wide/9,-0.5*high,int2str(allcls{end}{i}));
end
text(-3.45*par,(m+0.5)*high,'分类水平');
text(-3.1*par,0.5*high,num2str(lamda(m)));
text(x0+0.3*par,(m+0.5)*high,'分类数');
text(x0+1.35*par,0.5*high,num2str(m));
for i=1:m
    line([(i-1)*wide (i-1)*wide],[0 high]);
end
for i=2:m
    text(-3.1*par,(i-0.5)*high,num2str(lamda(m+1-i)));
    text(x0+1.35*par,(i-0.5)*high,num2str(m+1-i));
    line([lbx{i-1}(rr(m+1-i)) lbx{i-1}(rr(m+1-i)+1)],[i-1 i-1]*high);
    for j=1:length(lbx{i})
        line([lbx{i}(j) lbx{i}(j)],[i-1 i]*high);
    end
end

2.function Fuzzy_method_all

function Fuzzy_method_all
% 本程序通过一个较为复杂的例子演示了如何调用相关子程序实现模糊聚类与模糊综合评价。
clear,clc,format short g
% 原始数据矩阵.
%   人口    GDP      工业   财政收入  财政支出 固定资产 居民储蓄 环境治理
A0=[1198 78702835 18218601 11171514 12968389 33715013 86984521  388750
    949  43591500 22927300  4170479  5431219 18497987 28110200  799933
    940  20266320  8751622   773736  1282462 10994463 15532427   27727
    349  10136482  3786335   753306   954238  5011273 11839536  284793
    216   9000845  2837701   455202   754985  5498228  4530845   46078
    704  25196338 10088171  1756336  2687649 17903457 19219758  138814
    572  25696699 10588385  1961357  2665249 14694899 18664886  263991
    739  17411922  6840335   715908  1466748  9504180 11745969  284745
    980  20940751  5894994  1182934  1946839  8098657 15638334  346446
   1368 103663700 46701100 15760724 17955660 39250884 94802800 3108523
    607  27737800 11819400  2464392  2624614 16135518 18348800  850030
    666  34415068 15590895  3013888  2754809 14607422 24736900  141461
    560  28744435 14259477  2573799  2926969 15027686 17520388  185041
    470  10737600  4093600   794046  1028869  8167825  6245600   15407
    623  16640515  6601028  1156184  1100589  7323412 13113823  319000
    160  11680229  5603061  1361653  1591253  6620984  6808196  299504
    484  11838973  4481541   681075   933749  6430198  7338474  255279
    603  21850856  8615053  1284388  1469762 10167663 11825655  442928
    749  32065800 15274900  2259904  2367875 14856894 15676197  566252
    692  20134777  9440244  1760266  1876587 10319917 16207242   77486
    819  25907569 10007433  1786021  2578102 13252827 18877333  104608
    631  17989572  6007131  1328345  1671873 10898081 10930400   18404
    761  60738277 22210421  4270831  5067899 16963824 55623554  265927
    197  58135624 28866206  5008827  5714231 12736693 37447000   10793
    672   8701481  2212890   566221   930781  4472211  6814522  178113
    177   3501246   798078   211288   322454  1582308  4069576   10570
   3199  34915700 12341200  3177165  5942543 24518351 29490500  846064];
A0=[2027.960	101.220	393.420	77.610	266.600	723.880	177.370	89.280	24.830	15.740	17.140	6.580	10.860	3.180	5.030	18.850	9.970	70.660	553.106	.251	408.028	2.060	18.210	1.830	33.753	1119.853	.430	23.604	22.019	9.480	3.195	.388	2.559	.248	.394	17.678	5.305	7.873	.365	4.135	208.175	237.668	110.150	127.517	226.500	3.560	5.860	38.660	25.918	182.930	123.600	4.510	78.400	.110	24.070	.780	.260
2128.820	64.430	140.620	71.940	39.260	1560.970	32.380	11.130	24.110	20.690	4.600	9.420	22.820	1.320	1.300	15.450	13.570	74.130	626.478	.062	224.367	9.930	4.750	.770	30.904	762.525	.464	26.875	23.361	13.806	4.889	.453	3.881	.555	.394	27.455	8.511	11.558	1.420	5.967	205.000	229.136	113.498	115.638	228.800	3.950	5.190	44.050	25.986	81.620	98.300	3.830	77.500	.163	26.070	.650	-1.250
8397.280	108.070	222.350	173.080	67.540	7472.280	55.790	75.340	13.180	19.600	7.840	7.820	18.170	2.760	6.270	31.210	16.510	79.870	585.046	.315	157.939	8.080	2.960	1.050	19.303	266.640	.409	21.685	20.373	10.794	4.764	.354	4.254	.156	.394	164.993	19.977	115.732	6.482	22.802	256.190	273.758	132.209	141.549	257.600	3.910	7.160	35.990	28.997	83.130	105.400	5.600	71.800	.170	25.500	1.090	-.620
2144.680	79.390	133.830	158.740	156.720	1182.230	93.230	89.360	46.700	21.940	6.550	15.790	20.750	2.530	6.080	18.910	17.450	72.530	529.823	.097	79.685	3.770	5.230	.550	15.534	72.905	.266	10.698	8.638	4.482	3.412	.287	2.850	.102	.172	26.968	4.183	16.087	2.232	4.465	189.722	237.766	109.316	128.450	203.300	3.290	7.110	28.610	23.721	137.970	174.700	3.260	53.000	.174	25.980	1.840	-.370
1844.000	52.280	145.090	164.050	102.430	816.080	86.830	69.540	18.640	33.670	16.460	30.480	21.690	3.380	5.300	47.770	29.740	166.900	585.613	.041	120.606	9.490	3.770	1.440	31.536	143.513	.396	17.618	14.486	10.275	.637	.234	4.534	.403	.394	6.650	1.980	4.313	.358	6.203	209.663	195.460	99.585	95.875	212.900	3.640	6.650	32.000	24.084	515.460	254.200	2.990	65.600	.270	26.330	.880	-.330
3434.170	68.010	102.420	75.780	80.600	2932.760	18.010	19.390	23.170	12.630	2.650	10.520	7.900	2.460	2.270	10.000	9.770	43.200	536.643	.075	46.186	2.830	2.210	.000	36.774	115.943	.275	10.671	15.173	6.838	2.203	.254	1.850	.100	.394	7.727	1.056	5.765	.907	6.203	244.385	223.817	108.798	115.019	246.100	3.290	9.310	26.430	27.376	202.240	172.000	2.640	71.900	.193	25.160	1.810	-.160
2391.160	65.100	267.760	239.200	208.970	1096.280	74.060	89.560	18.190	46.980	11.360	18.030	33.840	1.860	2.630	26.360	36.750	107.210	487.172	.131	60.767	5.820	7.740	.540	25.591	433.751	.176	9.214	5.619	3.468	.623	.523	.555	.068	.394	9.865	3.171	6.480	.214	6.203	209.861	303.950	142.437	161.513	211.400	3.180	8.140	25.980	26.438	63.610	168.800	4.780	71.500	.141	25.610	2.050	-.380
1950.760	72.090	345.870	44.230	176.020	962.010	150.730	42.630	6.430	20.840	9.800	16.220	19.380	2.580	3.420	15.930	10.110	31.630	558.546	.181	241.397	5.710	13.550	2.510	50.434	1305.595	.415	15.241	22.489	8.483	5.949	5.283	4.534	.549	.116	115.555	11.630	73.211	18.099	12.615	198.849	196.990	94.336	102.654	226.500	2.920	6.470	34.990	25.620	213.090	181.100	6.410	59.600	.260	26.850	.800	-.510
2262.720	72.890	113.940	110.610	110.530	1334.190	95.180	42.800	7.070	33.490	15.150	24.820	26.610	5.190	8.970	32.610	34.740	160.440	700.828	.512	240.843	13.230	4.120	1.100	16.869	424.108	.666	30.114	24.362	20.490	4.907	.423	3.022	.904	.558	58.541	11.063	33.172	2.771	11.535	193.690	194.925	98.701	96.224	203.400	3.740	5.880	34.580	23.761	186.620	138.100	5.310	78.000	.130	23.810	1.440	-.380
1364.140	87.520	114.290	130.870	126.710	477.500	88.140	60.500	5.750	25.310	4.000	10.560	20.050	1.210	2.970	39.040	18.530	125.880	545.305	10.250	44.203	2.450	2.300	.240	10.427	459.569	.326	9.476	16.688	4.631	12.307	.527	11.140	.485	.156	28.748	10.367	10.274	5.746	2.361	167.202	161.421	79.379	82.041	181.200	3.650	6.670	27.160	19.676	255.440	200.800	4.590	71.700	.200	27.100	2.170	-1.120
2355.690	94.420	111.670	141.580	186.520	1150.090	158.360	83.160	33.210	23.930	7.170	13.440	18.070	3.530	4.620	29.090	17.710	255.190	542.662	.076	7.787	9.290	8.610	1.900	14.260	91.468	.279	6.075	4.543	2.517	26.851	.270	26.335	.246	.394	25.575	17.198	23.717	7.403	.975	209.563	237.891	113.952	123.939	210.200	3.530	5.500	38.240	24.527	177.830	118.800	3.410	58.400	.102	28.030	12.150	3.870
2556.790	63.320	71.680	69.350	47.890	2127.910	36.840	13.980	17.060	19.040	3.960	16.050	14.600	5.170	5.910	6.210	3.290	15.490	493.460	.065	32.343	6.080	5.330	1.130	21.080	132.216	.197	12.059	7.169	3.897	.696	.100	.545	.052	.394	2.480	1.679	23.717	.801	6.203	247.659	262.155	124.661	137.493	261.100	3.430	8.540	30.580	27.614	191.950	187.700	2.400	63.300	.243	26.570	2.040	.010
1416.110	54.300	110.630	80.650	72.320	621.250	41.250	28.390	16.760	37.670	15.110	21.500	36.360	4.020	2.500	39.270	28.770	167.680	606.204	.015	65.324	4.300	.830	1.150	28.076	99.881	.441	14.385	9.822	7.330	10.863	.695	10.050	.117	.394	40.759	9.138	3.088	18.331	10.202	197.857	212.237	110.421	101.816	203.400	3.860	4.340	23.750	23.353	159.970	148.000	4.670	68.100	.160	27.530	1.040	-1.570
1237.810	71.270	56.410	104.500	64.280	677.780	39.090	51.430	22.230	15.730	1.700	12.670	11.660	1.120	3.370	14.760	12.040	62.030	599.829	.060	140.257	5.730	4.120	1.630	41.577	991.046	.360	14.657	13.941	7.809	6.313	.661	5.065	.462	.126	134.638	7.123	113.258	11.158	3.098	191.508	255.335	120.444	134.892	193.900	3.390	5.400	35.900	24.060	209.110	136.300	4.600	66.200	.255	25.410	1.190	-.570
2177.910	85.200	223.120	226.600	172.690	817.570	96.020	100.860	34.220	35.080	10.720	23.910	41.260	2.100	7.540	48.430	31.440	186.070	524.613	.068	52.792	6.230	3.630	2.060	25.743	157.997	.219	11.901	25.417	5.511	.211	.098	4.534	.112	.394	9.718	.472	8.500	.746	6.203	179.107	208.933	105.755	103.178	214.900	3.190	8.570	25.090	25.012	159.310	174.500	2.900	67.700	.213	25.530	1.980	-.010
1553.500	73.340	110.720	110.490	112.050	679.250	46.950	40.080	27.080	23.880	3.970	14.310	19.170	3.560	12.910	42.030	28.130	181.690	583.374	.083	60.660	9.030	7.280	2.380	13.648	529.969	.237	11.214	10.086	9.157	4.556	.335	4.121	.100	.394	8.190	2.161	1.233	3.176	1.620	204.008	189.275	92.785	96.491	205.600	3.300	4.920	41.760	22.346	119.170	109.300	3.790	71.800	.135	26.110	1.330	-.340
1713.650	107.190	95.610	89.150	100.970	806.560	62.420	44.850	26.380	32.110	5.490	20.140	28.520	4.670	10.120	47.160	24.220	175.990	548.833	.056	59.424	5.880	5.110	.880	17.174	129.581	.359	15.336	15.730	8.701	.711	.414	.151	.145	.394	43.812	2.423	37.441	1.880	2.068	212.738	271.504	128.705	142.799	238.200	3.430	8.660	27.510	26.276	446.640	264.100	2.800	71.500	.330	25.400	1.180	-.250
2398.380	32.450	71.490	52.060	38.220	2097.610	24.310	10.010	11.420	10.850	1.590	4.620	4.440	1.340	2.760	6.000	3.830	14.540	513.817	.112	40.228	3.600	5.590	.520	27.077	158.870	.226	7.381	5.388	5.245	.416	.115	4.534	.301	.394	6.516	.741	5.310	.465	6.203	226.032	265.773	125.611	140.162	226.600	3.270	8.030	28.210	26.338	196.010	208.400	2.600	63.100	.160	25.520	2.870	.210
2463.600	72.940	96.130	91.110	68.700	1438.080	165.660	68.050	32.070	41.570	4.440	8.760	30.100	1.920	2.350	29.220	26.230	244.690	544.462	.072	115.704	5.560	4.270	.130	30.408	202.962	.380	17.426	13.700	9.454	3.821	.171	2.879	.177	.592	31.265	.954	24.877	1.670	3.764	205.794	220.333	107.601	112.732	214.900	3.570	6.810	31.540	23.441	173.090	168.800	6.320	67.400	.162	27.190	.800	-1.510
2273.630	73.490	157.320	131.450	177.910	754.890	9.890	100.700	12.620	42.040	17.170	29.960	44.740	3.030	2.790	76.890	59.410	537.280	559.332	.024	23.523	3.510	.920	.440	12.439	89.770	.282	12.677	8.115	8.155	1.545	.229	1.178	.279	.138	9.626	6.891	1.434	1.302	6.203	193.194	227.338	114.537	112.801	209.100	3.810	5.170	40.480	22.933	307.140	334.300	3.150	59.500	.232	27.090	1.960	-.430
6346.830	69.490	180.030	194.110	107.730	5144.810	62.220	16.360	53.640	41.460	13.660	16.120	34.280	2.540	2.230	92.410	41.040	233.250	563.794	.050	89.282	15.510	2.930	2.380	18.123	194.262	.379	16.192	13.613	7.515	7.847	.644	6.437	.632	.135	47.220	1.897	36.037	2.627	6.658	205.794	259.110	121.825	137.285	216.900	3.560	6.780	31.990	26.948	147.660	106.100	4.740	60.400	.108	25.180	1.210	.000
2566.610	110.520	207.260	251.110	237.700	863.990	197.320	70.860	43.770	58.800	11.860	66.930	66.980	6.680	5.250	44.610	43.580	220.580	488.712	.074	74.027	6.490	7.730	.770	21.824	417.665	.284	16.442	12.155	7.846	4.289	.205	3.449	.101	.534	13.800	.737	12.733	.330	6.203	224.147	226.399	105.994	120.405	234.700	3.650	5.970	39.360	25.674	106.610	115.800	3.320	57.400	.147	25.940	1.520	-.070
2380.810	120.860	138.150	159.470	131.540	1341.120	106.710	78.870	35.510	42.270	11.400	18.610	35.540	10.490	6.550	33.950	18.550	48.950	543.574	.097	172.626	4.080	5.200	.390	16.406	427.028	.573	29.704	24.257	24.295	9.968	.690	7.375	.388	1.515	44.748	6.612	24.496	.876	12.764	207.679	212.564	99.539	113.025	208.800	3.390	6.910	30.230	23.383	278.750	219.100	3.840	77.500	.233	26.650	1.380	-.420
1638.830	58.600	160.810	148.580	59.230	797.550	88.560	79.400	34.070	30.320	6.190	11.220	8.600	1.860	2.470	24.420	17.460	79.160	525.820	.033	144.881	8.360	4.600	1.700	15.066	144.729	.283	8.751	14.417	8.206	2.935	.321	2.396	.217	.394	14.380	5.840	4.251	1.020	3.270	201.825	244.512	113.682	130.830	203.300	3.610	7.270	27.980	25.815	517.450	237.400	2.990	76.700	.247	25.970	.900	-.290
1409.700	73.280	130.810	115.850	150.570	479.170	88.030	57.020	41.970	26.150	4.380	15.940	20.850	1.460	3.430	36.190	16.050	122.400	537.084	.064	49.643	2.870	2.480	.160	14.280	140.946	.351	11.502	9.324	5.373	2.129	.196	1.426	.219	.288	30.211	.328	24.614	2.697	2.572	150.337	156.038	80.300	75.738	194.600	3.380	8.530	22.810	18.515	288.690	251.300	4.100	58.500	.220	27.100	1.520	-.920
851.170	59.000	95.660	74.470	47.830	147.700	22.810	28.720	14.940	19.490	4.300	18.420	23.330	5.760	3.980	37.290	33.450	194.540	587.293	.416	58.469	7.150	1.400	.820	32.026	82.359	.317	7.348	3.778	3.383	2.086	.095	1.832	.159	.394	13.917	5.316	6.011	1.815	.775	173.353	197.377	89.120	108.257	195.700	3.680	4.580	42.740	19.758	793.470	245.500	3.350	68.300	.230	28.000	1.090	-.830
1116.610	51.450	132.550	79.060	48.360	418.010	29.250	28.580	9.960	28.740	5.400	10.470	32.750	7.440	4.760	27.460	24.440	149.190	528.331	.091	34.190	6.230	1.390	1.260	23.035	592.199	.265	8.897	10.310	4.711	1.569	.166	1.147	.257	.394	15.981	5.942	5.173	4.865	6.203	196.667	213.216	113.469	99.747	206.900	3.370	6.970	29.670	23.329	282.090	148.700	3.510	59.500	.200	28.790	2.330	-1.230];

A=[];
for i=2:size(A0,2)
    A=[A A0(:,i)./A0(:,1)];
end
Dynamic_clustering(A,1);
B=round(A); % A取整为B,B将写入正文.
[m n]=size(A);
% 标准化原始数据矩阵 A,标准化矩阵仍记为 A.
%A=Standard(A);
% 将标准化矩阵 A 写入 excel,文件名为 result.xls.
you=fopen('result.xls','w');
% 第一个问题:模糊聚类
R=Fuzzy_similarity_matrix(A,1);
% 再应用平方法计算R的传递闭包,仍记闭包为R.
R=Transtive_closure(R);
% 将传递闭包矩阵R写入resultforbook.xls.
fprintf(you,'       \n');
fprintf(you,'模糊等价矩阵的传递闭包\n');
geshi2=[]; % 写入excel的书写格式.
for i=1:m-1
    geshi2=[geshi2 '%f\t'];
end
geshi2=[geshi2 '%f\n'];
for i = 1:m
    fprintf(you,geshi2,R(i,:));
end;
% 求R的lamda截矩阵,这需要分类水平lamda:
lamda=Classification_level(R);
% 将分类水平 lamda 写入 result_forbook.xls.
fprintf(you,'       \n');
fprintf(you,'分类水平\n');
fprintf(you,geshi2,lamda);
for p=1:length(lamda)
    %显示 lamda(p) 水平上分类情况.
    M(:,:,p)=R>=lamda(p);
    [cls nmb]=Computing_cls(M(:,:,p));
    disp(' '); %显示一空行是为了方便阅读显示内容.
    disp(['在分类水平 ',num2str(lamda(p)),'上分为 ',num2str(nmb),'类: ']);
    cities{1}='品种1 '; cities{2}='品种2 '; cities{3}='品种3 '; cities{4}='品种4 '; cities{5}='品种5 '; 
    cities{6}='品种6 '; cities{7}='品种7 '; cities{8}='品种8'; cities{9}='品种9 '; cities{10}='品种10 '; 
    cities{11}='品种11 '; cities{12}='品种12 '; cities{13}='品种13 '; cities{14}='品种14 '; cities{15}='品种15 '; 
    cities{16}='品种16 '; cities{17}='品种17 '; cities{18}='品种18 '; cities{19}='品种19 '; cities{20}='品种20 '; 
    cities{21}='品种21 '; cities{22}='品种22 '; cities{23}='品种23 '; cities{24}='品种24 '; cities{25}='品种25 ';
    cities{26}='品种26 '; cities{27}='品种27 '; 
    for i=1:nmb
        if length(cls{i})==1
            cities{cls{i}}(end)=[];
            disp(['"',cities{cls{i}},'"','  自成一类']);
        else
            ct=[];
            for t=1:length(cls{i})
                ct=[ct cities{cls{i}(t)}];
            end
            ct(end)=[];
            disp([' "',ct,'"',' 归为一类']);
        end
    end
    %求水平lamda(p)上的 F 值--F(p).
    F(p)=F_statistic(M(:,:,p),A);
end
disp('       ')
[v,ind]=max(F);
disp(['在所有分类中分为 ' num2str(ind) ' 类是最合理的.'])
    
function S=Standard(x)
% 本子程序功能是通过标定和压缩对原始数据进行标准化。
n=size(x,2); mu=mean(x); sig=std(x);
for j=1:n
    R(:,j)=(x(:,j)-mu(j))/sig(j);
    S(:,j)=(R(:,j)-min(R(:,j)))/(max(R(:,j))-min(R(:,j)));
end

function r=Fuzzy_similarity_matrix(x,method)
% 本子程序功能是分别用7种方法建立模糊相似矩阵。
[m,n]=size(x);
for i=1:m
    for j=1:m
        ed(i,j)=sqrt(sum((x(i,:)-x(j,:)).^2));
        hd(i,j)=sum(abs(x(i,:)-x(j,:)));
    end
end
for i=1:m
    for j=1:m
        K(i,j)=sum(x(i,:).*x(j,:));
    end
end
for i=1:m
    K(i,i)=0;
end
E=max(max(ed)); H=max(max(hd)); M=max(max(K));
disp('  ');
if method==1
    disp('夹角余弦法建立模糊相似矩阵');
elseif method==2
    disp('相似系数法建立模糊相似矩阵');
elseif method==3
    disp('欧氏距离法建立模糊相似矩阵');
elseif method==4
    disp('海明距离法建立模糊相似矩阵');
elseif method==5
    disp('指数距离法建立模糊相似矩阵');
elseif method==6
    disp('数量积法建立模糊相似矩阵');
elseif method==7
    disp('算术平均最小法建立模糊相似矩阵');
end
for i=1:m
    for j=1:m
        if method==1
            r(i,j)=sum(x(i,:).*x(j,:))/sqrt(sum(x(i,:).^2)*sum(x(j,:).^2));
        elseif method==2
            r(i,j)=sum(abs((x(i,:)-mean(x(i,:))).*(x(j,:)-mean(x(j,:)))))/sqrt(sum((x(i,:)-mean(x(i,:))).^2*sum((x(j,:)-mean(x(j,:))).^2)));
        elseif method==3
            r(i,j)=1-ed(i,j)/E;
        elseif method==4
            r(i,j)=1-hd(i,j)/H;
        elseif method==5
            r(i,j)=exp(-sum(abs(x(i,:)-x(j,:))));
        elseif method==6
            if i==1
                r(i,j)=1;
            else
                r(i,j)=sum(x(i,:).*x(j,:))/M;
            end
        elseif method==7
            mink=[];
            for k=1:n
                mink=[mink min(x(i,k),x(j,k))];
            end
            r(i,j)=2*sum(mink)/sum(x(i,:)+x(j,:));
        end
    end
end

function b=FSM_square(x)
% 本子程序功能是求模糊相似矩阵的平方;输入参数x必须是一个模糊方阵,输出为x的平方。
[m,n]=size(x);
if m~=n
    disp('输入矩阵不是方阵,请重新输入一个方阵');
    return
end
for i=1:n
    for j=1:m
        for k=1:n
            b(k)=min(x(i,k),x(k,j));
        end
        xsquare(i,j)=max(b);
    end
end
b=xsquare;

function bibao=Transtive_closure(x)
% 本子程序功能是求模糊相似矩阵x的传递闭包,以用于模糊聚类。
[m,n]=size(x);
if m~=n
    disp('输入矩阵不是方阵,请重新输入一个方阵');
    return
end
xsquare=FSM_square(x);
while sum(sum(x==xsquare))~=n^2;
    x=xsquare;
    xsquare=FSM_square(x);
end
bibao=xsquare;

function A=Classification_level(x)
% 本子程序功能是求动态模糊聚类的分类水平。
[m,n]=size(x);
B=sort(reshape(x,1,m*n));
M=B(end);
i=1;
while B(i)<M
    L=find(B==B(i));
    B(L(2:end))=[];
    i=i+1;
end
L=find(B==M);
B(L(2:end))=[];
A=B;

function [cls,nmb]=Computing_cls(x)
% 本子程序功能是求存储分类的数组cls及相应分类下包含对象的数量。
m=size(x,1);
for i=1:m
    J=[];
    for j=1:m
        if x(i,j)==1
            J=[J,j];
        end
    end
    cl{i}=J;
end
for i=1:m
    for j=i+1:m
        if length(cl{j})==length(cl{i})
            if sum(cl{j}==cl{i})==length(cl{i})
                cl{i}=[];
            end
        end
    end
end
cls=[];
for i=1:m
    if ~isempty(cl{i})
        cls=[cls cl(i)];
    end
end
nmb=length(cls);
        
function F=F_statistic(r,x)
% 本子程序功能是求F统计量,以确定最佳分类数目。
[cls,nmb]=Computing_cls(r);
m=size(r,1);
ubar=mean(x);
for k=1:nmb
    s=cls{k};
    if length(s)==1
        ukbar(k,:)=x(s,:);
    else
        ukbar(k,:)=mean(x(s,:));
    end
    fz(k)=length(s)*sum((ukbar(k,:)-ubar).^2)/(nmb-1);
    for l=1:length(s)
        fm0(k,l)=sum((x(s(l))-ukbar(k,:)).^2);
    end
    fm(k)=sum(fm0(k,:)/(m-nmb));
end
F=sum(fz)/sum(fm);

function Dynamic_clustering(x,method)
% 本子程序功能是动态聚类图,返回的是原始数据的动态聚类图。
close all;
A=Standard(x);
B=Fuzzy_similarity_matrix(A,method);
R=Transtive_closure(B);
lamda=Classification_level(R);
m=length(lamda);
for i=1:m
    M(:,:,i)=R>=lamda(i);
    cls=Computing_cls(M(:,:,i));
    allcls{i}=cls;
end
rr=[];
for i=1:length(allcls)-1
    ind=[];
    for j=1:i
        ind0=[];
        for k=1:i+1
            ind0=[ind0 ismember(allcls{i+1}{k}(1),allcls{i}{j})];
        end
        ind=[ind;ind0];
    end
    tmp=allcls{i+1};
    for r=1:i
        if sum(ind(r,:))==2
            tp=ind(r,:);
            s=find(tp==1);
            allcls{i+1}{r}=tmp{s(1)};
            allcls{i+1}{r+1}=tmp{s(2)};
            t1=[1:r-1,r+1:i];
            t2=[1:r-1,r+2:i+1];
            for t=1:i-1
                allcls{i+1}{t2(t)}=allcls{i}{t1(t)};
            end
            rr=[rr r];
        end
    end
end
seq0=allcls{end};
seq=[];
for i=1:m
    seq=[seq allcls{end}{i}];
end
figure,hold on,title('动态聚类图');
wide=30; 
high=25;
x0=wide*(m-1); 
y0=high*m;
par=(wide+high)/2;
axis([-2.5*par,x0+1.5*par,-par,y0+par]);
axis('off');
lbx{1}=wide*[0:m-1];
for i=2:m
    lbx{i}(rr(m+1-i))=(lbx{i-1}(rr(m+1-i))+lbx{i-1}(rr(m+1-i)+1))/2;
    for u=1:rr(m+1-i)-1
        lbx{i}(u)=lbx{i-1}(u);
    end
    for v=rr(m+1-i)+1:m+1-i
        lbx{i}(v)=lbx{i-1}(v+1);
    end
end
for i=1:m
    text((i-1)*wide-wide/9,-0.5*high,int2str(allcls{end}{i}));
end
text(-3.45*par,(m+0.5)*high,'分类水平');
text(-3.1*par,0.5*high,num2str(lamda(m)));
text(x0+0.3*par,(m+0.5)*high,'分类数');
text(x0+1.35*par,0.5*high,num2str(m));
for i=1:m
    line([(i-1)*wide (i-1)*wide],[0 high]);
end
for i=2:m
    text(-3.1*par,(i-0.5)*high,num2str(lamda(m+1-i)));
    text(x0+1.35*par,(i-0.5)*high,num2str(m+1-i));
    line([lbx{i-1}(rr(m+1-i)) lbx{i-1}(rr(m+1-i)+1)],[i-1 i-1]*high);
    for j=1:length(lbx{i})
        line([lbx{i}(j) lbx{i}(j)],[i-1 i]*high);
    end
end