灰色预测模型实例※※

灰色预测模型.所谓灰色系统是介于白色系统和黑箱是灰箱系统.一般地说，社会系统、经济系统、生态系统都是灰色系统.灰色系统理论认为对既含有已知信息又含有未知或非确定信息的系统进行.尽管过程中灰色预测模型只需要较少的观测数据即可，这和时间序列分析，多元回归分析等需要较多数据的统计模型不一样.因此，对于只有少量观测数据的项目来说，灰色预测是一种有用的工具.一、GM(1,1)模型1981年提出来的，是一种对含有不确定因素系统进行预测的方法.目前使用最广泛的灰色预测模型是关于数列预测的一个变量、一阶微分的GM(1,1)模型.GM(1,1)模型是基于灰色系统的理论思想，将离散变量连续化，用微分方程.方程可以建立抽象系统的发展模型.GM(1,1)功的.GM(1,1)模型的建立 GM(1,1) 设时间序列X有n个观察值，X0 x0x02, ,x0n为了使其成为有规律的时间序列数据，对其作一次累加生成运算，即令n1

x0nXXx,x2,,xnX一般近似地服从指数规律.则生成的离散形式的微分方程具体的形式为dxdtaxu即表示变量对于时间的一阶微分方程是连续的.求解上述微分方程，解为ux(t)cea(t1)au当t=1时，x(t)x(1)，即cx(1)a，则可根据上述公式得到离散形式微分方程的具体形式为xxueatu aa aa其中，axxdx的背景值，也称初始值；a，u是待识别的灰色参数，adt为发展系数，反映x的发展趋势；u为灰色作用量，反映数据间的变化关系.按白化导数定义有dxlimdt t

x(tt)t1时，当t1时，则上式可记为dxdtlim(x(tt)x(t))t1dx这表明 是一次累减生成的，因此该式可以改写为dtdxdtx(1)(t1)x(1)(t)当txx(tx(ttx(t与1xtt)的平均值作为当tx(1)2x(1)t)x(1)t1将值带入式子，整理得1x(0)t1)2ax(1)t)x(1)t1u1由其离散形式可得到如下矩阵： 1 2x(1)(1)x(1)(2) x(0)(2) 1 x(0)(3)a 2x(1)(2)x(1)(3) u

x(0)(n) 1 x(1)(n1)x(1)(n) 2 (0)令 Y(0)2)x(0(3, x, T( )(0)11 2x(1)(1)x(1)(2) 111 1B

2x(1)(2)x(1)(3) 1 1 x(1)(n1)x(1)(n) 1 2 auT称YB为数据矩阵，为参数向量.则上式可简化为线性模型：YB由最小二乘估计方法得

au

BT

1

BTY GM(1,1参数a,uBTB1BTY事实上是数据矩阵B的广义逆矩阵.将求得的au值代入微分方程的解式，则(1) t)(x(1)u)eat (1) a a其中，上式是GM(1,1)模型的时间响应函数形式，将它离散化得 u uˆ(1)t)x(0)(1)aea(t1)a ˆt再作累减生成可进行预测.即ˆ(0)t)ˆ(1)t)ˆ(1)t1) u x(0)(1)a

1ea

ea(t1) 上式便是GM(1,1)模型的预测的具体计算式.ux(tceat求导还原得auˆ(0)t)a(x(0)(1))ea(t1)aGM(1,1)模型的检验GM(1,1)GM(1,1)模型的检验包括残差检验、关联度检验、后验差检验三种形式.对残差分布的统计特性进行检验，衡量灰色模型的精度.残差检验残差大小检验，即对模型值和实际值的残差进行逐点检验.设模拟值的残差序列为e(0)(t)，则e(0)(t)x(0)(t)ˆ(0)(t)令(t)为残差相对值，即残差百分比为 (t)x(0)t)ˆ(0)t) 令

1n

(t).

x(0)(t) t1

1n

2设残差的方差为S2，则S2 e(t)e .故后验差比例C为CS /S，2 2 n 2 1t1e误差频率P为PPe(t) 0.6745Se1对于C,P检验指标如下表：检验指标好合格勉强不合格P>0.95>0.80>0.70<0.70C<0.35<0.50<0.65>0.65表1灰色预测精确度检验等级标准一般要求20%，最好是10%，符合要求.关联度检验和实际值越接近.ˆ 设 X(0)(t) ˆ(0)(1),ˆ(0)(2),,ˆ(0)(n) X(0)(t) x(0)(1),x(0)(2),,x(0)(n)序列关联系数定义为minˆ(0)(t)x(0)(t)maxˆ(0)(t)x(0)(t)

,t0 ˆ(0)t)x(0)t)max ˆ(0)(t)x(0)t1

,t0式,ˆ(0)(t)x(0)(t)为第t个点x(0)和ˆ(0)的绝对误差(t)为第t个数据的关联系0,一般取0.5.x(0)t)ˆ(0)t)的关联度为

r1n

t精度等级关联度均方差比值小误差概率好（1级）精度等级关联度均方差比值小误差概率好（1级）合格（2）勉强（3）不合格（40.900.800.700.700.350.500.650.650.950.800.700.70表2精度检验等级60后验差检验后验差检验，即对残差分布的统计特性进行检验.检验步骤如下：0 1

x(0)(1),x(0)(2)x(0)(n)的均值和方差

1n t

x(0)(t) , S1

1n t

x(0)(t)x 2、计算残差数列e(0)

e(0)(1),e(0)(2),,e(0)(n)

的均值e和方差s22 1n

1n

2e nt1

e(0)(t) , S2

nt1

e(0)(t)e其中e(0)t)x(0)t)ˆ(0)(t),t1,2,,n为残差数列.3、计算后验差比值CS S2 14、计算小误差频率e ePP e(0)(t)

0.6745S1令S=0.67450

(t)e(0)(te|PPS1 0C0

0，当CC0

时，称模型为方差比合格模型；若对0P00P0C模型精度>0.95>0.80>0.70<0.70<0.35<0.5<0.65>0.65优合格不合格

0P

时，称模型为小残差概率合格模型.表3后验差检验判别参照表残差GM(1,1)模型X(0)GM(1,1)GM(1,1)差模型来修正.GM(1,1)残差模型来提高精度.X(0)GM(1,1)模型u uˆ(1)t1)[x(0)(1)aeataX(1)的预测值，定义残差序列e(0)(k)x(1)(k)ˆ(1)(k).k=,t+1,n，则对应的残差序列为 e(0)(k) e(0)(1),e(0)(2), ,e(0)(n)计算其生成序列e(1)(kGM(1,1)模型得修正模型

uˆ(1)(t1)[e(0)(1)eekae

ueae u

ux(1)(t1)x(0)(1)aeaka(kt)(ae)e(0)(1)ee ae1 kt其中(kt)0 kt为修正参数.应用此模型时要考虑：1、一般不是使用全部残差数据来建立模型，而只是利用了部分残差.2(ktt的取值有关.GM(1,1)模型的适用范围定理：当GM(1,1)发展系数|a|2时，GM(1,1)模型没有意义.我们通过原始序列X0与模拟序列X0进行误差分析，随着发展系数的增i i大，模拟误差迅速增加.当发展系数a0.3时，模拟精度可以达到98%以上；发展系数a0.595%以上；发展系数a170%；发展系数a1.550%.进一步对预测误差进行考虑，当发展系数a0.3时，1以上，2590%以上，1080%；数a0.8时，170%.通过以上分析，可得下述结论：1、当a0.3时，GM(1,1)可用于中长期预测；2、当0.3a0.5时，GM(1,1)可用于短期预测，中长期预测慎用；3、当0.5a0.8时，GM(1,1)作短期预测应十分谨慎；4、当0.8a1GM(1,1)模型；5、当a1GM(1,1)模型.模型实例分析例：预测学生后两个学期的成绩.学期1成绩学期2成绩学期3成绩4某学生7974.82574.2976.98则该学生成绩时间序列如下：X(0)x(0)(1),x(0)(2),x(0)(3),x(0)(4)79,74.825,74.29,76.98X作一次累加后的数列为X(1)x(1)(1),x(1)(2),x(1)(3),x(1)(4)79,153.825,228.115,305.095X做紧邻均值生成.令Z(1)(k0.5x(1)(k0.5x(1)(k1)，得Z(1)z(1)(2),z(1)(3),z(1)(4)116.4125,151.47,150.1925则数据矩阵B及数据向量Y为z(1)(2) 116.4125 x(0)(2) 74.825

Bz(1)(3) 151.47 ，Yx(0)(3)74.29z(1)(4) 150.1925 x(0)(4) 76.98对参数列a进行最小二乘估计，得76.61 (BTB)1BTYBTY0.0144[a,u]T 即 a0.0144，u76.61则GM(1,1)模型为时间响应式为

0.0144xdt

76.61当k1

ˆ(1)(k1)5399.1389e0.01445320.1389ˆ(1)(1)ˆ(0)(1)x(0)(0)79X(0ˆ(0)(k)ˆ(1)(k)ˆ(1)(k1)，取k2,3,4，得ˆ(0)ˆ(0)(1),ˆ(0)(2),ˆ(0)(3),ˆ(0)(4)79,74.281,74.3584,76.4513通过预测，得到实际值与预测值如下表：实际值预测值相对误差k第一学期79790第二学期74.82574.28100.73%第三学期74.2974.35840.0921%第四学期76.9876.45130.7051%表4四学期的实际值与预测值的误差表因为10%高，可进行预测和预报.我们对学生未来两个学期（也就是第五、六个学期）的成绩进行预测，分别为77.5602分和78.6851分.例：某大型企业1999年至2004年的产品销售额如下表，试建立GM(1,1)预测模型，并预测2005年的产品销售额。1999200020012002200320042.673.133.253.363.563.72年份销售额(亿元)X(0)(k={2.67，3.13，3.25，3.36，3.56，3.72}年份销售额(亿元)第1步 构造累加生成序列X(1)(k)={2.67，5.80，9.05，12.41，15.97，19.69}2

B和数据向量Yn 1 2x(1)(1)x(1)(2) 1

4.235 x(1)(2)x(1)(3) 2 1

7.425 1B2x(1)(3)x(1)(4) 1

1， 1

14.19 2 x(1)(4)x(1)(5) 2

17.83 1 2x(1)(5)x(1)(6) x(0)(2) 3.13 x(0)(3) 3.25Y x(0)(4)3.36n x(0)(5)x(0))a

3.563.723

计算=(BTB)1BTYb b BTB=707.46375 54.41 5 0.008667 0.094319(BTB)1=0.094319 1.226382 0.043879=(BTB)1BTY =n

2.9256634

得出预测模型dx(1)dt

0.043879x(1)=2.925663ˆ(1)(k1)=69.3457e0.043879k66.6757b（x(0)(1)=2.67；=－66.6757）a第5步 残差检验(1)(k)，得(1)(k)k=0,1,…,6)(2)(0)(kk=1,2,6(0)(k)＝｛2.67，3.11，3.25，3.40，3.54，3.71｝原始序列：X(0)(k)＝｛2.67，3.13，3.25，3.36，3.56，3.72｝(3)计算绝对残差和相对残差序列绝对残差序列：(0)＝｛0，0.02，0，0.04，0.02，0.01｝相对残差序列：＝｛0，0.64%，0，1.19%，0.56%，0.27%｝相对残差不超过1.19%，模型精确度高。第6步 进行关联度检验x(0)ˆ(0)的绝对残差序列(0))(0)＝｛0,0.02,0,0.04,0.02,0.01｝min{(0)(k)}=min｛0,0.02,0,0.04,0.02,0.01｝=0max{(0)(k)}=max｛0,0.02,0,0.04,0.02,0.01｝=0.04计算关联系数（大差。(k)min{(k)}Pmax{(k)}(k1,...,6,P0.5)(k)Pmax{(k)}求得(k)＝{1,0.5,1,0.33,0.5,0.67}计算关联度r1ni n k

(k)＝0.67r=0.67是满足P=0.5时的检验准则r>0.6的第7步 后验差检验 1(1)x(0)=[2.67+3.13+3.25+3.36+3.56+3.72]=3.286X(0)序列的均方差：[x(0)(k)x(0)]2S=( 2=0.36711计算残差的均值：6

n1[(k)]=0.015计算残差的均方差：[(k)]2S =( )1/2=0.01522 n1S(5)CC1=0.0152/0.3671=0.0414SS2(6)S=0.67450.3671=0.27460e (k)k

{0.15,0.005,0.015,0.025,0.005,0.005}所有e都小于S，故小残差概率P{e S }=1，而同时C=0.0414<0.35，故模型i 0 i 0x(1)(k1)=69.3457e0.043879k66.6757合格.第8步 预测：k=7，x(0)(8)=x(1)(8)x(1)(7)=4.23即2005年的产品销售额预测值为4.23亿元.灾变预测例：某地区平均降水量（单位：毫米）的原始数据为：Xx1,x2,...,x24={386.6,514.6,434.1,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(毫米)为旱灾年，试作旱灾预测。首先作灾变映射。x390(毫米)为异常值,则有Xx[q(1)],x[q(2)], ,x[q(6)]386.6,254.5,384.5,242.5,374.7,387.1x9,xxxx作异常值x[q(i)]到出现灾变点q(i)的映射Q0:x[q(i)]q(i)，得灾变日期序列Q0为

Q0q(1),q(2),q(3),q(4),q(5),q(6)111,2据此对Q0建立灾变日期序列的GM(1,1)模型。对Q0作一次累加生成，得Q(1){q(1)(1),q(1)(2),q(1)(3),q(1)(4),q(1)(5),q(1)(6)} 。 T

1 0.188422

1

a,

BTB

BTY9.54872

。记

的紧邻生成序 列为ZGM(1,1)为q(k0.188422z(1)(k)9.54872GM(1,1)序号响应式为b ˆ(1)(k1)(q(1) )eak b a a51.6772e0.188422k50.6773从而ˆ(k1)ˆ(1)(k1)ˆ(1(k)8.87478e0.188422k由此可得Q0的模拟序列ˆ(0)ˆ(k),k2,3,4,5,6}10.7,12.9,15.6,18.8,22.7}由(0)(k)x(0)(k)ˆ(0)(k)，得绝对残差序列及相对残差序列

(0){(0)(k),k2,3,4,5,6}}，(0)(i),i2,...,6}{0.19,0.14,0.025,0.044,0.013} i i q(i) 平均相对残差0.10，故可用

165 i2

0.08进行预测.

ˆ(k1)8.87478e0.188422kˆ(61)27，ˆ(7)ˆ(6)5即从最近一次旱灾发生的时间算起，五年之后可能发生旱灾.二、GM(1,N)模型GM(1,N)模型的建立GM(1,N)来建模分析.设X(0)(X(0)(1),X(0)(2),,X(0)(n))为系统特征数据序列，而1 1 1X(0)

1(X(0)(1),X(0)(2),,X(0)(n))2 2 2 2X(0)(X(0)(1),X(0)(2),,X(0)(n))3 3 3 3 X(0)

(X(0)(1),X(0)(2),,X(0)(n))N N N NX(1)X(0)的1-AGO序列(i1,2,NZ(1)X(1)的紧邻均i值生成序列，则称

i 1 1x(0)(k)az(1)(k)Nbx(1)(k)GM(1,N)模型.

1 1 iii2GM(1,N)模型中，a称为系统发展数据，bx(1)(k

称为驱动ii iNˆa,,b2,,bT称为参数列.N设X(0)为系统特征数据系列，X(0)(i2,3,,N)为相关因素数据序列，X(1)为1 i iX(0)1-AGOZ(1)X(1)的紧邻均值生成序列，则i 1 1z(1)(2) x(1)

x(1)(2) x(0)(2) 1

N 1 z(1)(3) x(1)

x(1)(3)

x(0)(3)B 1 2

N ，Y1 z(1)(n) x(1)(n) x(1)(n) x(0)(n)1 2 N 11则参数列[ab1

2,,

]T的最小二乘估计满足(BTB)1BTY1设b1

, ,b2

T，则称dx(1) ax(1)bx(1)bx(1)

x(1)dt 1

22 33 N NGM(1,N)x(0)(kaz(1)(kbx(1)(kbx(1)(kb

x(1)(k)的白化方程，也1称影子方程.

1 22 33 N N1 2 由ˆa,b,b ,,b1 2

(BTB)1BTY，则dx(1)1、白化方程 ax(1)

bx(1)按差分法离散，得到解为dt 1

i1i2N Nx(1)(t)eat[ bx(1)(t)eatdtx(1)(0) bx(1)(0)dt]1 ii 1 iiNi2 iNeat[x(1)(0)tNbx (0) bx (tedt](1) (1) at1 ii iit2 i2N2、当X(1)(i2,3,,N)变化幅度很小时，可视 bx(1)(k)为灰常量，则GM(1,N)i i1i2x(0)(k)az(1)(k)Nbx(1)(k)x(1)(0)x(0)(1))1 1 ii 1 1i2ˆ(1)(k1)(x(1)(0)1

bx(1)(k1))e

1

bx(1)(k1)13、累减还原式

1 a ii2

a iii2ˆ(0)(k1)a(1)ˆ(1)(k1)ˆ(1)(k1)ˆ(1)(k)1 1 1 14、GM(1,N)差分模拟式为ˆ(0)(k)az(1)(k)Nbˆ(1)(k)1 1 ii2GM(1,1)模型与GM(1,N)模型的比较GM(1,1)GM(1,NGM(1,N).GM(1,N)模型实例分析设系统特征数据序列为X(0)(2.874,3.278,3.307,3.390,3.679) ,相关因素数1X(0)7.04,7.645,8.075,8.53,8.774GM(1,2)模型.2GM(1,2)白化方程为dx1

axbx1 1 1dt 1 2X作一次累加后的序列为 X(1)

x(1)(1),x(1)(2),x(1)(3),x(1)(4),x(1)(5)

(2.874,6.125,9.459,12.849,16.528)1 1 1

1 1 1X(1)

x(1)(1),x(1)(2),x(1)(3),x(1)(4),x(1)(5)

(7.04,14.685,22.76,31.29,40.064)2 2 2 2 2 2X做紧邻均值生成序列为 Z(1)

z(1)(2),z(1)(3),z(1)(4),z(1)

(4.513,7.8055,11.154,14.6885)1 1 1 1 1则数据矩阵B及数据向量Y为z1(2)

2 4.513

14.685

2 3.278 1 2

2 z1(3)

22

7.8055

22.76

22

3.307 B

，Y

z1(4)

4

11.154

31.29

4

3.390 1

14.6885

40.064

3.679 z1(5)1

2

5 2对参数列a进行最小二乘估计，得2GM(1,2)模型为

2.2273 (BTB)1BTYBTY dx1

2.2273x0.9068x时间响应式为

1 1 1dt 1 2ˆ(k1)(x0(1)bxk)eakbxk1 1 a 2 a 22.8740.407xke2.2273k0.407xk2模拟数据，见下表实际数据模拟数据

2((k)序号 x(0)(k

ˆ(0)(k

(k)x(0)(k)ˆ(0)(k) k

x(0)(k)23.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(0)

x(0)(1),x(0)(2),,x(0)(n)

为原始时间序列，对它进行一次累加生成运算，得生成列为 X(1)

x(1)(1),x(1)(2),,x(1)(n)x(1)(ki1

x0(i)k1,2,n，生成的时间序列构成一灰色模块，建立GM(2,1)模型

d2dx(1) a bx(1)u

(1)dt2 dt若a0，令a1/abbauu/a，则上式变为d2x(1) dx(1)a bx(1)

(2)dt2 dt上两个式子从微分方程角度看，没有本质区别，但从拟合角度看，(1)式是拟合d2

dx(1) dx(1) 项，而(2)式是拟合 项， 实际上就是原始时间序列X(0)的近似，dt2 dt dtd2

dx(0) dx(0) 近似 ，而 的变化较大，拟合效果不太好，故从拟合角度看dt2 dt dt式比(1)式好。GM(2,1)模型比GM(1,1)模型多考虑了X(0)的影响，因此预测效果更好一些.按灰色系统理论，(2)式的离散形式为 a b x(0)

k 2x(1)(k)x(1)(k1)x(1)(k2)x(1)

k

2x

k x(1)(k1u (3)x(1)(0)x(1)(1)0, k2,3,4,,n数据矩阵B和数据向量Y为1

1 2

x(1)(2)x(1)(1)x(1)(0)x(1)(1)

x(1)(2)x(1)(1) 1111 x(1)(3)x(1)(2)x(1)(1)x(1)(0) x(1)(3)x(1)(2) 1B 2 1

2 1 2

x(1)(n)x(1)(n1)x(1)(n2)x(1)(n3)

x(1)(n)x(1)(n1)2

1 Y[x(0)(2),x(0)(3),,x(0)(n)]T，[a,b,u]T令k1,2,n，则(3)式可表示为矩阵形式Y上式的最小二乘解为

(BTB)1BTY当N1,NM(3)式得预测公式ˆ(1)ku2abˆ(1)(k1)aˆ(1)(k2)x(1)k3ab2ˆ(0)(0)ˆ(1)(k)ˆ(1)(k1) k2,3,4,,,n1,nmL时刻及以后的拟合值不满足要求，即(0)(k)x(0)(k)ˆ(0)(k)，LL1,n0XGM(2,1)模型，可求出ˆˆ

ˆ1

(L),ˆ

(L1),ˆ

(n),,ˆ

(nm)对ˆ进行一次累减可得ˆ0.残差変识可进行多次，直到满足要求为止，最后，ˆ(0)(ni)ˆ(0)(ni)作为时间序列第ni时刻的预测值.四、灰色模型程序GM(1,1) MATLAB程序function[]=greymodelshili(y)%本程序主要用来计算根据灰色理论建立的模型的预测值。%应用的数学模型是GM(1,1)。%原