数模自学笔记——灰色模型

参考资料：
《数学建模算法与应用（第2版）》（司守奎主编）
https://wenku.baidu.com/view/7864cb58866fb84ae45c8d65
https://wenku.baidu.com/view/eb12b9d233d4b14e8524689f.html

灰色系统，顾名思义，是白色系统和黑色系统的中间态。

白色系统是指一个系统的内部特征是完全已知的，即系统的信息是完全充分的。
而黑色系统是指一个系统的内部信息对外界来说是一无所知的，只能通过它与外界的联系来加以观测研究。
灰色系统内的一部分信息是已知的，另一部分信息时未知的，系统内各因素间具有不确定的关系。

灰色模型可用于对已知信息少的系统预测未来的发展，适合于“已知过去十年的指标数据，求预测未来十年的指标数据”此类题目。

GM(1,1) 标准模型

模型引入

灰色模型（Gray Model）这里仅讨论常用的 GM(1,1) 标准模型。

已知序列

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

求预测 $$x^{(k)}$$

这个序列的元素即为 $x^{(0)}(i)$ 的预测值。（这其实和拟合有一点 相似，不过拟合可能是已知上百个点，而灰色模型只需要几个点）

首先通过计算累加序列将序列变为光滑、递增的序列：

$$x^{(1)}=\left(x^{(0)}(1), \sum_{i=1}^{2} x^{(0)}(i), \cdots, \sum_{i=1}^{n} x^{(0)}(i)\right)$$

这个过程称为一次累加生成（1-AGO, Accumulating Generation Operator）

接下来建立基本模型 $GM(1,1)$，即白化方程（顾名思义，将灰色模型变成白色模型的方程）是一阶微分方程，即

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

其中，$a$，$b$ 为使用最小二乘法算出的常数，其公式为：

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

其中

$$z^{(1)}(k)=0.5 x^{(1)}(k)+0.5 x^{(1)}(k-1)$$

解上述微分方程可得：

$$
\hat{x}^{(1)}(t+1)=\left[x^{(0)}(1)-b / a\right] \mathrm{e}^{-a t}+b / a$$

然后对累计序列进行还原可得：

$$\hat{x}^{(0)}(t+1)=\hat{x}^{(1)}(t+1)-\hat{x}^{(1)}(t)=\left(1-\mathrm{e}^{a}\right)\left[x^{(0)}(1)-b / a\right] \mathrm{e}^{-a t}$$

模型检验

模型检验有三种方法：残差检验、后验差比值、小误差概率。

残差检验

残差检验需要先计算绝对残差序列：

$$e^{(0)}(i)=\left|x^{(0)}(i)-\hat{x}^{(0)}(i)\right| \quad i=1,2, \ldots, n$$

然后计算相对残差：

$$\varepsilon_i=\frac{e^{(0)}(i)}{X^{(0)}(i)} \times 100 \% \quad i=1,2, \ldots, n$$

求平均得平均相对残差：

$$\bar{\varepsilon}=\frac{1}{n} \sum_{i=1}^{n} \varepsilon_{i}$$

给定 $\alpha$，当 $\overline\varepsilon<\alpha$，且 $\varepsilon_n<\alpha$ 成立，称模型为残差合格模型。

当 $\alpha$ 没有给出时，可取 $\alpha = 1%$。

后验差比值

计算原始序列和残差序列的标准差，分别记作 $S_1$ 和 $S_2$。

$$\begin{aligned}
S_{1}^{2} &=\frac{1}{n} \sum_{k=1}^{n}\left[x^{(0)}(k)-\bar{x}^{(0)}\right]^{2} \\
S_{2}^{2} &=\frac{1}{n} \sum_{k=1}^{n}[e^{(0)}(k)-\bar{e}^{(0)}]^{2}\end{aligned}$$

其中，

$$\bar{x}^{(0)} =\frac{1}{n} \sum_{k=1}^{n} x^{(0)}(k), \bar{e}^{(0)}=\frac{1}{n} \sum_{k=1}^{n} e^{(0)}(k)$$

然后计算后验差比值：

$$C=\frac{S_2}{S_1}$$

后验差比值的精度标准将在后面与小误差概率一并给出。

小误差概率

计算绝对残差序列的均值：

$$\bar{e}^{(0)}=\frac{1}{n} \sum_{k=1}^{n} e^{(0)}(k)$$

然后计算以下概率：

$$p=P\left\{\left|e^{(0)}(k)-\bar{e}^{(0)}\right|<0.6475 S_{1}\right\}$$

下表为灰色模型预测精度等级。

等级精度	后验差比值 $C$	小概率误差 $p$
好	$\leqslant 0.35$	$\geq 0.95$
合格	$0.35<c \leqslant 0.50$	$0.80 \leqslant p<0.95$
勉强	$0.50<c \leqslant 0.65$	$0.70 \leqslant p<0.80$
不合格	$>0.65$	$<0.70$

基于背景值优化的改进 GM(1,1) 模型

该部分由 @Anne_kj 撰写。

设原始非负数据序列为
$$X^{(0)}=\left\{x^{(0)}(1), x^{(0)}(2), \cdots, x^{(0)}(n)\right\}$$
作一次累加 (1 - AGO) 得：

$$X^{(1)}=\left\{x^{(1)}(1), x^{(1)}(2), \cdots, x^{(1)}(n)\right\}$$

其中，$x^{(1)}(k)=\sum_{i=1}^{k} x^{(0)}(i), \quad k=1,2, \cdots, n$

对生成序列 $X^{(1)}(k)$ 建立灰色白化微分方程为：

$$\frac{d x^{(1)}(t)}{d t}+a x^{(1)}(t)=u$$

将上式在区间 $[k, k+1]$ 上积分，得：

$$x^{(1)}(k+1)-x^{(1)}(k)+a \int_{k}^{k+1} x^{(1)}(t) d t=u$$

其中，$k=1,2, \cdots, n-1$。则有：

$$x^{(0)}(k+1)+a \int_{k}^{k+1} x^{(1)}(t) d t=u$$

设 $z^{(1)}({k}+1)=\int_{k}^{k+1} x^{(1)}(t) dt$ 为 $x^{(1)}(t)$ 在区间 $[k, k+1]$ 的背景值，则有

$$x^{(0)}(k+1)+a z^{(1)}(k+1)=u$$

然而函数 $x^{(1)}(\mathrm{t})$ 并不知道，但根据一阶微分方程的解为指数函数，故可令：

$$x^{(1)}(t)=ce^{bt}$$

并假设该曲线过点 $\left(k, x^{(0)}(k)\right)$ 及 $\left(k+1, x^{(0)}(k)\right),$ 则有

$$x^{(1)}(k)=c e^{b k}, \quad x^{(1)}(k+1)=c e^{b(k+1)}$$

由以上两式可解得：

$$\begin{array}{l}
b=\ln x^{(1)}(k+1)-\ln x^{(1)}(k) \\
c=\frac{x^{(1)}(k)}{e^{b k}}=\frac{\left[x^{(1)}(k)\right]^{k+1}}{\left[x^{(1)}(k+1)\right]^{k}}
\end{array}$$

因此背景值为：

$$\begin{aligned}
\mathrm{z}^{(1)}(k+1) &=\int_{k}^{k+1} x^{(1)}(t) d t=\int_{k}^{k+1} c e^{b t} d t \\
&=\frac{x^{(1)}(k+1)-x^{(1)}(k)}{\ln x^{(1)}(k+1)-\ln x^{(1)}(k)}
\end{aligned}$$

待辨参数 $a, u$ 可用最小二乘法求出

$$\hat{\Phi}=[\hat{a} \hat{u}]^{T}=\left(B^{T} B\right)^{-1} B^{T} Y$$

其中

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

$z^{(1)}(k+1)$ 为由前面分析确定的背景值。

求出待辨识参数后，可得白化方程的离散解为：

$$\hat{x}^{(1)}(k+1)=\left[x^{(1)}(1)-\frac{\hat{u}}{\hat{a}}\right] e^{-\hat{a} k}+\frac{\hat{u}}{\hat{a}}$$

进行一次累減还原，可得原始数据序列的模拟值为：

$$\begin{aligned}
\hat{x}^{(0)}(k+1) &=\hat{x}^{(1)}(k+1)-\hat{x}^{(1)}(k) \\
&=\left(1-e^{a}\right)\left[x^{(1)}(1)-\frac{\hat{u}}{\hat{a}}\right] e^{-\hat{a} k}
\end{aligned}$$

GM(1, 1) 新陈代谢模型

若要进一步提高预测精度，可采用 GM(1, 1) 新陈代谢模型。首先采用原始序列建立 一个 GM(1, 1) 模型，按上述方法求出一个预测值，然后将该预测值补入已知数列中，同时 去除一个最旧的数据；在此基础上再建立 GM(1, 1) 模型，求出下一个预测值，以此类推。通过预测灰数的新陈代谢，逐个预测，依次递补，可以得到之后几期的数据，对原始数据数量进行有效扩充。

例题

例题可见本文开头的链接。