灰色预测之GM(1,1)，教你快速上手数学建模！

灰色预测之$GM(1,1)$，教你快速上手数学建模！

前言：在参加数学建模比赛时经常需要大家做预测，而我们常用的预测模型有回归分析预测模型、自回归移动平均模型、灰色系统预测模型、$BP$神经网络预测模型等。今天小编将为大家讲解最传统的灰色预测模型，让大家快速学会预测，上手数学建模。（内附Python完整代码！！）

一、灰色预测简介

小编先简单为大家介绍一下灰色系统理论：

灰色系统理论由中国学者邓聚龙教授于1982年创立，是一种专门用于研究“部分信息已知，部分信息未知”的不确定性系统问题的方法。目前灰色预测模型已经广泛应用于城市环境、交通管理、能源分析等众多领域。

二、$GM(1,1)$理论

接下来小编就给大家讲述最传统的灰色预测理论–$GM(1,1)$,我们一起来看看叭！

（一）、灰色累加和累减生成

定义1(累加生成)：

设原始序列 X ( 0 ) = { X ( 0 ) ( 1 ) , X ( 0 ) ( 2 ) , . . . , X ( 0 ) ( r ) } X^{(0)}=\{{X^{(0)}(1),X^{(0)}(2),...,X^{(0)}(r)}\} X(0)={X(0)(1),X(0)(2),...,X(0)(r)},经过一次累加过后的序列记为 X ( 1 ) = { X ( 1 ) ( 1 ) , X ( 1 ) ( 2 ) , . . . , X ( 1 ) ( r ) } X^{(1)}=\{{X^{(1)}(1),X^{(1)}(2),...,X^{(1)}(r)}\} X(1)={X(1)(1),X(1)(2),...,X(1)(r)}，其中
$$X^{(1)}(k)=\sum _{j=1}^k{X}^{(0)}(j),k=1,2,...,r$$
上述序列中，$X^{(0)}$称为$X^{(1)}$的一阶累加生成序列，上式称为一阶累加生成，简称1-AGO序列。

定义2(累减生成)：

设原始序列的1-AGO序列为 X ( 1 ) = { X ( 1 ) ( 1 ) , X ( 1 ) ( 2 ) , . . . , X ( 1 ) ( r ) } X^{(1)}=\{{X^{(1)}(1),X^{(1)}(2),...,X^{(1)}(r)}\} X(1)={X(1)(1),X(1)(2),...,X(1)(r)}，则
$$X^{(k)}=X^{(1)}(k)-X^{(1)}(k-1)$$
称为$X^{(1)}$的一阶累减生成序列。

上式可用于计算原始序列的预测值$\widehat{X^{(0)}}$。

(二)、GM(1,1)模型

1. 白化方程

   常微分方程
   $$\frac{d{X}^{(1)}(t)}{dt}+a{X}^{(1)}(t)=b$$
   称为GM(1,1)模型的白化方程。

2. 差分方程

   $$X^{(0)}(k)+a{Z}^{(1)}(k)=b$$

   称为GM(1,1)模型的基本形式。其中$a$称为发展系数，$b$称为灰作用量。$Z^{(1)}$称为$X^{(1)}$的紧邻均值序列（背景值序列），满足
   $$Z^{(1)}(k)=\frac{1}{2}[X^{(1)}(k)+X^{(1)}(k-1)],k=2,3,...,r$$

3. GM(1,1)差分方程的参数的最小二乘估计

   GM(1,1)差分方程的参数的最小二乘估计满足
   $$[a,b{]}^T=(B^{T}B{)}^{-1}{B}^T{Y}_R$$
   其中
   $$B={\left[\begin{array}{cccc}-Z^{(1)}(2)& -Z^{(1)}(3)& ...& -Z^{(1)}(r)\\ 1& 1& ...& 1\end{array}\right]}^T$$

   $$Y_R={\left[\begin{array}{cccc}{X}^{(0)}(2)& X^{(0)}(3)& ...& X^{(0)}(r)\end{array}\right]}^T$$

4. 对$GM(1,1)$的白化方程求解

   求解可得
   $$\widehat{X^{(1)}}(t)=(X^{(0)}(1)-\frac{a}{b})e^{-a(t-1)}+\frac{a}{b}$$
   其离散形式
   $$\widehat{X^{(1)}}(k+1)=(X^{(0)}(1)-\frac{a}{b})e^{-a(k)}+\frac{a}{b}$$
   上式可用于计算原始序列预测值的1-AGO序列$\widehat{X^{(1)}}$。

5. 实例分析

   下面我们一起来看看实例吧！

   原始序列 X ( 0 ) = { 4.1 ， 4.05 ， 4.09 ， 4.05 ， 4.02 ， 3.9 ， 3.8 ， 3.62 } X^{(0)}=\{{4.1，4.05，4.09，4.05，4.02， 3.9，3.8，3.62}\} X(0)={4.1，4.05，4.09，4.05，4.02，3.9，3.8，3.62}

   我们用前6个数据进行建模，后面的数据进行预测。并计算真实值和预测值的相对误差，具体代码如下：

   import numpy as np
   import matplotlib.pyplot as plt
   ####编写GM1_1函数
   def GM1_1(y,r)
       m = r+2 #m=r+rf(rf为预测的个数）
       ####生成一次累加序列
       x11 = np.cumsum(y)
       ####生成背景值序列
       z = 0.5 * x11[0:r - 1] + 0.5 * x11[1:r]
       z = z.reshape(-1, 1)
       ####构造B,Y
       E = np.ones(r - 1)
       E = E.reshape(-1, 1)
       B = np.hstack((-z, E))  # 矩阵的拼接
       Y = y[1: r]
       ####参数求解get_parameter
       u = np.matmul(np.linalg.pinv(B), Y)  # 参数的最小二乘估计
       ####求解
       x11_solve = np.zeros(m)
       y_solve = np.zeros(m)
       x11_solve[0] = y_solve[0] = y[0]
       for i in range(1, m):
           x11_solve[i] = (y[0]-u[1]/u[0])*np.exp(-u[0]*i)+u[1]/u[0]#白化方程的解
       ####累减还原
       y_solve = np.diff(x11_solve)
       y_solve = np.r_[y[0], y_solve ]
       return y_solve, x11_solve, u
   
   ####读取数据，运行函数
   y = np.array([4.1，4.05，4.09，4.05，4.02，
   3.9，3.8，3.62])
   r = 6 #用于建模的点数
   y_solve, x11_solve, u = GM1_1(y, r)
   print(f"预测值：{y_solve}")
   
   ####误差估计
   error = np.zeros(len(y))
   for i in range(len(y)):
       error[i] = (y[i] - y_solve[i])/y[i]
   print("Errors: ", error
         
   ####plot
   x = range(1, len(y)+1)
   plt.title('Result Analysis')
   plt.plot(x, y, color='green', label='training accuracy')
   plt.plot(x, y_solve, color='red', label='testing accuracy')
   plt.legend()  # 显示图例
   plt.show()
   

下表为参数求解的结果：

$a$	$b$
0.00912861	4.15189105

下表为真实值和误差值的相对误差：

$Errors(7)$	$Errors(8)$
$-0.02973711$	$-0.07111693$

下图为真实值和预测值的结果比较：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-eHPtJJBx-1611136166137)(灰色预测之GM(1,1)]，教你快速上手数学建模！.assets/Figure_1.png)

以上这些就是关于$GM(1,1)$的介绍啦，你学会了吗？当你数学建模需要预测的时候，你就可以使用它喔，如果有任何问题可以联系我们，之后会为大家介绍灰色预测之$DGM(1,1)$,敬请期待。

链图片转存中…(img-eHPtJJBx-1611136166137)]，教你快速上手数学建模！.assets/Figure_1.png)

以上这些就是关于$GM(1,1)$的介绍啦，你学会了吗？当你数学建模需要预测的时候，你就可以使用它喔，如果有任何问题可以联系我们，之后会为大家介绍灰色预测之$DGM(1,1)$,敬请期待。
更多请关注微信公众号：Python希望社