灰色预测模型——Python

GM(1,1) 预测模型

原理步骤

Step1: 数据检验和处理

设参考数列为 $x^{(0)} = (x^{(0)}(1),x^{(0)}(2),\cdots,x^{(0)}(n))$，计算序列的级比

$$\lambda(k)=\frac{x^{(0)}(k-1)}{x^{(0)}(k)}, k=2,3,\cdots,n.$$

如果所有 $\lambda(k)$ 都落在可容覆盖 $\Theta = (e^{-\frac{2}{n+1}},e^{\frac{2}{n+2}})$，则序列可以作为模型 $GM(1,1)$ 的数据进行灰色预测.
反之，则需要对序列 $x^{(0)}$ 做必要的变换处理，使其落入可容覆盖内，即取适当的常数 $c$，做平移变换

$$y^{(0)}(k) = x^{(0)}(k) + c, k = 1,2,\cdots,n,$$

使序列 $y^{(0)} = (y^{(0)}(1),y^{(0)}(2),\cdots,y^{(0)}(n))$ 的级比

$$\lambda_y(k)=\frac{y^{(0)}(k-1)}{y^{(0)}(k)} \in \Theta, k=2,3,\cdots,n.$$

Step2: 建立模型

（1）对原始数据进行一次累加

$$x^{(1)} = (x^{(0)}(1),x^{(0)}(1)+x^{(0)}(2),\cdots,x^{(0)}(1)+\cdots+x^{(0)}(n))$$

（2）建立灰微分方程和白化微分方程

灰微分方程为

$$x^{(0)}(k) + az^{(1)}(k) = b,\ k=2,3,\cdots,n,$$

其中 $z^{(1)}(k) = 0.5x^{(1)}(k) + 0.5x^{(1)}(k-1)$，相应的白化微分方程为

$$\frac{\mathrm{d}x^{(1)}(t)}{\mathrm{d}t} + ax^{(1)}(t) = b$$

（3）参数估计

记

$$B = \left[ \begin{array}{cc} -z^{(1)}(2) & 1\\ -z^{(1)}(3) & 1\\ \vdots & \vdots\\ -z^{(1)}(n) & 1\\ \end{array} \right],\\ Y = [x^{(0)}(2),x^{(0)}(3),\cdots,x^{(0)}(n)]^T$$

令 $u=[a,b]^T$，则求

$$\arg\min\limits_{u} (Y-Bu)^T(Y-Bu)$$

利用最小二乘法可得

$$\hat{u} = [\hat{a},\hat{b}]^T=(B^TB)^{-1}B^TY$$

（4）模型求解

将参数估计得到的结果代入求解微分方程，得

$$\hat{x}^{(1)}(k+1)=\bigg(x^{(0)}(1)-\frac{\hat{b}}{\hat{a}}\bigg)e^{-\hat{a}k} + \frac{\hat{b}}{\hat{a}}, \quad k=0,1,\cdots,n-1,\cdots$$

且

$$\hat{x}^{(0)}(k+1) = \hat{x}^{(1)}(k+1)-\hat{x}^{(1)}(k)$$

Step3. 模型得检验

（1）相对误差检验.

计算相对误差

$$\delta(k) = \frac{|x^{(0)}(k) - \hat{x}^{(0)}(k)|}{x^{(0)}(k)}$$

这里 $\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) = 1 - \frac{1-0.5a}{1+0.5a}\lambda(k)$$

若 $|\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）建立灰微分方程和白化微分方程

灰微分方程为

$$x^{(0)}(k) + az^{(1)}(k) = b,\ k=2,3,\cdots,n,$$

其中 $z^{(1)}(k) = 0.5x^{(1)}(k) + 0.5x^{(1)}(k-1)$，相应的白化微分方程为

$$\frac{\mathrm{d}x^{(1)}(t)}{\mathrm{d}t} + ax^{(1)}(t) = b$$

（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]
于是

$$\hat{a} = 0.0023, \hat{b} = 72.6573$$

（4）模型求解

将参数估计得到的结果代入求解微分方程，得

$$\hat{x}^{(1)}(k+1)=\bigg(x^{(0)}(1)-\frac{\hat{b}}{\hat{a}}\bigg)e^{-\hat{a}k} + \frac{\hat{b}}{\hat{a}}, \quad k=0,1,\cdots,n-1,\cdots$$

且

$$\hat{x}^{(0)}(k+1) = \hat{x}^{(1)}(k+1)-\hat{x}^{(1)}(k)$$

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)$

$$\hat{x}^{(0)}(k+1) = \hat{x}^{(1)}(k+1)-\hat{x}^{(1)}(k)$$

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

拟合图像：