灰色预测 GM(1,1) 模型及Python、R语言实现

GM(1,1)模型

1 GM(1,1)模型概述

灰色预测经常用来解决数据量较少且不能直接发现规律的数据。对于包含不确定信息的序列，灰色预测方法通过对原始数据进行处理，使之转化为灰色序列，并建立微分方程模型模型

GM(1,1)模型是灰色预测理论的基本模型，也是灰色系统理论中运用最广泛的一种动态预测模型，模型由一个单变量的一阶微分方程构成。GM(1,1)模型适合对“少数据，贫信息”的对象进行中短期预测

2 GM(1,1)数据处理方法

在灰色系统中，能获得的信息非常有限，且不易发现内部规律，因此我们需要一些处理，减弱序列的随机性，使得一般规律可以显现。处理方式主要有以下三种：累加生成、累减生成和加权邻值生成

3 GM(1,1)模型方法的可行性检验

（1）级比检验

由原始序列
$$X^{\left(0\right)}=\left(x^{\left(0\right)}\left(1\right),x^0\left(2\right),...,x^0\left(n\right)\right)$$

设 $\lambda \left(k\right)$ 为序列$X^{\left(0\right)}$ 的级比

$$\lambda \left(k\right)=\frac{x^{\left(0\right)}\left(k-1\right)}{x^{\left(0\right)}\left(k\right)}\text{，}k=2,3,\cdots ,n$$若满足：

$$\lambda \left(k\right)\in \left(e^{\frac{-2}{\text{n}+1}},e^{\frac{2}{n+1}}\right)$$
其中 n 为原始序列的数据个数，则$X^{\left(0\right)}$ 可建立GM(1,1) 模型，并进行灰色预测。

否则，需对数据进行处理，使其落在有效范围内，如平移变换：

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

---

（2）光滑度检验

设数据序列 $\left\{x^{\left(0\right)}\left(k\right)\text{，}k=1,2,\cdots ,n\right\}$ 为非负数据序列，则该数据序列为光滑离散序列的充要条件为 $p(k)$ 是 $k$ 的递减函数，即
$$p\left(k\right)=\frac{x^{\left(0\right)}\left(k\right)}{x^{\left(1\right)}\left(k-1\right)}\text{，}k=2,3,\cdots ,n$$
只有原始数据序列满足了这一条件才可以用于建模预测。

4 GM(1,1)预测模型的构建

设原始序列$X^{\left(0\right)}=\left(x^{\left(0\right)}\left(1\right),x^0\left(2\right),...,x^0\left(n\right)\right)$满足上述条件，其中$$x^{\left(0\right)}\left(k\right)\ge 0$$

（1）序列的累加处理
将原始序列进行一次累加生成后处理后，生成 $X^{\left(0\right)}$的1-AGO 序列 (累加生成序列)：
$$X^{\left(1\right)}=\left(x^{\left(1\right)}\left(1\right),x^{\left(1\right)}\left(2\right),...,x^{\left(1\right)}\left(n\right)\right)$$
其中
$$x^{\left(1\right)}\left(k\right)=\sum _{i=1}^k{x}^{\left(0\right)}\left(i\right)k=1,2,......,n$$

（2）计算紧邻均值序列$$Z^{\left(1\right)}=\left(z^{\left(1\right)}\left(2\right),z^{\left(1\right)}\left(3\right),...,z^{\left(1\right)}\left(n\right)\right)$$
其中$$z^{\left(1\right)}\left(k\right)=\frac{1}{2}\left(x^{\left(1\right)}\left(k\right)+x^{\left(1\right)}\left(k-1\right)\right)$$

（3）建立相关模型
建立一阶微分线性方程，即灰色微分方程，得到 GM(1,1)模型的均值形式：$$x^{\left(0\right)}\left(k\right)+a{z}^{\left(1\right)}\left(k\right)=b$$

通过 GM(1,1) 模型相应的白化微分方程：
$$\frac{{\text{dx}}^{\left(1\right)}}{dt}+a{x}^{\left(1\right)}=b$$

其中 a 表示发展系数，反映了${\hat{x}}^{\left(1\right)}$ 和${\hat{x}}^{\left(0\right)}$ 的发展态势； b 表示灰色作用量，或者内生控制灰数，是从行为序列中挖掘出来的数据，反映的是数据变化的关系，其确切内涵是灰的。

（4）计算 a、b

设 $\hat{a}$ 为待估参数向量，令 $\hat{a}=\left[\begin{array}{c}\text{a}\\ b\end{array}\right]$，利用最小二乘法求解，可得：
$$\hat{a}={\left(B^{T}B\right)}^{-1}{B}^{T}Y$$
其中B、Y 分别为
$$\text{B}=\left[\begin{array}{cc}-z^{\left(1\right)}\left(2\right)& 1\\ -z^{\left(1\right)}\left(3\right)& 1\\ ⋮& ⋮\\ -z^{\left(1\right)}\left(n\right)& 1\end{array}\right]\text{，Y}=\left[\begin{array}{c}{x}^{\left(0\right)}\left(2\right)\\ x^{\left(0\right)}\left(3\right)\\ ⋮\\ x^{\left(0\right)}\left(\text{n}\right)\end{array}\right]$$

GM(1,1)模型中的参数 $-a$ 为发展系数，$b$ 为灰色作用量。

$-a$ 反应了 $X^{(1)}$ 及原始序列 $X^{(0)}$ 的发展趋势。

作用量在物理学中是一个比较抽象的概念，它表示一个动力系统内在的演化趋势。

对原始数据序列进行计算求出的GM(1,1)模型中的 $b$ ，反应了原始数据的变化规律。

（5）建立 GM(1,1) 的时间响应式，即预测模型
$$\{\begin{array}{l}{\hat{x}}^{\left(1\right)}\left(k+1\right)=\left(x^{\left(0\right)}\left(1\right)-\frac{b}{a}\right)e^{-ak}+\frac{b}{a},\text{k}=1,2\cdots ,\text{n}\\ \\ {\hat{x}}^{\left(0\right)}\left(k+1\right)={\hat{x}}^{\left(1\right)}\left(k+1\right)-{\hat{x}}^{\left(1\right)}\left(k\right)\end{array}$$

5 GM(1,1) 模型的检验

检验灰色预测模型效果主要包括三种方法：残差检验、关联度检验和后验差检验方法，以此来检验模型的准确性。

（1） 残差检验

按预测模型计算 ${\hat{x}}^{\left(1\right)}\left(k\right)$，并将 ${\hat{x}}^{\left(1\right)}\left(i\right)$累减生成 ${\hat{x}}^{\left(0\right)}\left(k\right)$，然后计算原始序列 $x^{\left(0\right)}\left(k\right)$ 与 ${\hat{x}}^{\left(0\right)}\left(k\right)$的绝对误差序列即相对误差序列。

$$\{\begin{array}{l}\text{残差：}\epsilon \left(k\right)=x^{\left(0\right)}\left(k\right)-{\hat{x}}^{\left(0\right)}\left(k\right)\\ \\ \text{相对误差：}{\Delta }_k=\frac{|\epsilon \left(k\right)|}{x^{\left(0\right)}\left(k\right)}\\ \\ \text{平均相对误差：}\Delta =\frac{1}{n-1}\sum _{k=2}^n{\Delta }_k\end{array}$$

若对相应序列残差有 ${\Delta }_k<10\%$ ，则代表用 GM(1,1) 模型的模拟精度较高，若 ${\Delta }_k<20\%$ 则代表用 GM(1,1) 模型的模拟精度一般。

GM(1,1) 模型建模精度：
$$p=\left(1-\Delta \right)\times 100\%$$

（2）关联度检验法

绝对误差：
$${\Delta }^{\left(0\right)}\left(k\right)=\left|x^{\left(0\right)}\left(k\right)-{\hat{x}}^{\left(0\right)}\left(k\right)\right|$$

得出最大值、最小值

$$\{\begin{array}{l}{\Delta }_{\max }^{\left(0\right)}=\max {\Delta }^{\left(0\right)}\left(k\right)\\ \\ {\Delta }_{\min }^{\left(0\right)}=\min {\Delta }^{\left(0\right)}\left(k\right)\end{array}$$
计算关联系数
$$\eta \left(k\right)=\frac{{\Delta }_{\max }^{\left(0\right)}+\rho {\Delta }_{\min }^{\left(0\right)}}{{\Delta }^{\left(0\right)}\left(k\right)+\rho {\Delta }_{\max }^{\left(0\right)}}\text{，}\rho =0.5\text{，}k=1,2,\cdots ,n$$
得出关联度
$$r=\frac{1}{n}\sum _{k=1}^n\eta \left(k\right)$$
根据经验，当 p =0.5 时，关联度 r >0.6 时可认为模型的预测效果较好。

（3）后验差检验
原始序列的均值
$${\bar{x}}^{\left(0\right)}=\frac{1}{n}\sum _{k=1}^n{x}^{\left(0\right)}\left(k\right)$$

计算的原始序列的标准差：
$$S_1=\sqrt{\frac{\sum {\left[x^{\left(0\right)}\left(k\right)-{\bar{x}}^{\left(0\right)}\right]}^2}{n-1}}$$

残差序列的均值：
$${\Delta }^{\left(0\right)}=\frac{1}{n}\sum _{k=1}^n\left[x^{\left(0\right)}\left(k\right)-{\hat{x}}^{\left(0\right)}\left(k\right)\right]$$
计算绝对误差序列的标准差：
$$S_2=\sqrt{\frac{\sum {\left[{\Delta }^{\left(0\right)}\left(k\right)-{\Delta }^{\left(0\right)}\right]}^2}{n-1}}$$
计算后验差比值为：：
$$C=\frac{S_2}{S_1}$$
计算小误差概率：
$$P=p\left\{\left|{\Delta }^{\left(0\right)}\left(k\right)-{\Delta }^{\left(0\right)}\right|<0.6745S_1\right\}$$

GM(1,1) 通常用后验差检验方法来评价预测结果的好坏，主要根据其中的后验差比值（C）和小误差概率（P）这两个数值来检验模型。根据 C 和 P 的大小可以将灰色模型的预测精度分为以下几个等级，如下表

后验差比值 C	小误差概率 P	预测精度等级
<0.35	>0.95	好（一级）
<0.50	>0.80	合格（二级）
<0.65	>0.70	勉强合格（三级）
>=0.65	<=0.65	不合格 （四级）

6 GM(1,1) 模型的适用范围

建立于灰色理论基础上的GM(1,1)模型具有所需样本量小，无需考虑数据分布规律，模型预测精度较高，且易于计算检验，决定了GM(1,1)模型被广泛的应用于经济生活各个领域。但这并非说明GM(1,1)模型应用具有随意性与其他任何数学模型相同GM(1,1)模型也有其适用范围，若超出其适用范围使用GM(1,1)模型往往很难得到满意的的结果。

GM(1,1)模型的参数与GM(1,1)模型适用性之间的关系得到如下结论

（1）当 -a ≤ 0.3时，GM(1,1)模型可用于中长期预测；

（2）当 0.3 < -a ≤ 0.5时，GM(1,1)模型可用于短期预测，中长期预测需谨慎；

（3）当 0.5 < -a ≤ 0.8时，用GM(1,1)模型作短期预测慎用；

（4）当 0.8 < -a ≤ 1.0时，需采用残差修正GM(1,1)模型；

（5）当 1 < -a 时，不宜使用GM(1,1)模型。

也即

模型预测时长判断标准

7 GM(1,1) 残差模型

当原始数据序列$x^{\left(0\right)}\left(k\right)$建立的GM(1,1)模型检验不合格时，可以用GM(1,1)残差模型来修正。如果原始序列建立的GM(1,1)模型不够精确，也可以用GM(1,1)残差模型来提高精度。

若用原始序列 $x^{\left(0\right)}\left(k\right)$ 建立的GM(1,1)模型

$${\hat{x}}^{\left(1\right)}\left(k+1\right)=\left(x^{\left(0\right)}\left(1\right)-\frac{b}{a}\right)e^{-ak}+\frac{b}{a},\text{k}=1,2\cdots ,\text{n}$$

可获得生成序列$x^{\left(1\right)}\left(k\right)$的预测值，定义残差序列$e^{\left(0\right)}\left(k\right)=x^{\left(1\right)}\left(k\right)-{\hat{x}}^{\left(1\right)}\left(k\right)$。则对应的残差序列为：$e^{\left(0\right)}\left(k\right)=\left\{e^{\left(0\right)}\left(1\right),e^{\left(0\right)}\left(2\right),\cdots ,e^{\left(0\right)}\left(n\right)\right\}$

计算其生成序列$e^{\left(1\right)}\left(k\right)$，并据此建立相应的GM(1,1)模型：

$${\hat{e}}^{\left(1\right)}\left(k+1\right)=\left(e^{\left(0\right)}\left(1\right)-\frac{b_e}{a_e}\right)e^{-ak}+\frac{b_e}{a_e},\text{k}=1,2\cdots ,\text{n}$$
得修正模型
$${\hat{x}}^{\left(1\right)}\left(k+1\right)=\left(x^{\left(0\right)}\left(1\right)-\frac{b}{a}\right)e^{-ak}+\frac{b}{a}+\delta \left(k-i\right)\left(-a_e\right)\left[e^{\left(0\right)}\left(1\right)-\frac{b_e}{a_e}\right]e^{-ak}$$
其中
$$\delta \left(k-i\right)=\{\begin{array}{l}1\text{，}k\ge i\\ 0\text{，}k<i\end{array}$$

为修正系数。

应用此模型时要考虑：

（1）一般不是使用全部残差数据来建立模型，而只是利用了部分残差。

（2）修正模型所代表的是差分微分方程，其修正作用与$\delta \left(k-i\right)$ 中的 i 的取值有关。

R语言实现

灰色系统模型GM(1,1)的R语言实现

灰色预测GM11，R语言源程序

python

'''灰色预测函数'''
def GM11(x0): #自定义灰色预测函数
    import numpy as np
    x1 = x0.cumsum() # 生成累加序列
    z1 = (x1[:len(x1)-1] + x1[1:])/2.0 # 生成紧邻均值（MEAN）序列，比直接使用累加序列好，共 n-1 个值
    z1 = z1.reshape((len(z1),1))
    B = np.append(-z1, np.ones_like(z1), axis = 1)    # 生成 B 矩阵
    Y = x0[1:].reshape((len(x0)-1, 1))    # Y 矩阵
    [[a],[u]] = np.dot(np.dot(np.linalg.inv(np.dot(B.T, B)), B.T), Y)    #计算参数
    f = lambda k: (x0[0]-u/a)*np.exp(-a*(k-1))-(x0[0]-u/a)*np.exp(-a*(k-2))    #还原值
    delta = np.abs(x0 - np.array([f(i) for i in range(1,len(x0)+1)]))    # 计算残差
    C = delta.std()/x0.std()
    P = 1.0*(np.abs(delta - delta.mean()) < 0.6745*x0.std()).sum()/len(x0)
    return f, a, u, x0[0], C, P #返回灰色预测函数、a、b、首项、方差比、小残差概率