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

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

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

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

# GM(1,1) 标准模型

# 模型引入

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

已知序列

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

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

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

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

x(1)=(x(0)(1),∑i=12x(0)(i),⋯,∑i=1nx(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)GM(1,1)，即白化方程（顾名思义，将灰色模型变成白色模型的方程）是一阶微分方程，即

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

其中，aa，bb 为使用最小二乘法算出的常数，其公式为：

[a,b]T=(BTB)−1BTYB=[−z(1)(2)1−z(1)(3)1⋮⋮−z(1)(n)1],Y=[x(0)(2)x(0)(3)⋮x(0)(n)]\begin{array}{cc} {[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.5x(1)(k)+0.5x(1)(k−1)z^{(1)}(k)=0.5 x^{(1)}(k)+0.5 x^{(1)}(k-1)

解上述微分方程可得：

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

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

x^(0)(t+1)=x^(1)(t+1)−x^(1)(t)=(1−ea)[x(0)(1)−b/a]e−at\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)=∣x(0)(i)−x^(0)(i)∣i=1,2,…,ne^{(0)}(i)=\left|x^{(0)}(i)-\hat{x}^{(0)}(i)\right| \quad i=1,2, \ldots, n

然后计算相对残差：

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

求平均得平均相对残差：

εˉ=1n∑i=1nεi\bar{\varepsilon}=\frac{1}{n} \sum_{i=1}^{n} \varepsilon_{i}

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

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

# 后验差比值

计算原始序列和残差序列的标准差，分别记作 S1S_1 和 S2S_2。

S12=1n∑k=1n[x(0)(k)−xˉ(0)]2S22=1n∑k=1n[e(0)(k)−eˉ(0)]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}

其中，

xˉ(0)=1n∑k=1nx(0)(k),eˉ(0)=1n∑k=1ne(0)(k)\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=S2S1C=\frac{S_2}{S_1}

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

# 小误差概率

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

eˉ(0)=1n∑k=1ne(0)(k)\bar{e}^{(0)}=\frac{1}{n} \sum_{k=1}^{n} e^{(0)}(k)

然后计算以下概率：

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

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

等级精度 后验差比值 CC 小概率误差 pp
好 ⩽0.35\leqslant 0.35 ≥0.95\geq 0.95
合格 0.35<c⩽0.500.35<c \leqslant 0.50 0.80⩽p<0.950.80 \leqslant p<0.95
勉强 0.50<c⩽0.650.50<c \leqslant 0.65 0.70⩽p<0.800.70 \leqslant p<0.80
不合格 >0.65>0.65 <0.70<0.70

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

该部分由 @Anne_kj (opens new window) 撰写。

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

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

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

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

dx(1)(t)dt+ax(1)(t)=u\frac{d x^{(1)}(t)}{d t}+a x^{(1)}(t)=u

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

x(1)(k+1)−x(1)(k)+a∫kk+1x(1)(t)dt=ux^{(1)}(k+1)-x^{(1)}(k)+a \int_{k}^{k+1} x^{(1)}(t) d t=u

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

x(0)(k+1)+a∫kk+1x(1)(t)dt=ux^{(0)}(k+1)+a \int_{k}^{k+1} x^{(1)}(t) d t=u

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

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

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

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

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

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

由以上两式可解得：

b=ln⁡x(1)(k+1)−ln⁡x(1)(k)c=x(1)(k)ebk=[x(1)(k)]k+1[x(1)(k+1)]k\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}

因此背景值为：

z(1)(k+1)=∫kk+1x(1)(t)dt=∫kk+1cebtdt=x(1)(k+1)−x(1)(k)ln⁡x(1)(k+1)−ln⁡x(1)(k)\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,ua, u 可用最小二乘法求出

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

其中

Y=[x(0)(2)x(0)(3)⋮x(0)(n)],B=[−z(0)(2)1−z(0)(3)1⋮⋮z(0)(n)1]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)z^{(1)}(k+1) 为由前面分析确定的背景值。

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

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

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

x^(0)(k+1)=x^(1)(k+1)−x^(1)(k)=(1−ea)[x(1)(1)−u^a^]e−a^k\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) 模型，求出下一个预测值，以此类推。通过预测灰数的新陈代谢，逐个预测，依次递补，可以得到之后几期的数据，对原始数据数量进行有效扩充。

# 例题

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