第8章 灰色预测方法

1、 1982年我国学者邓聚龙先生创立了灰色系统理论，目前许多国家及国际组织的知名学者从事灰色系统的理论和应用研究工作。 灰色系统研究的是“部分信息明确，部分信息未知”的“小样本，贫信息”不确定性系统，它通过对已知“部分” 信息的生成去开发了解、认识现实世界。着重研究“外延明确，内涵不明确”的对象。第八章 灰色预测方法8.1 灰色系统基本原理与灰数灰色系统基本原理与灰数1、差异信息原理：、差异信息原理： 差异即信息，凡信息必有差异。2、解的非唯一性原理：、解的非唯一性原理：信息不完全、不确定的解是非唯一的。该原理是灰色系统理论解决实际问题所遵循的基本法则。3、最少信息原理：、最少信息原理：灰色系统

2、理论的特点是充分利用已占有的“最少信息”。4、认知根据原理：、认知根据原理：信息是认知的根据。5、新信息优先原理：、新信息优先原理：新信息对认知的作用大于老信息。6、灰性不灭原理：、灰性不灭原理： “信息不完全”是绝对的。一、原理一、原理二、灰数及其运算二、灰数及其运算 1、灰数：只知道大概范围而不知道其确切的数，、灰数：只知道大概范围而不知道其确切的数，通常记为：通常记为：“ ”。例如：例如： （1）多少层的楼房算高楼，中高楼，低楼。多少层的楼房算高楼，中高楼，低楼。 （2 2）多么大的苹果算大苹果，小苹果）多么大的苹果算大苹果，小苹果。2、灰数的种类、灰数的种类（1）仅有下界的灰数。）仅有

3、下界的灰数。 有下界无上界的灰数记为：有下界无上界的灰数记为： a, 、 (a)（2）仅有上界的灰数。）仅有上界的灰数。 有上界无下界的灰数记为：有上界无下界的灰数记为： - ,a a（3）区间灰数）区间灰数 既有上界又有下界的灰数：既有上界又有下界的灰数： a, a （4）连续灰数与离散灰数）连续灰数与离散灰数 在某一区间内取有限个值的灰数为离散灰数，在某一区间内取有限个值的灰数为离散灰数，取值连续地取满整个区间的灰数为连续灰数。取值连续地取满整个区间的灰数为连续灰数。（5）黑数与白数黑数与白数 当 （- ， ）或 (1, 2)，（即当 的上界、下界皆为无穷或上、下界都是灰数时，称为黑数黑数

4、， 当 a,a且a=a,时，称 为白数。为白数。（6）本征灰数与非本征灰数 本征灰数是指不能或暂时还不能找到一个白数作为其“代表”的灰数； 非本征灰数是凭借某种手段，可以找到一个白数作为其“代表”的灰数。 则称此白数为相应灰数的白化值，记为 并用 (a) 表示以a为白化值的灰数。(100)100如：托人代买一件价格为100元左右的衣服，可将100作为预测衣服价格 (100)的白化数，记为a（7）信息型灰数 因暂时缺乏信息而不能肯定其取值的数。但到一定的时间，通过信息补充，灰数可以完全变白。从本质上看，灰数可分为信息型、概念型和层次型灰数。（8）概念型灰数，也称意愿型灰数 指由人们的某种概念、意

5、愿形成的灰数。（9）层次型灰数 指由层次的改变形成的灰数。(宏观白,微观灰)a3、区间灰数的运算设灰数1 a, b, 2 c,d (ab,cd) (1) 1 + 2 a+c,b+d (2) -1 -a, -b (3) 1 - 2 =1 +(- 2) a-d,b-c (4) 1 2 minac,ad,bc,bd,maxac,ad,bc,bd (5) 1/ 2 mina/c,a/d,b/c,b/d,maxa/c,a/d,b/c,b/a (6)若k为正实数， 则： k1 ka, kb（2）对一般的区间灰数 ，将白化取值为4.灰数白化与灰度（1）有一类灰数是在某个基本值附近变动的，这类灰数白化比较容易

6、，可将其基本值为主要白化值。可记为其中 为忧动灰元。此灰数的白化值为(1)ab ()aaa0,1a()aa,a b 定义：形如 的白化称为等权白化等权白化。(1)ab 0,1定义：在等权白化中 而得到的白化值称为等权均等权均值白化。值白化。21在区间灰数取值的分布信息缺乏时，常采用等权均值白化。在灰数的分布信息已知时，常采用非等权均值白化。如：如：某人2000年的年龄可能是40岁到60岁，根据了解，此人受初中级教育12年，且20世纪60年代中期考入大学，故此人的年龄到2000年为58左右的可能性较大。或者在56岁到60岁的可能性较大。 4 0 , 6 0 注：白权化函数被用来描述一个灰数对其取

7、值范围内不同数值的“偏爱”程度。定义：设区间灰数1 a, b, 2 c,d (ab,cd)当 时，称 1与2取数一致；当 时，称1与2取数不一致。12(1)0 ,1,(1)0 ,1abab定义：起点，终点确定的左升、右降连续函数称为典型的白化权函数。f(x)10 x1x2x3x4L(x)R(x)xa定理1：区间灰数不能相消、相约。即：灰数自差一般不能等于0，仅当减数与被减数的取数一致时，灰数的自差才等于0。如： 2，5， =0 取数一致 -3，3 取数不一致 =1 取数一致 2/5，5/2 取数不一致 如： / 灰度：是灰数的测度。 灰度在一定程度上反映了人们对灰色系统之行为特征的未知程度。它

8、与相应定义信息域的长度及其基本值有关。8.2 灰灰 色色 预预 测测 概概 念念 一、灰色预测的概念 （1）灰色系统、白色系统和黑色系统 白色系统是指一个系统的内部特征是完全已知的，即系统的信息是完全充分的。 黑色系统是指一个系统的内部信息对外界来说是一无所知的，只能通过它与外界的联系来加以观测研究。 灰色系统内的一部分信息是已知的，另一 部分信息是未知 的，系统内各因素间有不确定的关系。 用灰色数学来处理不确定量，使之量化。 （2）灰色系统特点 充分利用已知信息寻求系统的运动规律。 关键：如何使灰色系统白化、模型化、优化 灰色系统视不确定量为灰色量，提出了灰色系统建模的具体数学方法，它能用时

9、间序列来确定微分方程的参数。灰色系统理论能处理贫信息系统。 （只要求较短的观测资料即可） 灰色预测法是一种对含有不确定因素的系统进行预测的方法。 灰色预测是对既含有已知信息又含有不确定信息的系统进行预则，就是对在一定范围内变化的、与时间有关的灰色过程进行预测。 （3）灰色预测法 灰色预测通过鉴别系统因素之间发展趋势的相异程度，即进行关联分析，并对原始数据进行生成处理来寻找系统变动的规律，生成有较强规律性的数据序列，然后建立相应的微分方程模型，从而预测事物未来发展趋势的状况。 灰色预测法用等时距观测到的反映预测对 象特征的一系列数量值构造灰色预测模型， 预测未来某一时刻的特征量，或达到某一 特征

10、量的时间。 （4）灰色预测的四种常见类型 灰色时间序列预测 即用观察到的反映预测对象特征的时间序列来构造灰色预测模型，预测未来某一时刻的特征量，或达到某一特征量的时间。 畸变预测 即通过灰色模型预测异常值出现的时刻，预测异常值 什么时候出现在特定时区内。 系统预测 通过对系统行为特征指标建立一组相互关联的灰色预测模型，预测系统中众多变量间的相互协调关系的变化。拓扑预测 将原始数据做曲线，在曲线上按定值寻找该定值发生的所有时点，并以该定值为框架构成时点数列，然后建立模型预测该定值所发生的时点。四种方法共同点： （1）允许少数据预测； （2）允许对灰因果律事件进行预测，如：1）灰因白果律事件：粮食

11、预测，影响因素很多，是灰因；然而粮食产量是具体，是白果。2）白因灰果律事件：项目开发预测，投入是具体，为白因；而收益暂时不清楚，为灰果。3）具有可检验性：含建模可行性的级比检验（事前检验），建模精度检验（模型检验），预测的滚动检验（预测检验）。 二、生成列 为了弱化原始时间序列的随机性，在建立灰色预测模型之前，需先对原始时间序列进行数据处理，经过数据处理后的时间序列即称为生成列。 灰色系统常用的数据处理方式有累加和累减两种。 （1）数据处理方式累加是将原始序列通过累加得到生成列。累加的规则: 将原始序列的第一个数据作为生成列的第一个数据，将原始序列的第二个数据加到原始序列的第一个数据上，其和作

12、为生成列的第二个数据，将原始序列的第三个数据加到生成列的第二个数据上，其和作为生成列的第三个数据，按此规则进行下去，便可得到生成列。记原始时间序列为： nXXXXX00000,.3,2,1生成列为： nXXXXX11111,.3,2,1(1)(0)(1)(0)(0)(1)(0)1(1)()(1)(0)1(1)()(1)(0)1(1)(1)(2)(1)(2)(1)(2)(1)( )(2)(1)( )( )(1)( )niniXXXXXXXXnXiXnXnXnXiXnXn 对非负数据，累加次数越多则随机性弱化 越多，累加次数足够大后，可认为时间序 列已由随机序列变为非随机序列。 一般随机序列的多次

13、累加序列，大多可用 指数曲线逼近。同理，可作m次累加： kimmiXkX11累减 将原始序列前后两个数据相减得到累减生成列 累减是累加的逆运算，累减可将累加生成 列 还原为非生成列，在建模中获得增量信息。一次累减的公式为： 1001kXkXkX(0)(0)(0)(0)(0)(1),(2),(3)(18)(1962,2089,1325,995,1768,1982,2658,2653,3000,2669,2577,2768,2643,2938,3650,2597,1687,1678)XXXXX(1)X =(1962,4051,5376,6371,8139,10121,12779,15432,184

14、32,21101,23678,26446,29089,32027,35677,38274,39961,41639)例 原始数据为三、关联度 关联度分析是分析系统中各因素关联程度的方法，在计算关联度之前需先计算关联系数。（1）关联系数设 nXXXkX0000,.,2,1 nXXXkX0000,.,2,1则关联系数定义为： kXkXkXkXkXkXkXkXk00000000maxmaxmaxmaxminmin)(式中： kXkX00 kXkX00minmin kXkX00maxmax为第k个点 称为分辨率，00.950.800.700.70 C0.350.500.650.65 好 合格 勉强合格

15、不合格例例 某矿某年某矿某年3-7月份的轻伤事故情况如表所示月份的轻伤事故情况如表所示:原始数据列为：原始数据列为： 34,33,31,29,26)()0(ix累加生成数列为累加生成数列为： 153,119,86,55,26)()1(ix T34,33,31,29nx,.,3x,2xyT000N月份月份34567轻伤轻伤人次人次2629313334 表表1 1 轻伤事故人次轻伤事故人次113615.10215.7015.401153119211119862118655211552621B34333129113615 .10215 .7015 .4011111365 .1025 .705 .40

16、113615 .10215 .7015 .40yBBBa 1NT1T1275 .1136645 .3495 .34975.3561211275 .11366754257848. 10172161127. 00172161127. 060001970370. 051038016914.2790531754295. 0所以所以 0.0532a 27.1038u 705.509au26)1 ()0(x其中其中称为发称为发展灰数；展灰数；称为称为内生控内生控制灰数。制灰数。所以，所以， aueau1x1kx ak1175.509705.509260532.0ke75.509705.5350532. 0

17、ke即事故预测公式为：即事故预测公式为： 75.509705.5351kx 0532. 01ke 1)1()1()0(kxkxkx生成数列的预测值、原始数列的还原值分别如表所示。生成数列的预测值、原始数列的还原值分别如表所示。 为了得到原始数列的预测值，需要将生成数列的预测为了得到原始数列的预测值，需要将生成数列的预测值作累减还原为原始值，即根据下式求得：值作累减还原为原始值，即根据下式求得：表表2 2 生成数列的预测值与误差检验生成数列的预测值与误差检验1)1(kx1)1(kx1kq-0.03153.0315340.30118.701193-0.1486.14862-0.2755.27551

18、026260k表表3 3 原始数列的还原值与误差检验原始数列的还原值与误差检验 kx0 kx0 kq12626022929.27-0.2733130.870.1343332.560.4453434.33-0.33平均值平均值30.630.606-0.006k数据方差和残差方差分别为：数据方差和残差方差分别为： 24. 86 .30346 .30336 .30316 .30296 .3026512222221S078424. 0006. 033. 0006. 044. 0006. 013. 0006. 027. 0006. 00512222222S后验差比值为：后验差比值为：0976.024.8

19、078424.012SSC小误差频率小误差频率 9362. 1006. 06745. 01kqPSqkqPP 006. 0006. 00006. 01q 264. 0006. 027. 0006. 02q 136. 0006. 013. 0006. 03q 446. 0006. 044. 0006. 04q 324. 0006. 033. 0006. 05q所以所以 1P根据根据 和和 的评价标准（表的评价标准（表3 3），本例题），本例题的预测结果的评价等级为的预测结果的评价等级为“好好”。35. 0C95. 0P可对可对8月的轻伤事故进行预测月的轻伤事故进行预测 75.509705.535

20、1kx 0532.01ke 20.18975.509705.53515x 1kx 50532.011e 17.3603.15320.18951x 61x 60 x 即根据预测，如果不能采取更有效的事故预防措施的话，即根据预测，如果不能采取更有效的事故预防措施的话，下一月份的轻伤事故人次将是下一月份的轻伤事故人次将是3636人。人。采用采用三、三、GM(1,1)模型应用实例的模型应用实例的MATLAB实现实现解解 (0)X(1)(1)累加生成数列为累加生成数列为： 年份199920002001200220032004销售额2.673.133.253.363.563.72建立建立GMGM（1 1，

21、1 1）预测模型，并预测）预测模型，并预测20052005年产品销售额年产品销售额 原始数据列为：原始数据列为： 2.67,3.13,3.25,3.36,3.56,3.722.67,3.13,3.25,3.36,3.56,3.72 (1)X 2.6700 5.8000 9.0500 12.4100 15.9700 19.6900X0=2.67,3.13,3.25,3.36,3.56,3.72;X1(1)=X0(1)for k=2:6 X1(k)=X1(k-1)+X0(k)end(2)(2)构造数据矩阵构造数据矩阵B B和数据向量和数据向量Y Y： (1)(1)(1)1( )( )(1)2ZkX

22、kXkz =z = 0 4.2350 7.4250 10.7300 14.1900 17.8300load hslitifor k=2:6 z(k)=(1/2)*(X1(k)+X1(k-1)endB = -4.2350 1.0000 -7.4250 1.0000 -10.7300 1.0000 -14.1900 1.0000 -17.8300 1.0000Y = 3.1300 3.2500 3.3600 3.5600 3.7200B=(-z(2:6) ones(5,1)Y=(X0(2:6)(3)(3)计算系数计算系数 1()TTB BB Y alfa alfa = = -0.0440 -0.0

23、440 2.9256 2.9256alpha=inv(B*B)*B*Y(4)(4)得出预测模型得出预测模型 11d0.0442.9256dXXt 1011atXkXeaa0.04469.345766.6757teu=alpha(2)/alpha(1)v=X0(1)-uu = -66.5503v = 69.2203(5)(5)进行参差检验进行参差检验 (1)X得得 1011akXkXeaa0.04469.345766.6757keu=alpha(2)/alpha(1)v=X0(1)-ufor n=0:6 X2(n+1)=v*exp(-alpha(1)*n)+uendX2u = -66.5503v

24、 = 69.22031 1）根据预测公式，计算）根据预测公式，计算 X2 = 2.6700 5.7809 9.0315 12.4283 15.9777 19.6867 23.5623(0)X得得X3(1)=X2(1)for m=1:6 X3(m+1)=X2(m+1)-X2(m)end2 2）累减生成序列）累减生成序列 X3 = 2.6700 3.1109 3.2507 3.3968 3.5494 3.7089 3.8756而原始数据为而原始数据为(0)X2.67,3.13,3.25,3.36,3.56,3.722.67,3.13,3.25,3.36,3.56,3.72 3 3）计算绝对参差和相

25、对参差序列）计算绝对参差和相对参差序列 绝对参差序列绝对参差序列 daita0 =0.0000 0.0191 0.0007 0.0368 0.0106 0.0111daita0=abs(X0-X3(1:6)(0)0,0.0191,0.0007,0.0368,0.0106,0.0111相对参差序列相对参差序列 kesi = 0.0000 0.0061 0.0002 0.0109 0.0030 0.0030kesi=daita0./X0平均相对参差平均相对参差 meankesi=mean(kesi)meankesi = 0.00390.60.5的检验准则的检验准则meanaita=mean(ait

26、a) = 0.6745（7 7）进行后验差检验）进行后验差检验 1 1）计算）计算X0X0均值、均方差均值、均方差X0mean=mean(X0) =0.2817X0mean=mean(X0) =0.2817X0std=std(X0) =0.3671X0std=std(X0) =0.3671daita0mean=mean(daita0)= 0.0130daita0mean=mean(daita0)= 0.0130daita0std=std(daita0)= 0.0137daita0std=std(daita0)= 0.0137C = 0.0372C = 0.03724 4）计算小参差概率）计算小

27、参差概率010.6745SS2 2）计算参差均值、均方差）计算参差均值、均方差3 3）计算）计算C=daita0std/X0stdC=daita0std/X0stdS0=0.6745*X0stdS0 = 0.2476|( )|kk e = 0.0130 0.0061 0.0124 0.0237 0.0025 0.0020e=abs(daita0-daita0mean)对所有的对所有的e都小于都小于S0，故小参差概率，故小参差概率0()10.95kPS P=length(find(eS0)/length(e)C = 0.03720.35,C = 0.03720.35,故预测模型是合格的。故预测模

28、型是合格的。而同时而同时(8)(8)预测预测 (1)X得得 1011atXkXeaa0.04469.345766.6757ke即即2005年的产品销售额预测值为年的产品销售额预测值为4.0498亿元。亿元。u = -66.5503v = 69.2203X2006= 4.0498X2005=X3(7)X2(8)=v*exp(-alpha(1)*7)+uX3(8)=X2(8)-X2(7)X2006=X3(8)即即2005年的产品销售额预测值为年的产品销售额预测值为3.8756亿元。亿元。即即2006年的产品销售额预测值为年的产品销售额预测值为4.0498亿元。亿元。四、四、GM(1,1)GM(1,

29、1)参差模型参差模型(1)X可获得生成序列的预测值可获得生成序列的预测值 1011akXkXeaa若用原始序列建立的若用原始序列建立的GM(1,1)GM(1,1)模型模型 对于参差序列对于参差序列(0)(1)(1)( )( )( )jXjXj(1,2, )jn若存在若存在 ，使得当，使得当 时，时， 的符号一致，且的符号一致，且 0k0kk(0)( )k04nk则称参差序列则称参差序列(0)(0)(0)(0)00(),(1),( )Ekkn为可建模参差尾部。为可建模参差尾部。的累加生成序列的累加生成序列 010()01akkkkeaa其其GM(1,1)时间响应式为时间响应式为计算参差序列计算参

30、差序列(0)(0)(0)(0)00(),(1),( )Ekkn得修正模型得修正模型(1)E 1011akXkXeaa 00()00()akkakkkeaa其中其中0001()0kkkkkk正负号取值与参差正负号取值与参差尾部符号一致尾部符号一致 0(0)(0)(0)1111(1),(2),.,( )XXXXn如果考虑得系统由若干个相互影响的因素组成：为系统特征数据序列，而其相关因素序列有 1Z 为 的紧邻均值生成序列 1X五、五、GM(1,N)GM(1,N)模型模型设设 0(0)(0)(0)2222(1),(2),.,( )XXXXn 0(0)(0)(0)(1),(2),.,( )NNNNXX

31、XXn 0iX 1iX为为的累加生成列，的累加生成列， 1(0)(1)112( )( )( )NiiiXkaZkb Xk则GM（1，N）的灰微分方程模型为： 2 ,TNa bb1TTB BB Y可利用最小二乘法求解。解得：其中(1)(1)(1)12(1)(1)(1)12(1)(1)(1)12(2)(2)(2)(3)(3)(3)( )()()NNNZXXZXXBZnXNXN(0)1(0)1(0)1(2)(3)( )XXYXn 1111ddXaXt称微分方程： 为灰色微分方程 的白化方程，也称影子方程。（1）白化方程的解（响应序列形式）为(1)2Niiib X 1(0)(1)112( )( )(

32、)NiiiXkaZkb Xk 1011111XkXa(1)2(1)Nakiiib Xke(1)21(1)Niiib Xka （2）累减还原值 011Xk 111Xk 11Xk （1）建立模型常用数据有：科学试验数据、经验数据、生产数据、决策数据；（2）序列生成数据是建立灰色模型的基础数据； （3）一般非负序列累加生成后，得到准光滑序列，对于满足光滑条件的，即可建立GM模型； （4）模型精度可以通过不同的灰数生成方式，数据的取舍，序列的调整、修正以及不同级别的参差GM模型补充得到提高；建模的思路：建模的思路： （5）灰色系统理论采用参差检验、关联度检验、后验差检验三种方法检验，判断模型的精度。

33、8.4 灾变预测(1), (2),., ( )Xxxx n 灰色灾变预测的任务是给出下一个或几个异常值出现的时刻，以便人们提前预防，采取政策，减少损失。定义：设原始序列为给定上限异常值则称满足 ( )x q i的序列 (1), (2),., ( )Xx qx qx q m为上灾变序列。且称且称(0)(1), (2),., ( )Qqqq m为灾变日期序列。为灾变日期序列。命题：设(0)(1), (2),., ( )Qqqq m且灾变日期序列且灾变日期序列GM(1,1)序号响应式为：序号响应式为：其累加序列为(1)(1)(1)(1)(1),(2),.,( )Qqqqm的紧邻生成序列为(1)Q(1

34、)Z则称 1( )( )q kazk为灾变GM（1，1）模型。 111akqkqeaa1q k 11qk 1 qk即 01akqea1q k 0(1)1a kqea 0(1)1aakeqea(1), (2),., ( )Xxxx n定义：设原始序列为n为现在，给定异常值(0)(1), (2),., ( )Qqqq m相应的灾变日期序列相应的灾变日期序列其中其中( )()q mn为最近一次灾变发生的日期，为最近一次灾变发生的日期，则称则称1q m为下一次灾变发生的预测日期，为下一次灾变发生的预测日期， q mk为未来第为未来第k次灾变的预测日期。次灾变的预测日期。 例 某地区平均降水量(单位：m

35、m)的原始数据为X=386.6,514.6,434.0,484.1,647.0,399.7, 498.7,701.6,254.5,463.0,745.0,398.3, 554.5,471.1,384.5,242.5,671.7,374.7, 458.9,511.3,530.8,586.0,387.1,454.4 规定年降水量为 灾害年，试作灾害预测。390()mm解解按照按照( )390()x tmm为异常值，则有为异常值，则有 (1), (2),., (6)Xx qx qx q=386.6，254.5 ，384.5，242.5，374.7，387.1 (1), (9), (15), (16)

36、, (18), (23)xxxxxxx=386.6,514.6,434.0,484.1,647.0,399.7 . 498.7,701.6,254.5,463.0,745.0,398.3 . 554.5,471.1,384.5,242.5,671.7,374.7 . 458.9,511.3,530.8,586.0,387.1,454.4; for n=1:24 if x(n)=390 n, X=x(n) end endn = 1X = 386.6000n = 9X = 254.5000n = 15X = 384.5000n = 16X = 242.5000n = 18X = 374.7000n = 23X = 387.1000(0)(1), (2),., (6)Qqqq得灾变日期序列：得灾变日期序列：1,9,15,16,18,23作一次累加生成作一次累加生成(1)(1)(1)(1)(1),(2),.,(6)Qqqq1,10,25,41,59,82求得参数向量求得参数向量1TTB BB Y0.18849.5487于是得