模糊聚类分析例子

1、1. 模糊聚类分析模型环境区域的污染情况由污染物在4个要素中的含量超标程度来衡量。设这 5个环境区域的污染数据为 X1 =(80, 10, 6, 2),X2=(5O, 1,6, 4),X3=(90, 6, 4, 6),X4=(4O, 5, 7, 3),X5=(1O, 1,2, 4).试用模糊传递闭包法对X进行分类。8050由题设知特性指标矩阵为：X*9040101016516647224634数据规格化：最大规格化x' M其中:M jmax(X1j,X2j,., Xnj)0.8910.860.330.560.10.860.6710.60.5710.440.510.50.110.10.2

2、90.67Xo构造模糊相似矩阵：采用最大最小法来构造模糊相似矩阵R (rij)5 5,10.540.620.630.240.5410.550.700.53R 0.620.5510.560.370.630.700.5610.380.240.530.370.381利用平方自合成方法求传递闭包t(R)依次计算R2, R4,R8,由于R8 R4，所以t(R) R4R2R410.630.620.630.530.6310.560.700.530.620.5610.620.53 ,0.630.700.6210.530.530.530.530.53110.630.620.630.530.6310.620.70

3、0.530.620.6210.620.53 = R80.630.700.6210.530.530.530.530.531选取适当的置信水平值0,1,按 截矩阵进行动态聚类。把t(R)中的元素从大到小的顺序编排如下：1>0.70>0.63>062>053. 依次取 =1,0.70,0.63, 062, 053，得10t(R)100001000001000 00 00 0，此时X被分为 5类：x X2 , X3 , X4 , X51 00 1t(R)0.7100000 01 00 11 00 001010000，此时X被分为4类：X1, X2 , X4, X3, X501t

4、 (R) 0.6311010110100010011010000，此时 X被分为 3类：x1 , X2, X4, x3 , x501t(R) 0.6211110111101111011110000，此时 X被分为 2类：X1 , X2 , X4,X3, X5011t(R)0.5311，此时 X被分为 1 类：X1,X2,X3,X4,X51Matlab程序如下:%数据规格化MATLAB程序a=80 10 6 250 1 6 490 6 4 640 5 7 310 1 2 4;mu=max(a) for i=1:5for j=1:4r(i,j)=a(i,j)/mu(j);end end%采用最大最

5、小法构造相似矩阵r=0.88891.00000.85710.33330.55560.10000.85710.66671.00000.60000.57141.00000.44440.50001.00000.50000.11110.10000.28570.6667;b=r： for i=1:5for j=1:5R(i,j)=sum(mi n( r(i,:);b(:,j)')/sum(max(r(i,:);b(:,j)');endend%利用平方自合成方法求传递闭包t(R)0.62990.69850.55991.00000.3818矩阵合成的MATLAB函数fun cti onrha

6、t=hech(ij;n=len gth(r);for i=1: n for j=1: n rhat(i,j)=max(mi n(r(i,:);r(:,j)');end end求模糊等价矩阵和聚类的程序R= 1.00000.54090.62060.62990.24320.54091.00000.54780.69850.53390.62060.54781.00000.55990.36690.24320.53390.36690.38181.0000;R1=hech (R)R2=hech (R1)R3=hech (R2)bh=zeros(5);bh(fi nd(R2>0.7)=12. 模

7、糊综合评判模型某烟草公司对某部门员工进行的年终评定，关于考核的具体操作过程,以对一名员工的考核为例。如下表所示，根据该部门工作人员的工作性质,将18个指标分成工作绩效（U1 ）、工作态度（U2）、工作能力（U3）和学习成长（U4 ）这4各子因素集。员工考核指标体系及考核表一级指标二级指标评价优秀良好一般较差差工作绩效工作量0.80.150.500工作效率0.20.60.10.10工作质量0.50.40.100计划性0.10.30.50.050.05工作态度责任感0.30.50.150.050团队精神0.20.20.40.10.1学习态度0.40.40.10.10工作主动0.10.30.30.2

8、0.1性360度满意0.10.20.50.20.1度工作能力创新能力0.10.30.50.20自我管理0.20.30.30.10.1能力沟通能力0.20.30.350.150协调能力0.10.30.40.10.1执行能力0.10.40.30.10.1学习成长勤情评价0.30.40.20.10技能提咼0.10.40.30.10.1培训参与0.20.30.40.10工作提供0.40.30.20.10请专家设定指标权重，一级指标权重为:A040.3,0.2,0.1二级指标权重为:A 0.2,0.3,0.3,0.2A0.3,020.1,0.2,0.2A30.1,020.3,0.2,0.2A 0.3,0

9、.2,0.2,0.3对各个子因素集进行一级模糊综合评判得到:Bl0.39,0.39,0.26,0.04,0.01B2A2 G R20.21,0.37,0.235,0.125,0.06B3AiR30.15,0.32,0.355,0.125,0.06B4aJr40.27,0.35,0.24,0.1,0.021/7 1/513这样，二级综合评判为:0.390.390.260.040.01A0.210.370.2350.1250.060.150.320.3550.1250.060.270.350.240.10.2B air0.4,0.3,020.10.28,0.37,0.27,0.09,0.04根据最

10、大隶属度原则，认为该员工的评价为良好。同理可对该部门其他员工 进行考核。3. 层次分析模型你已经去过几家主要的摩托车商店，基本确定将从三种车型中选购一种，你选择的标准主要有：价格、耗油量大小、舒适程度和外观美观情况。经反复思考比较，构造了它们之间的成对比较判断矩阵。1131515A= 1718三种车型（记为a,b,c ）关于价格、耗油量、舒适程度和外表美观情况的成 对比较判断矩阵为a 1b 1/2c 1/3211/ 21/511/71/271舒适程度外表11/ 31/ 5311/ 4151/ 31/511/ 7根据上述矩阵可以看出四项标准在你心目中的比重是不同的，请按由重到轻顺序将它们排出。解

11、：用matlab求解层次总排序的结果如下表准则价格耗油量舒适程度外表总排序权值准则层权值0.58200.27860.08990.0495方案层a0.53960.10560.62670.18840.4091单排序b0.29700.74450.27970.73060.4416权值c0.16340.14990.09360.08100.1493Matlab程序如下:clc,clear n1=4;n2=3;a=11/311/8 1/5 1/3 1;b1=11/21/31/2 1；b2=11/51/21/7；b3=11/31/51/4 1 ;b4=11/51/31/71;x,y=eig(eval(char

12、( 'b' ,int2str(i);% 一致性指标RIri=0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45;x,y=eig(a);%x为特征向量，y为特征值 lamda=max(diag(y);n um=fi nd(diag(y)=lamda);w0=x(: ,n um)/sum(x(:, num);w0%准则层特征向量%准则层一致性比例CR0=(lamda-n 1)/( n1-1)/ri( n1) for i=1: n1Iamda=max(diag(y);%方案层的特征向量%方案层的一致性比例n um=fi nd(diag(y)=lamda);w

13、1(:,i)=x(:, num)/sum(x(: ,n um);CR1(i)=(Iamda-n2)/( n2-1)/ri( n2);end w1CR1, ts=w1*w0, CR=CR1*wO%ts为总排序的权值，CR为层次总排序的随机一致性比例 %当CR小于0.1时，认为总层次排序结果具有较满意的一致性并接受该结果，否则对判断矩阵适当修改4. 灰色预测GM ( 1,1 )模型某地区年平均降雨量数据如表某地区年平均降雨量数据年份122345678降雨量390.6412320559.2380.8542.4553310561年份10911121314151617降雨量300390.66325404

14、06.2313.8576587.6318.5规定hz=320，并认为x®(i)v=hz为旱灾。预测下一次旱灾发生的时间解:初始序列如下X® =(390.6,412,320,559.2,380.8,542.4,553,310,561,300,632,540,406.2,313.8,576,587.6,318.5)由于满足x(0)(i)<=320的x(0)(i)为异常值，易得下限灾变数列为x0z= (320,310,300,313.8,318.5)其对应的时刻数列为t = (3,8,10,14,17)建立GM (1,1 )模型对原始数据t做一次累加，即 t(1) = (3

15、,11,21,35,52)构造数据矩阵及数据向量计算a，b a=-0.2536 ，b=6.2585 建立模型 y=-24.6774+27.6774*ex p(.253610*t)模型检验年份原始值模型值残差相对误差级比偏差333.000887.98960.01040.00130.5161101010.2960-0.29600.0296-0.0324141413.26810.73190.05230.0783171717.0983-0.09830.0058-0.0627通过计算可以预测到第六个数据是22.0340由于22.034与17相差5.034，这表明下一次旱灾将发生在五年以后。计算的MATL

16、AB程序如下:clc,clear a=390.6,412,320,559.2,380.8,542.4,553,310,561,300,632,540,406.2,313.8,576,587.6,318.5'； x0=fi nd(av=320);x0=x0： n=len gth(x0) lamda=x0(1: n-1)./x0(2:n) ran ge=minm ax(lamda) x1=cumsum(x0) for i=2:n z(i)=0.5*(x1(i)+x1(i-1);endB=-z(2: n)', o nes( n-1,1);Y=x0(2: n)'u=BY x=d

17、solve( 'Dx+a*x=b' ,'x(0)=x0');x=subs(x, 'a' ,'b' ,'x0' ,u(1),u(2),x1(1);yuce1=subs(x, 't' ,0:n-1);digits(6),y=v pa(x) yuce=x0(1),diff(yuce1) ep sil on= x0-yuce rho=1-(1-0.5*u(1)/(1+0.5*u(1)*lamda yuce1=subs(x, 't' ,0:n);yuce=x0(1),diff(yuce1) 5

18、. Verhulst预测模型在实际问题中，常遇到原始数据本身呈 S形的过程，这时，可取原始数据为X 其一次累减生成（1 IAGO ）为x（0），建立Verhulst模型，直接对x进行预测（模拟）。现以中国道路交通事故死亡人数为例，建立交通事故死亡人数Verhualst预测模型。由中国交通年鉴、中国汽车工业年鉴等可得近年 来中国道路交通事故死亡人数统计资料，见表 14 0表14道路交通事故死亡人数统计年份19901991199219931994199519961997死亡人数（万人）4.935.335.876.356.637.157.377.39年份19981999200020012002200

19、32004死亡人数（万人）7.818.359.3910.5910.9410.4410.71解：19902003年中国道路交通事故死亡人数曲线见图 2，可见曲线呈S形，故可建立Verhulst模型进行预测，其建模过程如下。(1 )设X为19902003年死亡人数的原始数据序列，即(1) (1) (1) (1) (1)x (X ,x2 ,X3 .X14 ) (4.93, 5.33, 5.87, 6.35, 6.63, 7.15,7.37, 7.39, 7.81, 8.35, 9.39,10.59,10.94,10.44)(0)X(2 )对x(1)作一次累减生成(1 IAGO ),由(x!0) ,x

20、20) ,.Xi(4)(4.93, 0.4, 0.54, 0.48, 0.28, 0.52, 0.22,0.02, 0.42, 0.54,1.04,1.2, 0.35, -0.5)对X作紧邻均值生成，令zk10.5(xk1)xJdk2,3,.14(zj, z31),.z14)(5.13, 5.6, 6.11,6.49, 6.89, 7.26, 7.38,7.6, 8.08, 8.87, 9.99,10.765,10.69)对参数列(a,b)T进行最小二乘估计，得T 1 T(B B) BY0.1280.0089z1(5) Verhulst 模型为dx( 0.128Xdt0.0089( X)2(6

21、)模型精度检验(过程略)平均相对误差 = 3.74%，则模型精度为二级；同时算得绝对关联度 g为 0.9845，均方差比值C为0.2355，则模型精度为一级，可见模型精度较高，可用于事故预测。Matlab编程如下clc,clear x1=4.93 5.33 5.87 6.35 6.63 7.15 7.37 7.39 7.81 8.35 9.39 10.59 10.94 10.44;n=len gth(x1);nia n=1990:2003;plot(nian,x1,'o-');x0=diff(x1);x0=x1(1),x0 for i=2:n z1(i)=0.5*(x1(i)+

22、x1(i-1);endB=-z1(2:e nd)', z1(2:e nd)'A2Y=x0(2:e nd)' abhat=BYx=dsolve( 'Dx+a*x=b*x2','x(0)=x0' ); x=subs(x, 'a' ,'b' ,'x0' ,abhat(1),abhat(2),x1(1);yuce=subs(x, 't' ,0:n-1);digits® x=v pa(x) x1_all=x1;ep sil on= x1_all-yuce delta=abs(

23、e psil on. /x1_all) delta_mea n=mea n( delta) x1_all_0=x1_all-x1_all(1);yuce_0=yuce-yuce(1);s0=abs(sum(x1_all_0(1:e nd-1)+0.5*x1_all_0(e nd);s1=abs(sum(yuce_0(1:e nd-1)+0.5*yuce_0(e nd);tt=yuce_0-x1_all_0;s1_s0=abs(sum(tt(1:e nd-1)+0.5*tt(e nd);absdegree=(1+s0+s1)/(1+s0+s1+s1_s0) c=std(e psilo n,1)/

24、std(x1_all,1) yuce=subs(x, 't' ,0:n)6.GM (2,1 )模型19962001年上海市上网户数数据序列为x(0)(x(0)(1),x(0)(2),.x(0)(6)(0.33,0.9,10.24, 42.24,88.24,104.1)在互联网发展初期，增长势头十分强劲。因此，定理5引入的一阶缓冲算子弱化该序列的增长趋势，一阶缓冲序列仍记为x(0)，x(0)=(41,49,61,78,96,104 )，试以该序列为基础建立GM(2,1)模型解:x(0)的1-AGO序列和1-IAGO序列分别为x(1) =(41,90,151,229, 325, 4

25、29)(1)x(0) =(0, 8,12,17,18, 8)x(1)的紧邻均值生成序列z(1) =(0, 65.5,120.5,190, 277, 377)B=x(0) (2),x(0) (3),z(2),1z(3),1，Y= x(0) (2) x(0) (3)x(0 )(2)x(0 )(3)x(0) (1) x(0 )(2)x(0)(6),z(1)( 6),1 x(0) (6)x(0)x(0)a1 a2 b(BtB) 1BtY1.09220.195931.7983可得,GM(2,1 )时间响应式(1)x (k1) =-.532426*ex p(.865973*t)+203.850*ex p(

26、 .226223*t)-162.317所以预测的数据为(41, 51,63, 77, 92,104) 误差分析实际数据模拟数据残差相对误差4951.0148-2.01480.0411ep silo n= xO-xO_hat%计算残差6163.1412-2.14120.03517877.21110.78890.01019692.15483.84520.0401104104.4780-0.47800.0046Matlab程序如下clc,clearx0=41,49,61,78,96,104;n=len gth(x0);x1=cumsum(x0) %x1为累加数列a_x0=diff(x0);a_x0=

27、0,a_x0% a x0为累减数列for i=2:nz(i)=0.5*(x1(i)+x1(i-1);endB=-x0(2:e nd)',-z(2:e nd)', o nes( n-1,1);Y=a_x0(2:e nd)'u=BY %a1 , a2 , b 的值x=dsolve( 'D2x+a1*Dx+a2*x=b','x(0)=c1,x(5)=c2');x=subs(x, 'a1' ,'a2' ,'b' ,'c1' ,'c2' ,u(1),u(2),u(3),

28、x1(1),x1 (6);yuce=subs(x, 't' ,0:n-1);digits(6),x=v pa(x)%x为时间响应式delta=abs(e psilo n./xO)%计算相对误差7.波形预测模型上海证券交易所综合指数的波形预测。根据上海证券交易所综合指数的周收盘指数数据,从1997年2月21日到1998年10月31日的周收盘指数曲线如图所示:解:取9条等间隔的等高线，分别为1=1140，2=1170，3=1200，4=1230，5=1260，6=1290，7=1320，8 =1350，9=1380i的等高时刻序列分别为对应于1 =1140，(8.8,12.8,60

29、.2,71.8,72.7,73.6)(4.4,31.7,34.2,41,42.4,76.8,78.3)对应于2=1170 ，Q20) =(5.2,19.8,23,25.6,26.9,31.2,34.8,39.5,44.6,76,76.2,79.2对应于 3=1200，4=1230，5=1260，6=1290，7=1320，8=1350 ，9=1380分别为q30) =(5.9,19.5,24.8,25.2,26.5,30.3,46.2,53.4,55.4,75.5,79.7(6.5,19.2,28.3,29.5,49.7,50.8,56.2,76.4,82.9,85q50) =(7,14.2,

30、16.4,16.5,18.8,56.7,75.2)q60) =(8.3,13.4,16.9,56.2,74.6)Q80)=(9.6,12.5,61.8,69.8,70.9,71.8)q90) =( 10.8,12.4,64.1,69)对 Qi(0)(i=1,2,3,9)序列，进行GM (1,1 )预测，起响应时分别为y1=-109.738+114.138*ex p(0.214831*t) y2=-94.0581+99.2581*ex p(0.158430*t) y3=-94.6529+100.553*ex p(.166865*t) y4=-145.162+151.662*ex p( .159938*t) y5=-3.68695+10.6869*ex p( .446077*t) y6=-12.1881+20.4881*ex p(.550388*t) y7=-176.276+185.076*ex p( .191636*t) y8=-182.496+192.096*ex p( .185059*t) y9=-34.6983+45.4983*ex p(.488018*t)对在1998年11月到2000年3月这5个月进行预测，可得等高时刻的预测序列(0)Q1(99.3,123.1 )(0)Q2(97.3,114.1,133.6(0)Q3