灰色预测法GM(1,1)理论及应用

灰色预测法GM(1,1)理论及应用
一、概念
1. 灰色预测法是一种对含有不确定因素的系统进行预测的方法。

灰色系统是介
于白色系统和黑色系统之间的一种系统。

灰色系统内的一部分信息是已知的，另一部分信息时未知的，系统内各因素间具有不确定的关系。

2. 灰色预测，是指对系统行为特征值的发展变化进行的预测，对既含有已知信
息又含有不确定信息的系统进行的预测，也就是对在一定范围内变化的、与时间序列有关的灰过程进行预测。

尽管灰过程中所显示的现象是随机的、杂乱无章的，但毕竟是有序的、有界的，因此可以通过对原始数据进行生成处理来寻找系统变动的规律，生成有较强规律性的数据序列，然后建立相应的微分方程模型，从而预测事物未来发展趋势的状况。

灰色预测是利用这种规律建立灰色模型对灰色系统进行预测。

二、灰色预测的类型
1. 灰色时间序列预测；即用观察到的反映预测对象特征的时间序列来构造灰色
预测模型，预测未来某一时刻的特征量，或达到某一特征量的时间。

2. 畸变预测；即通过灰色模型预测异常值出现的时刻，预测异常值什么时候出
现在特定时区内。

3. 系统预测；通过对系统行为特征指标建立一组相互关联的灰色预测模型，预
测系统中众多变量间的相互协调关系的变化。

4. 拓扑预测；将原始数据作曲线，在曲线上按定值寻找该定值发生的所有时点，
并以该定值为框架构成时点数列，然后建立模型预测该定值所发生的时点 三、GM （1，1）模型的建立 1. 数据处理
为了弱化原始时间序列的随机性，在建立灰色预测模型之前，需先对原始时间序列进行数据处理，经过数据处理后的时间序列即称为生成列。

i. 设()()()()()()()()(){}
,,, (00000)
123X X X X X n = 是所要预测的某项指标的原始
数据，计算数列的级比()()()
(),,,,()
00123X t t t n X t λ-=
=。

如果绝大部分的级比
都落在可容覆盖区间(,)221
1
n n e
e
-++内，则可以建立GM(1,1)模型且可以进行灰色
预测。

否则，对数据做适当的预处理。

方法目前主要有数据开n 方、数据取对数、数据平滑。

预处理的数据平滑设计为三点平滑，具体可以按照下式处理
()()()()()()
()()/00001214X t X t X t X t ⎡⎤=-+++⎣⎦
()()()()
()()/00013124X X X ⎡⎤=+⎣⎦ ()()()()()()/000134X n X n X n ⎡⎤=-+⎣⎦
ii. 预处理后对数据作一次累加生成处理，即：将原始序列的第一个数据作为生
成列的第一个数据，将原始序列的第二个数据加到原始序列的第一个数据上，其和作为生成列的第二个数据。

按此规则进行下去，便可得到生成列。

根据()
()()()101
k
n X
k X n ==∑，得到一个新的数列
()()()()()()()()(){}
,,,...11111123X X X X X n =
这个新的数列与原始数列相比，其随机性程度大大弱化，平稳性大大增加。

2. 新数列的变化趋势近似地用下面的微分方程描述。

()
()11dX aX u dt
+= 其中：a 称为发展灰数；u 称为内生控制灰数。

3. 模型求解。

令()()()[(),(),,()]00023T n Y X X X n =⋯，ˆα
为待估参数向量，ˆa u α⎛⎫
= ⎪⎝⎭
， ()()()()()()(()())(()())(()())111111112 12123 12
11 12X X X X B X n X n ⎡⎤-+⎢⎥⎢⎥⎢⎥
-+⎢⎥=⎢⎥⋯⋯⎢⎥⎢⎥--+⎢⎥⎣⎦
， 于是模型可表示为
ˆn Y B α= 通过最小二乘法得到：
()ˆ1
T T n B B B Y α
-= 求解微分方程，即可得灰色预测的离散时间响应函数：
()()()()ˆ1011at u u X t X e a a -⎡⎤+=-+⎢⎥⎣⎦
,,,...,0121t n =- ()()ˆ11X
t +为所得的累加的预测值，将预测值还原即为： ()()()ˆˆˆ()()-()01111X
t X t X t +=+ 注：若数据经过预处理，则还需经过相应变换才能得到实际预测值。

4、模型检验
灰色预测检验一般有残差检验、关联度检验和后验差检验。

1) 残差检验
()()()ˆˆˆ()()-(-)0111X
t X t X t = ()()()()()()ˆ000t X
t X t ∆=- ()()
()(),,,,()
0012t t t n X t ε∆==
分别求出预测值、绝对误差值和相对误差值，计算出平均相对误差判断精度是否理想。

2) 关联度检验
i. 定义关联系数()t η
()()()()()
()()
()
min max ()max 0000t t t t t ρη∆+∆=
∆
+∆
其中：①()()0t ∆为第t 个点()0X 与()ˆ0X
的绝对误差； ②ρ称为分辨率，0<ρ<1,一般取ρ=0.5；
③对单位不一，初值不同的序列，在计算相关系数前应首先进行初始
化，即将该序列所有数据分别除以第一个数据。

ii. 定义关联度()1
1n
t r t n η==∑，称为()()0X t 与()()ˆ0X
t 的关联度 根据上述方法算出()()ˆ0X
k 与原始序列()()0X k 的关联系数，然后计算出关联度，根据经验，当ρ=0.5时，关联度大于0.6便满足检验标准。

3) 后验差检验
计算原始序列标准差和绝对误差序列的标准差分别为:
1S =
,2
S =计算方差比21
S C S =
，小误差概率()()(){}.00
106745P P t S =∆-∆<，令
()()()00t e t =∆-∆，.0106745S S =，则{}0t P P e S =<
检验指标P 和C 与灰色预测精度检验等级标准如下表所示： XXX 表 C <0.35
<0.5
<0.65
≥0.65
四、残差模型修正
若用原始经济时间序列()0X 建立的GM （1，1）模型检验不合格或精度不理想时，要对建立的GM （1，1）模型进行残差修正或提高模型的预测精度。

修正的方法是建立GM （1，1）的残差模型。

设))(),...,2(),1(()0()0()0()0(n εεεε=其中，
()()()0k x k ε=-()ˆ()1x k 为)1(X 的残差序列。

若存在k 0，满足
1.的符号一致；)(,)0(0k k k ε≥∀
2.40≥-k n ，则称|))(||,...,)1(||,)((|)0(0)0(0)0(n k k εεε+为可建模残差尾段，仍记为))(),...,1(),(()0(0)0(0)0()0(n k k εεεε+=
设))(),...,1(),(()0(0)0(0)0()0(n k k εεεε+=为可建模残差尾段，其一次累加序列
))(),...,1(),(()1(0)1(0)1()1(n k k εεεε+=的GM(1,1)模型的时间响应式为
0)]([0)0()1(,))(()1(ˆ0
k k a b
e a b k k k k a ≥+-=+--ε
εεεεεε
则残差尾段的模拟序列为
))(ˆ),...,1(ˆ),(ˆ(ˆ)0(0)0(0)0()0(n k k εεεε
+= 其中
0)]
([0)0()0(,))()(()1(ˆ0k k e a b k a k k k a ≥--=+--εε
εεεε
➢ 若用)0(ˆε
修正)1(ˆX 则称修正后的时间响应式 ⎪⎪⎩
⎪⎪⎨⎧≥-±+-<+-=+----0)]([0)
0()0(0
)
0()
1(,))(())1((,))1(()1(ˆ0k k e
a b k a a b e a b x k k a b e a b x k x k k a ak ak εεεεε 为残差修正GM(1,1)模型，简称残差GM(1,1)。

其中残差修正值)]([0)0()0(0))()(()1(ˆk k a e a b
k a k ----=+εε
εεεε
的符号应与残差尾段)0(ε的符号保持一致。

➢ 若)1()0()1()1()0())1()(1()1(ˆ)(ˆ)(ˆ----=--=k a a e a
b
x e k x k x k x
则相应的残差修正时间响应式
⎪⎪⎩
⎪⎪⎨⎧≥-±--<--=+----0)]([0)
0()0(0
)
0()0(,))(())1()(1(,))1()(1()1(ˆ0k k e
a b k a e a b x e k k e a b x e k x k k a ak a ak a εεεεε 称为累减还原式的残差修正模型。

取定k 后，按此模型，可对k>k0的模拟值进行休整，修正后的精度如下表：
就只有考虑采用其它模型或对原始数据序列进行适当取舍。

再用P 和C 检验预测效果。

五、GM(1,1)模型的适用范围
灰色GM(1,1)模型评价推广 ( 1) 灰色GM(1,1)模型优点
灰色GM(1,1)预测模型在计算过程中主要以矩阵为主, 它与MATLAB 的结合解决了它在计算中的问题. 由MATLAB 编制的灰色预测程序简单实用, 容易操作, 预测精度较高.
( 2) 灰色GM(1,1)模型的缺点
该模型是指运用曲线拟合和灰色系统理论对我国人口发展进行预测的方法,
因此它对历史数据有很强的依赖性, 而且GM (1,1)的模型没有考虑各个因素之间的联系. 因此, 误差偏大, 尤其是对中长期预测, 例如对中国人口总数变化情况做长期预测时, 误差偏大, 脱离实际. 下面我们来讨论GM(1,1)模型的适用范围.
GM(1,1)模型的白化微分方程：
(1)
(1)dX aX u dt
+= 其中a 为发展系数,
可以证明，当GM(1,1)的发展系数||2a ≥时，GM(1,1)模型无意义。

因此，
(,][,)22-∞-⋃+∞是GM(1,1)发展系数a 的禁区。

在此区间，GM(1,1)模型失去意
义。

一般地，当||2a <时，GM(1,1)模型有意义。

但是，随着a 的不同取值，预测效果也不同。

通过数值分析，有如下结论：
（1）当.03a -≤时，GM(1,1)的1步预测精度在98%以上，2步和5步预测精度都在97%以上，可用于中长期预测；
（2）当..0305a <-≤时，GM(1,1)的1步和2步预测精度都在90%以上，10步预测精度也高于80%，可用于短期预测，中长期预测慎用； （3）当..0508a <-≤时，GM(1,1)用作短期预测应十分慎重；
（4）当.081a <-≤时，GM(1,1)的1步预测精度已低于70%，应采用残差修正模型；
（5）当1a ->时，不宜采用GM(1,1)模型。

如果要考虑到多因素的联系和影响, 此时我们不妨建立GM( 1, n) 模型. GM( 1, N) 模型能模拟系统发展的动态过程, 不但吸收了传统的灰色模型的建立, 而且建立了多中改进的灰色模型, 提高了预测精度.
论文小结处：与传统的数理统计模型相比，该模型在…预测方面具有明显优点：① 无需典型的概率分布；② 减少时间序列的随机性；③ 小样本即可计算；④ 计算简便。

用灰色理论预测…理论可靠，方法较简单。

对原始数据系列长度要求不高，即使在系列较短的情况下也能取得令人满意的预报结果，弥补了其他方法无法进行短期系列观测资料的…的预测。

本文建立的模型经拟合精度检验（P= ，C=），模型判为…，预测精度高，能达到预测要求。

例1（数列预测）：设原始序列
)679.3,390.3,337.3,278.3,874.2())5(),4(),3(),2(),1(()0()0()0()0()0()0(==x x x x x X
试用GM(1,1)模型对)0(X 进行模拟和预测，并计算模拟精度。

解：第一步：对)0(X 进行一次累加，得
)558.16,897.12,489.9,152.6,874.2()1(=X
第二步：对)0(X 作准光滑性检验。

由
)
1()
()()1()0(-=k x k x k ρ
得5.029.0)5(,5.036.0)4(,54.0)3(<≈<≈≈ρρρ。

当k>3时准光滑条件满足。

第三步：检验)1(X 是否具有准指数规律。

由
)(1)
1()()()1()1()
1(k k x k x k ρσ+=-=
得29.1)5(,36.1)4(,54.1)3()1()1()1(≈≈≈σσσ
当k>3时，5.0],5.1,1[)()1(<=∈δσk ，准指数规律满足，故可对)1(X 建立GM(1,1)模型。

第四步：对)1(X 作紧邻均值生成，得
)718.14,184.11,820.7,513.4()1(=Z
于是
⎥
⎥
⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡----=679.3390.3337.3278.3)5()4()3()2(,1718.141184.111820.71513.41)5(1)4(1)3(1)2()0()0()0()0()1()
1()1()1(x x x x Y z z z z B 第五步：对参数列T b a ],[ˆ=α
进行最小二乘估计。

得 ⎥
⎦
⎤
⎢⎣⎡-==-0653.30372.0)(ˆ1Y B B B T T α 第六步：确定模型
0653.30372.0)1()
1(=-x dt
dx 及时间响应序列
402151.82276151.85))1(()1(ˆ0372.0)0()1(-=+-=+-k ak e a
b
e a b x k x
第七步：求)1(X 的模拟值
)5558.16,9422.12,4605.9,1060.6,8704.2())5(ˆ),4(ˆ),3(ˆ),2(ˆ),1(ˆ(ˆ)1()1()1()1()1()1(==x x x x x X
第八步：还原求出)0(X 的模拟值。

由
)(ˆ)1(ˆ)1(ˆ)1(ˆ)1()1()1()1()0(k x k x k x a k x
-+=+=+ 得
)
6136.3,4817.3,3545.3,2320.3,8740.2())5(ˆ),4(ˆ),3(ˆ),2(ˆ),1(ˆ(ˆ)0()0()0()0()0()0(==x x x x x X
第九步：检验误差。

由下表可算出残差平方和：
第十步：预测)1(+k x
......
8928.3)7(ˆ1991.24402151.82276151.85)7(ˆ7505.3)6(ˆ3063.20402151.82276151.85)6(ˆ)0(60372.0)1()0(50372.0)1(==-===-=⨯⨯x e x
x
e x
例2 （灾变预测）：某企业生产用原料属受自然灾害影响较大的农产品。

一般来说，自然灾害的发生有其偶然性，但对历史数据的整理，仍可发现一定的规律性。

为尽量减少生产不受自然灾害的影响，该企业希望了解影响原料供应的规律性并提前做好原料储备，所收集数据见下表，并规定每亩平均收获量小于320千克时为欠收年份，将影响原料的正常供应，现应用灰色灾变预测来预测下次发生欠收的年份。

形成初始序列)0(ω。

本例初始序列：)17,14,10,8,3()0(=ω
一次累加生成序列：)52,35,21,11,3()1(=ω
)1(ω的紧邻均值生成序列：)5.43,28,16,7()1(=Z
第二步：按)1(Z 建GM(1,1)模型。

⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣
⎡=⎥
⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡----=1714108)5()4()3()2(,15.43128116171)5(1)4(1)3(1)2()0()0()0()0()1()
1()1()1(ωωωωY z z z z B ⎥
⎦
⎤
⎢⎣⎡-==-25834.625361.0)(ˆ1Y B B B a T T
67702.2467702.27])([)1(25361.0)0()1(-=+
-=+-t at e a
b
e a b t t ωω 第三步：预测
当t=6时，684.73)6(ˆ)1(=ω
6848.21)6(ˆ)0(=ω
因此，下次发生收获量小于320千克的年份为：2011年至2012年，即四至五年后将出现欠收年份。

其他预测类型见参考书。

五、残差GM(1,1)模型
当GM(1,1)模型精度不符合要求时，可使用残差序列建立GM(1,1)模型，对原来模型进行修正，以提高精度。

定义4 设
))(),...,2(),1(()0()0()0()0(n εεεε=
其中，)()()0(k x k =ε-)(ˆ)1(k x
为)1(X 的残差序列。

若存在k 0，满足 1.的符号一致；)(,)0(0k k k ε≥∀
2.40≥-k n ，则称
|))(||,...,)1(||,)((|)0(0)0(0)0(n k k εεε+ 为可建模残差尾段，仍记为
))(),...,1(),(()0(0)0(0)0()0(n k k εεεε+=
命题1 设))(),...,1(),(()0(0)0(0)0()0(n k k εεεε+=为可建模残差尾段，其一次累加序列
))(),...,1(),(()1(0)1(0)1()1(n k k εεεε+=
的GM(1,1)模型的时间响应式为
0)]([0)0()1(,))(()1(ˆ0k k a b
e a b k k k k a ≥+-=+--ε
εεεεεε
则残差尾段的模拟序列为
))(ˆ),...,1(ˆ),(ˆ(ˆ)0(0)0(0)0()0(n k k εεεε+= 其中
0)]
([0)0()0(,))()(()1(ˆ0k k e a b k a k k k a ≥--=+--εε
εεεε
定义5 若用)0(ˆε
修正)1(ˆX 则称修正后的时间响应式 ⎪⎪⎩
⎪⎪⎨⎧≥-±+-<+-=+----0)]([0)
0()0(0
)
0()1(,))(())1((,))1(()1(ˆ0k k e
a b k a a b e a b x k k a b e a b x k x k k a ak ak εεεεε 为残差修正GM(1,1)模型，简称残差GM(1,1)。

其中残差修正值
)]([0)0()0(0))()(()1(ˆk k a e a b
k a k ----=+εεεεεε
的符号应与残差尾段)0(ε的符号保持一致。

定义6 若)1()0()1()1()0())1()(1()1(ˆ)(ˆ)(ˆ----=--=k a a e a
b
x e k x k x k x
则相应的残差修正时间响应式
⎪⎪⎩
⎪⎪⎨⎧≥-±--<--=+----0)]([0)
0()0(0
)
0()0(,))(())1()(1(,))1()(1()1(ˆ0k k e
a b k a e a b x e k k e a b x e k x k k a ak a ak a εεεεε 称为累减还原式的残差修正模型。

例题 湖北省云梦县油菜发病率数据为
)
15,17,5.15,18,14,21,35,45,40,25,40,20,6())
13(),...,8(),7(),6(),5(),4(),3(),2(),1(()0()0()0()0()0()0()0()0()0()0(==x x x x x x x x x X
建立GM(1,1)模型，得时间响应式为
999.573999.567)1(ˆ06486.0)1(+-=+-k e k x
作累减还原，得 )
4768.176478.188974.192307.216534.221719.247900.255192.273682.293308.314303.336704.35()}(ˆ{ˆ132
)0()0(，，，，，，，，，，，==k x X
检验其精度：列出误差检验表
0尾段
)
4768.2,6478.1,3974.4,2307.3,6534.8())13(),12(),11(),10(),9(()0()0()0()0()0()0(-----==εεεεεε
此为可建模残差尾段，去绝对值，得
11
)4768.2,6478.1,3974.4,2307.3,6534.8()0(=ε
建立GM(1,1)模型，得)0(ε的一次累加序列)1(ε的时间响应式：
7.3224)1(ˆ)9(16855.0)1(+-=+--k e k ε
其导数还原值为
)9(16855.0)9(16855.0)0(0452.4)24)(16855.0()1(ˆ----=--=+k k e e k ε
由k ak a e e a
b x e k x k x k x 06486.0)0()1()1()0(0614.38))1()(1()(ˆ)1(ˆ)1(ˆ--=--=-+=+可得累减还原式残差修正模型为
⎩
⎨⎧≥-<=+----9,0452.40614.389,0614.38)1(ˆ)9(16855.006486.006486.0)0(k e e k e k x k k k 其中，)1(ˆ)0(+k ε
的符号与原始残差序列的符号一致。

按此模型，可对k=10,11,12,13四个模拟值进行休整，修正后的精度如下表：
足建模要求，若对残差精度仍不满意，就只有考虑采用其它模型或对原始数据序列进行适当取舍。