模型简介
灰色预测模型使用的不是原始数据序列，而是生成的数据序列。其核心体系是灰色模型，即对原始数据作累加生成（或其他方法生成）得到近似指数规律再进行建模的方法。
灰色预测模型的优点有：
- 仅依靠少量数据就能解决预测问题
- 利用微分方程挖掘系统的本质，精度高
- 运算简便，易于检验
灰色预测模型有以下局限性：
- 只适用于中短期的预测
- 只适合指数增长的预测
GM(1,1)预测模型
GM(1,1)预测模型的简要原理是指：首先利用累加技术使数据具备指数规律，然后建立一阶微分方程并对其求解，将所求结果再累减还原，即为灰色预测值，从而对未来进行预测。下文以一个实例介绍该模型的原理与步骤。
问题
北方某城市1986-1992年道路交通噪声平均声级数据如下：
| 年份 | 1986 | 1987 | 1988 | 1989 | 1990 | 1991 | 1992 |
|---|---|---|---|---|---|---|---|
| 平均声级 | 71.1 | 72.4 | 72.4 | 72.1 | 71.4 | 72.0 | 71.6 |
请预测1993年噪声平均声级数据。
级比检验
由于我们研究的对象是每一年的平均声级数据，因此将建立平均声级时间序列如下：
计算序列的级比公式为
当每一个$k$对应的级比值均落在区间$(e^{-\frac{2}{n+1}},e^{\frac{2}{n+2}})$内时，则序列$x^{(0)}$可以作为模型GM(1,1)的数据进行灰色预测。否则，需要对$x^{(0)}$作平移变换，使其落入区间内。取合适的常数$c$，此时序列的级比为
常数$c$的取值需要不断尝试。级比检验代码如下
1 | # 级比检验 |
该例所有数据均通过级比检验。
建立模型
原始数据序列$x^{(0)}$并没有规律，很难预测出未来数据的值和走向，因此我们要挖掘数据内在规律。通常，我们求原始数据序列的前缀和。
由此得到新的时间序列$x^{(1)}=(71.1,143.5,215.9,288,359.4,431.4,503)$。多数情况下，新的数据序列$x^{(1)}$会呈现指数函数曲线走势，该问题给出的数据较为特殊，呈现直线走势。但是我们依旧可以将一段直线视为指数函数曲线的一部分，即可以利用函数$e^{f(t)}+c$拟合$x^{(1)}$的图像。
拟合函数$e^{f(t)}+c$形式特殊，不难联想到一阶常微分方程的通解。因此，我们可以对$x^{(1)}$建立一阶微分方程，方程的解就是$x^{(1)}$的拟合函数。
由于$x^{(1)}$并非连续函数，而是一个与时间有关的离散的序列，所以$dt=\Delta t=t-(t-1)=1$，$dx^{(1)}=\Delta x^{(1)}=x^{(1)}(t)-x^{(1)}(t-1)=x^{(0)}(t)$。因此，上述一阶微分方程可以转换为以下形式
很多有关文献中，为了提高模型的合理性，$x^{(1)}$被修正为均值生成序列$z^{(1)}=(z^{(1)}(0),z^{(1)}(2),…,z^{(1)}(k))$，其中$z^{(1)}(k)=0.5x^{(1)}(k)+0.5x^{(1)}(k-1)$.最后得到灰微分方程为
式中，$x^{(0)}$与$z^{(1)}$均是已知数据，只要求出$a$和$b$，即可带入微分方程求出$x^{(1)}(t)$的解，也就是数据序列的拟合函数。这里使用最小二乘法，使拟合值与真实值误差的平方和最小。记$\pmb{u}=[a,b]^T$，$\pmb{Y}=[x^{(0)}(2),x^{(0)}(3),…,x^{(0)}(n)]^T$，$\pmb{B}=\begin{bmatrix}-z^{(1)}(2)&1\\ -z^{(1)}(3)&1\\ …&..\\ -z^{(1)}(n)&1\\\end{bmatrix}$，由最小二乘法，求得使$J(\pmb{u})=(\pmb{Y}-\pmb{Bu})^T(\pmb{Y}-\pmb{Bu})$达到最小值的$\pmb{u}$的估计值为
1 | def leastSquare(x, z): |
x为原始数据序列$x^{(0)}$，z为均值生成序列$z^{(1)}$.最后求得$\widehat{a}=0.00234$，$\widehat{b}=72.7$.
最后求解微分方程，得到$\widehat{x}^{(1)}$的函数式
1 | # 求解微分方程，得到预测值序列 |
xEstSum为预测序列$\widehat{x}^{(1)}$，xEst为预测序列$\widehat{x}^{(0)}$.求解结果如下
检验预测值
对于求出来的预测值，需要检验其合理性，有两种检测方法。
残差检验
令残差$\varepsilon (k)$为
且$\widehat{x}^{(0)}(1)=x^{(0)}(1)$，如果$\varepsilon(k)<0.1$，认为预测模型达到要求；如果$\varepsilon(k)<0.2$，认为预测模型达到较高要求。如果残差不在上述范围，说明灰色预测模型不适合求解该问题。
1 | # 残差检验 |
级比偏差值检验
由参考数据$x^{(0)}(k-1)$，$x^{(0)}(k)$计算级比$\lambda(k)$，再计算级比偏差
同理，如果$\rho(k)<0.1$，认为预测模型达到要求；如果$\rho(k)<0.2$，认为预测模型达到较高要求。如果级比偏差不在上述范围，说明灰色预测模型不适合求解该问题。
1 | # 级比偏差值检验 |
结论
模型中各种检验指标值的计算结果如下：
可以看到除了个别检验指标在接受范围外，大部分指标值都落在接受范围内，说明模型的精度较高，可以进行预报。因此1993年平均噪声平均声级为
完整代码如下：1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93'''
-*- coding: utf-8 -*-
Copyright(c) 2022 Zhou Yee
@Author: Zhou Yee
@Date: 2022-09-12 17:08:43
@Email: zhou.yee@foxmail.com
'''
# 灰色预测模型GM(1,1)
import numpy as np
# 级比检验
def preinspect(x):
n = x.shape[0]
interval = (np.exp(-2/(n+1)), np.exp(2/(n+2)))
for k in range(1, n):
if not interval[0]<x[k-1]/x[k]<interval[1]:
return False
return True
# 对数据平移
def translation(x, c):
return x + c
# 求前缀和序列
def prefixSum(x):
xSum = x.copy()
for k in range(x.shape[0]):
xSum[k] = x[0:k+1].sum()
return xSum
# 均值生成序列
def averageArray(xSum):
z = xSum.copy()
z[0] = xSum[0]
for k in range(1, xSum.shape[0]):
z[k] = 0.5*xSum[k] + 0.5*xSum[k-1]
return z
# 最小二乘法求解灰微分方程系数
def leastSquare(x, z):
z = np.expand_dims(z, axis=0)
Y = np.expand_dims(x[1:].copy(), axis=0).T
B = np.insert(z[:,1:].T*(-1), 1, values=np.ones(z.shape[0]), axis=1)
u = np.linalg.inv(B.T@B) @ B.T @ Y
# return a,b
return u[0][0], u[1][0]
# 求解微分方程，得到预测值序列
def estimateArray(x, a, b):
n = x.shape[0]
xEstSum = np.zeros((n))
xEst = np.zeros((n))
for k in range(n):
xEstSum[k] = (x[0] - b/a) * np.exp(-a*k) + b/a
xEst[0] = xEstSum[0]
for k in range(1, n):
xEst[k] = xEstSum[k] - xEstSum[k-1]
return xEst
# 残差检验
def residualTest(x, xEst):
n = x.shape[0]
residual = np.zeros((n))
for k in range(1, n):
residual[k] = (x[k] - xEst[k])/x[k]
return residual
# 级比偏差值检验
def ratioTest(x, xEst, a):
n = x.shape[0]
la = np.zeros((n))
for k in range(1, n):
la[k] = x[k-1]/x[k]
ratio = np.zeros((n))
for k in range(1, n):
ratio[k] = 1 - ((1-0.5*a)/(1+0.5*a)) * la[k]
return ratio
if __name__ == '__main__':
x = np.array([71.1,72.4,72.4,72.1,71.4,72.0,71.6])
if preinspect(x):
xSum = prefixSum(x)
z = averageArray(xSum)
# print(xSum, z)
# print(np.expand_dims(z, axis=0))
a, b = leastSquare(x, z)
print('a:{}\nb:{}'.format(a, b))
xEst = estimateArray(x, a, b)
print('预测序列：{}'.format(xEst))
print('残差检验：{}'.format(residualTest(x, xEst)))
print('级比偏差值检验：{}'.format(ratioTest(x, xEst, a)))
print('1993年预测值：{}'.format(
((x[0] - b/a) * np.exp(-a*7) + b/a) -
((x[0] - b/a) * np.exp(-a*6) + b/a)
))
输出结果如下1
2
3
4
5
6
7
8
9a:0.0023437864785236795
b:72.65726960367881
预测序列：[71.1 72.40574144 72.23623656 72.0671285 71.89841633 71.73009912
71.56217595]
残差检验：[ 0.00000000e+00 -7.93016634e-05 2.26192594e-03 4.55915376e-04
-6.98062087e-03 3.74862332e-03 5.28268865e-04]
级比偏差值检验：[ 0. 0.02025481 0.00234104 -0.0018101 -0.00743993 0.01065487
-0.00323247]
1993年预测值：71.39464589292038