《灰色模型1.docx》由会员分享，可在线阅读，更多相关《灰色模型1.docx（19页珍藏版）》请在课桌文档上搜索。
1、灰色预测模型灰色预测是就灰色系统所做的预测.所谓灰色系统是介于白色系统和黑箱 系统之间的过渡系统，其具体的含义是：如果某一系统的全部信息已知为白色系 统，全部信息未知为黑箱系统，部分信息已知，部分信息未知，那么这一系统就 是灰箱系统.一般地说，社会系统、经济系统、生态系统都是灰色系统.灰色系统理论认为对既含有已知信息又含有未知或非确定信息的系统进行 预测，就是对在一定方位内变化的、与时间有关的灰色过程的预测.尽管过程中 所显示的现象是随机的、杂乱无章的，但毕竟是有序的、有界的，因此这一数据 集合具备潜在的规律，灰色预测就是利用这种规律建立灰色模型对灰色系统进行 预测.灰色预测模型只需要较少的观
2、测数据即可，这和时间序列分析，多元回归分 析等需要较多数据的统计模型不一样.因此，对于只有少量观测数据的项目来 说，灰色预测是一种有用的工具.一、GM(Ij)模型多色系统理论是邓聚龙教授在1981年提出来的，是一种对含有不确定因索 彳进行预测的方法.通.过鉴别系统因素之间发展趋势的相异程度,进行关联分 析，并通过对原始数据进行生成处理来寻找系统的变化规律，生成较强规律性数 据字歹U，然后建立相应微分方程模型，从而预测事物未来的发展&势和未来状态.GM(IjI：M，,.代忤差分方拧，按时间累加后所形成的新的时间序列呈现的规律可用一阶线性微 分方程的解来逼近，用生成数序列代替原始时间序列，弱化原始
3、时间序列的随机 性.这样可以对变化过程作较长时间的描述，进而建立微分方程形式的模型.其 建模的实质是建立微分方程的系数，将时间序列转化为微分方程，通过灰色微分 方程口以建立抽象系统的发展模型.经证明，经阶线性微分方程的解逼近所揭.始时间数.J .功的.1.1 GM(Ll)模型的建立成为有规律的时间序列数据，对其作一次累加生成运算，即令XSSn=l从而得到新的生成数列X，X=卜(1),*(2),数列 则生成的离散形式的微分方程具体的形式为即表示求解上述微分方程，解为x(t) = cea+-a当尸1时，x(r) = x(l), KPc = X(I)-,则可根据上述公式得到离散形式微分方程 a的具体
4、形式为x(r) = pl)-,)+其中，G项中的X为公的背景值，也称初始值；a， 是待识别的灰色参数, 出按白化导数定义有dx 1. x(t+t)-x(t)=Iimdt 3，显然，当时间密化值定义为1时，当41时，则上式可记为dx=lm(x( +，)- x(f)dt Af这表明虫是一次累减生成的，因此该式可以改写为dt= (r + l)-(r)dt当，足够小时，变量JV从Xa)到、(，+/)是不会出现突变的，所以取力，)与+，)的平均值作为当，足够小时的背景值，即/)=；/() +“(+ 1),匚1： 值带入式子，整理得?0)(r 1) = -g 力 + X (Z +1) + w由其离散形式可
5、得到如下矩阵：二x+ X f(2)2 L令y = ()(2),x (3),0n)r-扪+X一如(I)+ x S) 1 = ( P称Y为数据向量，8为数据矩阵，为参数向量.则上式可简化为线性模型：Y = Ba由最小二乘估计方法得a = pzl = (rl B上式即为GM(Ij)参数出的矩阵辨识算式,式中(3%厂/丫事实上是数据矩阵 8的广义逆矩阵.将求得的。，值代入微分方程的解式，则,r) = wi)-Vfl-0+-a a其中，上式是GM(IJ)模型的时间响应函数形式，将它离散化得xm(Z) = f0i(l)-L-0+-(0)对序列外)(。再作累减生成可进行预测.即(O = Jf(O-iO-l)
6、=卜用(IY上式便是GM(1,1)模型的预测的具体计算式.或对x。) = ce + 2求导还原得 aa1.2 GM(1,1)模型的检验GM(1,1)模型的检验包括残差检验、关联度检验、后验差检验三种形式.每种检验对应不同功能：残差检验属于算术检验，对模型值和实际值的误差 进行逐点检验；关联度检验属于几何检验范围，通过考察模型曲线与建模序列曲 线的几何相似程度进行检验，关联度越大模型越好；后验差检验属于统计检验， 对残差分布的统计特性进行检验，衡量灰色模型的精度.残差检验残差大小检验，即对模型值和实际值的残差进行逐点检验.设模拟值的残差序列为O,则e(o)(r) = x(o)(r)-x(o)(r
7、)令(/)为残差相对值，即残差百分比为) =x(o)(O-x(o)(O令A为平均残差，X = -f(t). /=I1 2设残差的方差为s；,则=一2卜(力-可.故后验差比例C为c = S2s, n /=I误差频率 P为 P= Pet)-e0.950.800.700.70C0.350.500.65表1灰色预测精确度检验等级标准一般要求V 20%,最好是V10%,符合要求.关联度检验关联度是用来定量描述各变化过程之间的差别.关联系数越大，说明预测值 和实际值越接近.设X(o)(r) = x(o)(l),x(o)(2),并5)X(o)(r) = x(o)(l),x(o)(2),.,xo)()序列关联
8、系数定义为min |x0) (r) -x(0) (r)| + max |x(0) (r) -x(0) (r)4(r) = |i(0) (r) - x(0) (r) + max |i(0) (r) - x(0) (r)1, r = 0式中，W()7()(U为第/个点X和髀的绝对误差，a)为第，个数据的关联系 数，P称为分辨率，即取定的最大差百分比，O0l,一般取0 = 0.5.X和土的关联度为精度等级关联度均方差比值小误差概率好(1级)0.900350.95合格(2级)0.800.500.80勉强(3级)0.700.650.70不合格(4级)0.65后验差检验后验差检验，即对残差分布的统计特性进
9、行检验.检验步骤如下：1、计算原始时间数列X()=U,9)(2),小)(研的均值和方差仁(U2 tl1 ro,=-yo) , s； /-12、计算残差数列*)=忖,e的均值0和方差s；e=-e(t) ,2n 1n /=I其中dQ) = ()一小)(f) , f = l,2,为残差数歹U3、计算后验差比值C = Sj S4、计算小误差频率P=P 1|e(0) (Z) - e I 0，当Co,当P时.，称模型为小残差概率合格模型.PC模型精度0.950.800.700.65勉强合格0.65不合格表3后验差检验判别参照表1.3 残差GM(1,1)模型当原始数据序列X建立的GM(1,1)模型检验不合格
10、时，可以用GM(1,1)残 差模型来修正.如果原始序列建立的GM(1,1)模型不够精确，也可以用GM(1,1) 残差模型来提高精度.若用原始序列X建立的GM(1,1)模型x(r+l) = 0l)-Kr+-a a可获得生成序列X的预测值，定义残差序列*)/)=/伏)-/(左).若取k=t, 什1,则对应的残差序列为) = e ,e22),.,e()()计算其生成序列e伏)，并据此建立相应的GM(1,1)模型e(t + 1) = ?0(1) - I + 殳ae 4得修正模型x( + l)= xo)(l)- e-ak + - + (k-t)(-ae) e(l) 4 eL。L ae其中次ZT)=F为修
11、正参数.0 kl时，模拟精度 低于70%；发展系数-a1.5时，模拟精度低于50%.进一步对预测误差进行考虑，当发展系数-av.3时，1步预测精度在98% 以上，2步和5步预测精度都在90%以上，10步预测精度亦高于80%；当发展系 数- 0.8时，1步预测精度已低于70%.通过以上分析，可得下述结论：1、当-av.3时，GM(1,1)可用于中长期预测；2、当0.3-q0.5时，GM(1,1)可用于短期预测，中长期预测慎用；3、当O.5v-aO.8时，GM(1,1)作短期预测应十分谨慎；4、当0.8(l) = x(0)(0) = 79还原求出X的模拟值.由文(O)(A) =戈(a)-*D(A-
12、l),取 = 2,3,4,得(。)=(。),S(O)(2),主(3), x(0) (4) = (79,74.281,74.3584, 76.4513)通过预测，得到实际值与预测值如下表：实际值预测值相对误差(&)第一学期79790第二学期74. 82574. 28100. 73%第三学期74. 2974. 35840. 0921%第四学期76. 9876. 45130. 7051%表4四学期的实际值与预测值的误差表因为e(Z)v10%,那就可得学生的预测值，与现实值进行比较得出该模型精度较 高，可进行预测和预报.我们对学生未来两个学期(也就是第五、六个学期)的成绩进行预测，分别 为 77. 5
13、602 分和 78. 6851 分.例：某大型企业1999年至2004年的产品销售额如下表，试建立GM(I1)预测模 型，并预测2005年的产品销售额。年份199920002001200220032004销售额 （亿元）2.673.133.253.363.563.72解：设 X伏)=2.67, 3.13, 3.25, 3.36, 3.56, 3.72)构造累加生成序列X0,)=2.67, 5.80, 9.05,12.41, 15.97,19.69构造数据矩阵B和数据向量，1111第3步BTB 二-4.235-7.425-10.73-14.19-17.83707.46375-54.41-54.4
14、1(BTB)T =0.0086670.0943190.0943191.226382d = (BB) BYn =-0.0438792.925663得出预测模型z,)dt0.043879 ”=2.925663xk +1) =69.3457 e- 66.6757(x(0)(1)=2.67; -=-66.6757)第5步残差检验(1)根据预测公式，计算划”),得Xm(k)= 2.67, 5.78, 9.03, 12.43, 15.97, 19.68, 19.69 (U ,6)(2)累减生成父(幻序列,2=1,2,6反(Q= 2.67, 3.11, 3.25, 3.40, 3.54, 3.71)原始序列
15、:Xw(k) = 2.67, 3.13, 3.25, 3.36, 3.56, 3.72)(3)计算绝对残差和相对残差序列绝对残差序列：=0, 0.02, 0, 0.04, 0.02, 0.01)相对残差序列:= 0, 0.64%, 0, 1.19%, 0.56%, 0.27%)相对残差不超过1.19%,模型精确度高。第6步进行关联度检验(1)计算序列X与i(0)的绝对残差序列(k)=0,0.02,0,0.04,0.02,0.01)min 伏) = min 0,0.02,0,0.04,0.02,0.01 = 0maxA)(Q ) = max 0,0.02,0,0.04,0.02,0.01) =
16、0.04(2)计算关联系数由于只有两个序列(即一个参考序列，一个被比较序列)故不再寻求第二级最小差和最 大差。(2) =min (&) + Pmax (左)()Pmx ()求得小k)=l, 0.5,1, 0.33, 0.5, 0.67)(3)计算关联度=-,()=067r=0.67是满足P=0.5时的检验准则r0.6的。第7步后验差检验(1)计算：x(0)=-2.67+3.13+3.25+3.36+3.56+3.72=3.28 6(2)计算序列的均方差：(3wd,2)2 =O 3671n-(3)计算残差的均值：=-(A)=0.015 6(4)计算残差的均方差：Y)-2 ,S,=( J-),z2
17、=0.0I522 H-IS(5)计算 C： C = -l =0.0152/0.3671=0.0414 S?(6)计算小残差概率：S0 =0.6745 0.3671=0.2746/=(A:) - = 0.15,0.005,0.015,0.025,0.005,0.005)所有曾都小于S。，故小残差概率P (2),为系统特征数据序列，而Xr)=(X。)(2),国。)()XT)=(X)(1), Xf)(2),XPs)X-=(X,X-(2),为相关因素序列，X,为Xf)的AGO序列(i = l,2,N), Zf)为Xj的紧邻均 值生成序列，则称4t(Q + z(Q = SbiXT (k)i=2为GM(1
18、,N)模型.在GM(1,N)模型中，。称为系统发展数据，人吧)(Q称为驱动项，儿称为驱动系数，A = ,尻岳,如了称为参数列.设引为系统特征数据系列，X：)(i = 2,3,N)为相关因素数据序列，X,为诸Xr)的LAGO序列，Zf)为Xf)的紧邻均值生成序列，则-Y)蜷播f(2)B =-z1(3)嫂XD，Y =X。)-Z；DS)EDs)即)5)*)(初则参数列a = ,b,b2,。Nr的最小二乘估计满足a = (BB)- N N砰)(& +1)= (0)-ZMD (八 e-ak +1 yto(D (& +1) a i=2a i=23、累减还原式举伏+ 1) = 染/+ 1) =茶伏+ 1)
19、-染4、GM(1,N)差分模拟式为铲出)=-对)(&) +1嫡I”) BtY设2 = ,。1,岳.,九，则称- + W=仇嫂 + 如？ + b斓at为GM(I, N)模型染) + z,(Q =仇嫂(女) +0/) +久母(Z)的白化方程，也 称影子方程.由A = 0,万，岳，,如=(BTB)T BrYf则z)N1、白化方程十知=兄尻针按差分法离散，得到解为dti=2引)=“忱 JOE(3力 + x1m(0)-次 J bix(O)dt/=2 Ji=2 JNN=e-,xfo(0) - B染(0) + X bixt)ealdtt=2i=2J2、当Xr (i = 2,3,N)变化幅度很小时，可视肝)为
20、灰常量，则GM(I,N)/=2模型染”)+材) = 琛伏)的近似响应时间式为(石(0)取为染)/=21.3 GM(I)模型与GM(1,N)模型的比较GM(1,1)是基本预测模型，具有全信息、.而GM(1,N)为分析模型、因子模型, 它不具有全信息，一般不适应于预测.然而，当有必要对多因子的系统作整体的、 全局的、动态的分析时，就需要使用GM(1,N).1.4 GM(LN)模型实例分析设系统特征数据序列为X1w = (2.874,3.278,3.307,3.390,3.679),相关因素数据序列为 Xf )=(7.04,7.645,8.075,8.53,8.774),试建立 GM(1,2)模型.
21、GM(1,2)白化方程为zr()包- + *)=M)dt对X作一次累加后的序列为Xf)=(XD(I)，X(2), Xh(3), /)(4), XD)=(2.874,6.125,9.459,12.849,16.528)X?=(引),引),(3),引)(4),(5)= (7.04,14.685,22.76,31.29,40.064)对X做紧邻均值生成序列为Z；D=(ZT (2), Zy ,z,(4), z(5) = (4.513,7.8055,11.154,14.6885) 则数据矩阵B及数据向量y为-z?引)-4.51314.685 口-3.278B =智)(3)蝴(3)-7.805522.76
22、，Y =斓3.3070(4)研-11.15431.29研3.390V)4(5),-14.688540.064穹(5)3.679对参数列,7 =力丁进行最小二乘估计，得a = (BBYBY = BtY= 0.9068则GM (1,2)模型为 Jl)+ 2.2273Xf) = 0.9068x5dt1-时间响应式为科)(4 + 1) = (Wo)(I) - 2. ( & +1 )e- + 2r? (Z +1)=(2.874 - 0.407 lx? (& + l)22273 + 0.407x? (女 +1)模拟数据，见下表序号实际数据 0)()模拟数据 (女)残差式 k) = ()(k)-S(k)相对
23、误差X(S(Q23. 2782. 7700. 50815. 5%33. 3073. 548-0. 2417. 3%43. 3905. 535-0. 1454. 3%53. 6793. 5820. 0972. 6%三、GM(2,1)模型设X(O)=W),/)(2),.,x(o)()为原始时间序列，对它进行一次累加生成 运算，得生成列为X=x,x,/)()其中”(Z) =为户，=生成的时间序列构成一灰色模块，建立 (k)x(l (k 1)k = 2,3,4,wz + l, ,nnz若第L时刻及以后的拟合值不满足要求，即*(Z) = X(S/)-”)伏)，(Z = L,L + 1,/)不在误差允许范
24、围内，对/)按照前面处理X的办法，建立残差GM(2,1)模型，可求出”)为金=白(L),/(L +1), ,)(),)( + w)对进行一次累减可得幺.残差变识可进行多次，直到满足要求为止，最后，可把( + i) +铲)( + i)作为时间序列第 + i时刻的预测值.四、灰色模型程序 GM(1,1) MATLAB 程序function =greymodelshili(y)%本程序主要用来计算根据灰色理论建立的模型的预测值。%应用的数学模型是GM(l,l)o %原始数据的处理方法是一次累加法。 y=2.67 3.13 3.25 3.36 3.56 3.72;%原始数据n=length(y);D=
25、y * 0;OneS(n-1);yy=ones(n,l);yyd)=y(l); for i=2:nyy(i)=yy(i-l)+y(i);%生成序列x(l), (1)在 X 的上方end B=ones(n-1,2); for i=k(n-l)B(i,l)=-(yy(i)+yy(il)2;B(i,2)=l;%数据矩阵end BT=B; for i=l:(n-l)z(i,l)=(yy(i)+yy(i+1)/2; %z(l), (1)在 Z 的上方 end C=ones(l,n-l)*z;E=y*0jz;F=z*z; for j=l:n-lYN=y(j+l); %数据向量 end YN=YN,;A=i
26、nv(BT*B)*BT*YN;a=A(l);%发展系数u=A(2);%灰作用量t=ua; t-test=l;%预测的个数，根据题目可以自己取i=l:t_test+n;yys(i+1 )=(y( 1 )-t).*exp(-a.*i)+t;% 白化响应式yys(l)=y(l);%生成序列的估计值，x(l)的估计值,(1)在X的上面xxs(i+ l)=yys(i+ l)-yys(i)xxs(l)=y(l)%原始数据的估计值，x(0)的估计值,(0)在X的上面for i=l:nDet(i)=abs(xxs(i)-y(i)y(i) end DET=Det*ones(n, 1 )n for j=n+t_t
27、est:-l:2ys(j)=yys(j)-yys(j-l);%生成序列估计值残差end for i=l:nerror(i)=xxs(i)-y(i); %计算残差endC=std(error)std(y)%调用统计工具箱的标准差函数计算后验差的比值Cx=l:n;xs=2:n+t_test;yn=ys(2:n+t_test);plot(x,y,r,xs,yn,*-b,);disp(预测值 为：,num2str(xxs( 1 :n+t_test)J);% 所有的估计值 GM(2,1) MATLAB 程序x=79 74.825 74.29 76.98 78.8 77.1;n=length(x);xl(
28、l)=x(l);for i=2:nxl(i)=xl(i-l)+x(i);endfor i=3:n-l;X(l)=-(xl(2)-x 1( 1)/2;X 二-(xl -xl -xl )/2;X(i)=-(x l(i+ l)-x l(i)-x l(i-l)+x l(i-2)2;endfor i=kn-l;Z(i)=-(xl(i)+xl(i+1)/2;endA=ones(n-l,l);B=IX, ZA;Y=x0(2mybeta=inv(B,*B)*B*Y;a=beta(l)b=beta(2)u=beta(3)m=input(m=);x2(0)=0;x2( 1 )=x 1 (1 );x 1 (-1 )=0; x 1(0)=0;for k=2:mx2(k)=(2*u+(2+a-b)*x2(k-1 )+a*(x2(k-2)-x l(k-3)(a+b+2);endx3(l)=x(l);for k=2:mx3(k)=x2(k)-x2(k-l);endx2x3
