灰色预测模型——Python
GM(1,1) 预测模型
原理步骤
Step1: 数据检验和处理
设参考数列为 \(x^{(0)} = (x^{(0)}(1),x^{(0)}(2),\cdots,x^{(0)}(n))\)，计算序列的级比
如果所有 \(\lambda(k)\) 都落在可容覆盖 \(\Theta = (e^{-\frac{2}{n+1}},e^{\frac{2}{n+2}})\)，则序列可以作为模型 \(GM(1,1)\) 的数据进行灰色预测.
反之，则需要对序列 \(x^{(0)}\) 做必要的变换处理，使其落入可容覆盖内，即取适当的常数 \(c\)，做平移变换
使序列 \(y^{(0)} = (y^{(0)}(1),y^{(0)}(2),\cdots,y^{(0)}(n))\) 的级比
Step2: 建立模型
（1）对原始数据进行一次累加
（2）建立灰微分方程和白化微分方程
灰微分方程为
其中 \(z^{(1)}(k) = 0.5x^{(1)}(k) + 0.5x^{(1)}(k-1)\)，相应的白化微分方程为
（3）参数估计
记
令 \(u=[a,b]^T\)，则求
利用最小二乘法可得
（4）模型求解
将参数估计得到的结果代入求解微分方程，得
且
Step3. 模型得检验
（1）相对误差检验.
计算相对误差
这里 \(\hat{x}^{(0)}(1) = x^{(0)}(1)\)，若 \(\delta(k) < 0.2\)，则可认为达到一般要求；若 \(\delta(k) < 0.1\)，则认为达到较高要求.
（2）级比偏差值检验
由参考数列 \(x^{(0)}(k-1),x^{(0)}(k)\) 计算出级比 \(\lambda(k)\)，再用发展系数 \(a\) 求出相应的级比偏差
若 \(|\rho(k)| < 0.2\)，则可认为达到一般要求；若 \(|\rho(k)| < 0.1\)，则认为达到较高要求.
Step4. 预测
根据所建立得 \(GM(1,1)\) 模型预测即可.
实例
某城市 1986-1992 年道路噪声平均声级如下表
| 序号 | 年份 | \(L_{eq}\) |
|---|---|---|
| 1 | 1986 | 71.1 |
| 2 | 1986 | 72.4 |
| 3 | 1986 | 72.4 |
| 4 | 1986 | 72.1 |
| 5 | 1986 | 71.4 |
| 6 | 1986 | 72.0 |
| 7 | 1986 | 71.6 |
试用 \(GM(1,1)\) 拟合.
Step1. 级比检验
定义算法：
from cmath import exp
def Grade_ratio_test(arr): # 传入 ndarray 型序列
n = len(arr)
t1 = exp(-2/(n+1))
t2 = exp(2/(n+2))
arr0 = arr[0:n]
arr1 = arr[1::]
lam = arr0 / arr1
for i in lam:
if i <= t1 or i >= t2:
return ("级比检验不通过！")
return("级比检验通过！")
测试：
import numpy as np
x0 = np.array((71.1, 72.4, 72.4, 72.1, 71.4, 72.0, 71.6))
Grade_ratio_test(x0)
'级比检验通过！'
因此级比检验通过，可以用该数据进行 \(GM(1,1)\) 模型的检验.
Step2. 建模
（1）对原始数据进行一次累加
l1 = []
for i in range(len(x0)):
l1.append(sum(x0[0:i+1]))
x1 = np.array(l1)
print(x1)
[ 71.1 143.5 215.9 288. 359.4 431.4 503. ]
于是 \(x^{(1)} = (71.1, 143.5, 215.9, 288, 359.4, 431.4, 503)\)
（2）建立灰微分方程和白化微分方程
灰微分方程为
其中 \(z^{(1)}(k) = 0.5x^{(1)}(k) + 0.5x^{(1)}(k-1)\)，相应的白化微分方程为
（3）参数估计
x11 = x1[0:len(x1)-1]
x12 = x1[1::]
z1 = (x11+x12)/2
B = np.empty((6,2))
B[:,1] = np.ones(6)
B[:,0] = -z1
Y = x0[1::]
u = np.dot(np.dot(np.linalg.inv(np.dot(B.T,B)),B.T),Y)
print(u)
[2.34378648e-03 7.26572696e+01]
于是
（4）模型求解
将参数估计得到的结果代入求解微分方程，得
且
def GM11(a,b,k):
result = round(np.real((x0[0]-b/a)*exp(-a*k) + b/a),4)
return result
x1_hat = np.empty(6)
a = 0.0023
b = 72.6573
for i in range(len(x1)-1):
x1_hat[i] = (GM11(a, b, i))
print(x1_hat)
[ 71.1 143.5105 215.7546 287.8327 359.7453 431.4926]
再将 \(\hat{x}(k+1)\) 还原为 \(\hat{x}^{(0)}(k+1)\)
x1_hat1 = x1_hat[0:len(x1_hat)-1]
x1_hat2 = x1_hat[1::]
x0_hat0 = x1_hat2 - x1_hat1
x0_hat = np.empty(7)
x0_hat[0] = x1_hat[0]
x0_hat[1::] = x0_hat0
print(x0_hat)
[71.1 72.4057 72.2363 72.0671 71.8984 71.7301 71.5622]
Step3. 模型的检验
def residual(x0, x0_hat):
'''残差'''
return x0 - x0_hat
def relative_error(x0,x0_hat):
'''相对误差'''
return np.abs(x0-x0_hat)/x0
def Grade_ratio_deviation(a, lam):
'''级比偏差'''
rho = 1 - (1-0.5*a)/(1+0.5*a)*lam
return rho
residual(x0,x0_hat)
array([ 0. , -0.0057, 0.1637, 0.0329, -0.4984, 0.2699, 0.0378])
relative_error(x0, x0_hat)
array([0.00000000e+00, 7.87292818e-05, 2.26104972e-03, 4.56310680e-04, 6.98039216e-03, 3.74861111e-03, 5.27932961e-04])
Grade_ratio_deviation(a,lam)
array([ 0.02025481, 0.00234104, -0.0018101 , -0.00743993, 0.01065487, -0.00323247])
于是可得如下表
| 序号 | 年份 | 原始值 | 预测值 | 残差 | 相对误差 | 级比偏差 |
|---|---|---|---|---|---|---|
| 1 | 1986 | 71.1 | 71.1 | 0 | 0 | |
| 2 | 1987 | 72.4 | 72.4057 | -0.0057 | 0.01% | 0.0203 |
| 3 | 1988 | 72.4 | 72.2362 | 0.1638 | 0.23% | 0.0023 |
| 4 | 1989 | 72.1 | 72.0671 | 0.0329 | 0.05% | -0.0018 |
| 5 | 1990 | 71.4 | 71.8984 | -0.4984 | 0.7% | -0.0074 |
| 6 | 1991 | 72.0 | 71.7301 | 0.2699 | 0.37% | 0.0107 |
| 7 | 1992 | 71.6 | 71.5622 | 0.0378 | 0.05% | -0.0032 |
拟合图像：
