基于改进的灰色NGM(1,1)预测模型研究及应用——以新发传染病为例

Research and Application of Improved Gray NGM(1,1) Based Prediction Model—Taking Emerging Infectious Diseases as an Example

摘要: 2020年突发传染病严重威胁人们的生命健康，精准的预测对于传染病的管理有重要意义。基于传染病发展的不确定性，建立传统灰色NGM(1,1)模型，而NGM(1,1)模型初始值的选取、背景值的构造以及参数估计方法是导致该模型精度不稳定的重要原因。首先，本文以三参数的紧邻均值背景值来优化传统NGM(1,1)模型背景值为基础，以平均相对误差最小为目标，利用粒子群优化算法求最优初始值的取值，然后将三参数的紧邻均值背景值改为加权的三参数背景值，利用粒子群优化算法寻求最优的加权权数，最后将参数估计方法最小二乘法改为加权最小二乘法求模型的灰参数，构成NGM(1,1)模型初始值的选取、背景值的构造和参数估计方法的组合优化。通过算例来验证优化后的灰色NGM(1,1)模型性能，并将改进的模型应用于传染病的预测。

Abstract: The outbreak of infectious diseases in 2020 is a serious threat to people’s life and health, and accu-rate prediction is important for the management of infectious diseases. Based on the uncertainty of infectious disease development, the traditional gray NGM(1,1) model is established, and the selec-tion of initial values of NGM(1,1) model, the construction of background values and the parameter estimation methods are important reasons for the unstable accuracy of this model. Firstly, this pa-per optimizes the background value of the traditional NGM(1,1) model based on the three-parameter tight-neighborhood mean background value, and uses the particle swarm optimi-zation algorithm to find the optimal initial value with the objective of minimizing the average rela-tive error, then changes the three-parameter tight-neighborhood mean background value to a weighted three-parameter background value, and uses the particle swarm optimization algorithm to seek the optimal weighting weights, as well as finally changes the parameter estimation method least squares Finally, the least squares method of parameter estimation is changed to the weighted least squares method to find the gray parameters of the model, which constitutes the optimization of the combination of the initial value selection, background value construction and parameter es-timation method of the NGM(1,1) model. The performance of the optimized gray NGM(1,1) model is verified by arithmetic examples, and the improved model is applied to the prediction of infectious diseases.

1. 引言

灰色预测模型是解决“少数据，贫信息”的小数据的重要方法，灰色预测模型以不需要大量的数据样本，短期预测效果好，运算过程简单的特点，且在已知数据少的情况下仍然可以较好的拟合与预测数据的发展趋势，所以灰色预测广泛应用于各个领域 [1] [2] [3] 。朱晓宵等 [4] 以武汉2020年新冠肺炎情况为例，首先构建滚动灰色GM(1,1)模型，用来预测疫情的变化趋势，并与传染病动力学模型预测精度比较，得出滚动GM(1,1)模型 > SEIH1H2RD传染病动力学模型，这说明对于具有不确定性的传染病的传播趋势，灰色预测模型有很好的预测效果。GM(1,1)模型是目前应用最为广泛的灰色预测模型，但是据学者研究，经典GM(1,1)有时也会出现较大的误差，导致经典的GM(1,1)模型的模拟精度不理想，而模拟精度主要与模型参数、模型结构和原始数据序列等有关，针对这一问题，目前学者们主要从模型参数优化和模型结构优化这两个方面研究了GM(1,1)模型的优化问题。如何提高灰色预测模型来预测事物的发展情况也是本文的研究目标。

在模型结构优化方面，考虑到经典的GM(1,1)预测模型参数估计来源于GM(1,1)模型的差分形式，时间响应函数来源于微分方程，造成二者的“非同源性”，为了实现模型参数估计和时间响应函数来源相同，谢乃明等提出了离散的GM(1,1)模型，即DGM(1,1)模型，并且得到了很好的应用 [5] [6] 。不管是GM(1,1)模型还是DGM(1,1)模型都只在近似齐次指数数据序列上可以得到较好的模拟效果，对于近似非齐次指数数据的模拟效果不佳，所以有了GM(1,1)模型和DGM(1,1)模型的很多拓展。张和平等 [7] 将可以描述趋势为S型过程的Verhulst模型应用于网络舆情的预测，得出模型的适用性；邱慧 [8] 等将Verhulst模型应用于人口老龄化的预测；石季雨等 [9] 等将对于近似非齐次指数数据模拟效果好的NGM(1,1)模型应用于中国能源消费情况的预测。

在模型参数优化方面，主要包括初始值优化、背景值优化和参数估计方法优化三个方面，这也是目前对于灰色预测模型优化最多的三个方面。经典的灰色预测模型是选择一阶累加生成序列的第一个数据作为模型的初始值，而这个选择并不能突出新信息的重要性。所以考虑到新信息的重要性，Zhao等 [10] 引入滚动机制；也有学者引入新陈代谢模型 [11] [12] ，用每天更新的新数据来替换最初始数据，给予新数据更高的优先级；丁松等 [13] 考虑新信息优先原理，引入权重对生成序列加权求和来优化GM(1,1)幂模型初始值，然后利用平均相对误差最小来确定权重值；郑坚等 [14] 同样考虑了生成序列的综合影响，按生成数据的大小对数据赋权值，并且引入了时间权重函数来反应各个误差平方和的重要性。经典的GM(1,1)模型的背景值是采用紧邻均值构造法，构造出的梯形面积就是经典的灰色模型的背景值，该值与实际的曲边梯形的面积之差就是经典GM(1,1)模型背景值构造的误差值。为减小经典背景值选取的误差，卢捷等 [15] 将背景值设为变量，然后用平均相对误差最小来确定背景值参数，证明了优化后的模型精度提高，同时也对初始值重新取值对模型进行优化；杨孝良等 [16] 考虑到原始序列中极端值的影响，提出了三参数紧邻均值来构造背景值，弱化了原始序列中极端值的影响；刘震等 [17] 对NGM(1,1)的微分方程进行积分，通过分析误差产生原因来重新推导背景值表达式。对于灰色预测模型的参数估计方法，经典GM(1,1)模型的参数估计方式是最小二乘法，最小二乘法要求数据总体为正态分布，而经典GM(1,1)模型的数据总体特征不明确，且最小二乘法存在稳健性差的问题。许多研究者对此进行了改进，比如何霞 [18] 利用加权最小二乘法来提高模型的精度和稳健性；周德强 [19] 利用LS-SVM来估计灰色Verhulst模型的参数，符合灰色Verhulst模型小样本的建模，提高了灰色Verhulst预测模型的精度。

以上研究都是从灰色预测误差主要来源的一个方面或者两个方面对灰色预测模型进行优化，在一定程度上提高了模型的精度和减少了预测模型误差，但是没有总体解决误差来源问题，所以基于以上研究，本研究首先构建NGM(1,1)模型，从背景值、初始值和参数估计方法三个方面对NGM(1,1)灰色预测模型进行组合优化。用优化后的灰色NGM(1,1)预测模型进行算例分析，计算出模型的相对平均误差、方差比(C)和小误差概率(p)三个指标值，通过与其他模型对比三个指标值的大小，比较模型的性能。最后将本文构建的组合优化的模型应用于武汉市新型冠状病毒肺炎新增感染人数模拟，分别采用不同模型建模并与本文模型进行对比，结果表明所提出的优化模型能够获得较高精度。

2. NGM(1,1)模型优化

NGM(1,1)模型是基于经典GM(1,1)模型构建的一种能够同时模拟近似齐次指数数据和近似非齐次指数数据的灰色预测模型，而GM(1,1)模型是针对近似齐次指数数据有较为理想的预测效果，由于现实生活中大多数数据都是近似非齐次指数的，近似齐次指数数据只是一种理想状态，所以基于现实数据，NGM(1,1)模型具有更好的预测效果。本节基于文献 [16] 依次从初始值取值、背景的构造和参数估计方法三个方面来优化NGM(1,1)模型。

2.1. 三参数背景值的NGM(1,1)模型

NGM(1,1)模型 [20] 原始数列为： $X^{(0)}=\left(x^{(0)}\left(1\right),x^{(0)}\left(2\right),\cdots ,x^{(0)}\left(n\right)\right)$ ，累加生成

$X^{(1)}=\left(x^{(1)}\left(1\right),x^{(1)}\left(2\right),\cdots x^{(1)}\left(n\right)\right)$ ，其中 $x^{(1)}\left(k\right)={\displaystyle \underset{i=1}{\overset{k}{\sum }}{x}^{(0)}\left(k\right)}$ ，其中 $k=1,2,\cdots ,n$ ，用一阶微分方程拟合 $X^{(1)}$ ，

得到NGM(1,1)的白化微分方程形式：

$\frac{\text{d}{x}^{(1)}}{\text{d}t}+a{x}^{(1)}=bt+c$ (1)

生成紧邻均值生成序列 $Z^{(1)}=\left(z^{(1)}\left(2\right),z^{(1)}\left(3\right),\cdots ,z^{(1)}\left(n\right)\right)$ ，其中

$z^{(1)}\left(k\right)=\frac{1}{2}\left(x^{(1)}\left(k\right)+x^{(1)}\left(k-1\right)\right),k=2,3,\cdots ,n$ ，作为模型的背景值，已知数据是离散的而不是连续的，将 $\frac{\text{d}{x}^{(1)}}{\text{d}t}\to \frac{\Delta x^{(1)}}{\Delta t}=\frac{x^{(1)}(k)-x^{(1)}\left(k-1\right)}{t-\left(t-1\right)}={\displaystyle \underset{i=1}{\overset{k}{\sum }}{x}^{(0)}\left(i\right)}-{\displaystyle \underset{i=1}{\overset{k-1}{\sum }}{x}^{(0)}\left(i\right)}=x^{(0)}\left(k\right)$ ， $k=2,3,\cdots ,n$ ，和背景值带入式(1)得到

NGM(1,1)的基本形式：

$x^{(0)}\left(k\right)+a{z}^{(1)}\left(k\right)=bk+c$ (2)

式(2) a、b和c为待估灰参数。根据杨孝良等 [16] 为了提高NGM(1,1)模型背景值的平滑效果，减少模型中极端值对模型精度的影响，提出三参数紧邻均值作为NGM(1,1)模型的背景值，即：

$z^{(1)}\left(k\right)=\frac{1}{3}\left(x^{(1)}\left(k\right)+x^{(1)}\left(k-1\right)+x^{(1)}\left(k-2\right)\right)$ (3)

其中 $k=3,4,\cdots ,n$ ，代替传统NGM(1,1)背景值 $z^{(1)}\left(k\right)=\frac{1}{2}\left(x^{(1)}\left(k\right)+x^{(1)}\left(k-1\right)\right)$ ，其中 $k=2,3,\cdots ,n$ ，将三

参数背景值带入式(2)，利用最小二乘法求解参数a、b和c。将 $x^{(0)}\left(k\right)=x^{(1)}\left(k\right)-x^{(1)}\left(k-1\right)$ 和式(3)带入式(2)得：

${\hat{x}}^{(1)}\left(k\right)=\frac{1-\frac{a}{3}}{1+\frac{a}{3}}{x}^{(1)}\left(k-1\right)-\frac{\frac{a}{3}}{1+\frac{a}{3}}{x}^{(1)}\left(k-2\right)+\frac{b}{1+\frac{a}{3}}k+\frac{c}{1+\frac{a}{3}}$ (4)

令 $a_1=\frac{1-\frac{a}{3}}{1+\frac{a}{3}}$ ， $a_2=-\frac{\frac{a}{3}}{1+\frac{a}{3}}$ ， $a_3=\frac{b}{1+\frac{a}{3}}$ ， $a_4=\frac{c}{1+\frac{a}{3}}$ ，得到离散的NGM(1,1)模型：

${\hat{x}}^{(1)}\left(k\right)=a_1{x}^{(1)}\left(k-1\right)+a_2{x}^{(1)}\left(k-2\right)+a_{3}k+a_4$ (5)

由于式(5)是三阶差分方程求解很困难，所以直接利用递推算法进行求解 ${\hat{x}}^{(1)}\left(k\right)$ ，然后利用 $x^{(0)}\left(k\right)$ 与

$x^{(1)}\left(k\right)$ 和 $x^{(1)}\left(k-1\right)$ 关系求 $x^{(0)}\left(k\right)$ 的预测值 ${\hat{x}}^{(0)}\left(k\right)$ 。

2.2. 优化三参数背景值NGM(1,1)模型的初始值

杨孝良等 [16] 中是将具有三参数背景值的NGM(1,1)模型的初始值 ${\hat{x}}^{(1)}\left(1\right)$ 和 ${\hat{x}}^{(1)}\left(2\right)$ 看为已知数据，即令 ${\hat{x}}^{(1)}\left(1\right)=x^{(1)}\left(1\right)$ ， ${\hat{x}}^{(1)}\left(2\right)=x^{(1)}\left(2\right)$ ，直接给预测值赋原值，虽然序列前两个值的没有预测误差，但是初始值不是参与模型运算得出可能会导致后面的预测值出现较大的误差，所以这里将初始值 ${\hat{x}}^{(1)}\left(1\right)$ 和 ${\hat{x}}^{(1)}\left(2\right)$ 设为未知，令：

${\hat{x}}^{(1)}\left(1\right)={\xi }_1$ (6)

${\hat{x}}^{(1)}\left(2\right)={\xi }_2$ (7)

利用式(5)和 $x^{(0)}\left(k\right)$ 与 $x^{(1)}\left(k\right)$ 和 $x^{(1)}\left(k-1\right)$ 关系求预测值 ${\hat{x}}^{(0)}\left(k\right)$ ，然后以平均相对误差式(8)最小为目标，利用粒子群优化算法确定 ${\hat{x}}^{(1)}\left(1\right)$ 和 ${\hat{x}}^{(1)}\left(2\right)$ 的最优值 ${\xi }_1$ 和 ${\xi }_2$ 。

$fitness=\frac{1}{n-1}{\displaystyle \underset{k=1}{\overset{n}{\sum }}\left|\frac{x^{(0)}\left(k\right)-{\hat{x}}^{(0)}\left(k\right)}{x^{(0)}\left(k\right)}\right|}$ (8)

2.3. 优化三参数背景值NGM(1,1)模型的背景值

杨孝良等 [16] 提出的三参数背景值与传统的NGM(1,1)模型的背景值构造一样是取的紧邻生成序列的均值，为了提高模型预测精度，将紧邻均值生成序列的均值改为未知的加权数，即将式(3)改为：

$z^{(1)}\left(k\right)=\alpha x^{(1)}\left(k\right)+\beta x^{(1)}\left(k-1\right)+\gamma x^{(1)}\left(k-2\right)$ (9)

并满足 $\alpha +\beta +\gamma =1$ 。将式(9)和 $x^{(0)}\left(k\right)$ 与 $x^{(1)}\left(k\right)$ 和 $x^{(1)}\left(k-1\right)$ 的关系式代入式(2)得：

$x^{(1)}\left(k\right)=\frac{1-\alpha \beta }{1+a\alpha }{x}^{(1)}\left(k-1\right)-\frac{\gamma }{1+a\alpha }{x}^{(1)}\left(k-2\right)+\frac{b}{1+a\alpha }k+\frac{c}{1+a\alpha }$ (10)

将寻求的最优初始值 ${\xi }_1$ 和 ${\xi }_2$ 带入模型中，再次以平均相对误差最小为目标(式(8))，通过粒子群优化

算法寻求最优的 $\alpha$ 、 $\beta$ 和 $\gamma$ 值，当 $\alpha =\beta =\lambda =\frac{1}{3}$ ，背景值变为杨孝良等 [16] 提出的紧邻三参数均值背景值

构造。

2.4. 优化三参数背景值NGM(1,1)模型的参数估计方法

传统的NGM(1,1)模型及其他的灰色预测模型都是利用最小二乘法求解模型的灰参数，即利用残差平方和最小为目标来求灰参数，即：

$s={\displaystyle \underset{k=3}{\overset{n}{\sum }}{\left(x^{(0)}\left(k\right)-{\hat{x}}^{(0)}\left(k\right)\right)}^2}$ (11)

式(11)没有考虑到新旧信息数据对于模型的贡献不同，对于序列长度为n的原始数据，根据数据值发生的时间的不同，给予不同重要程度，当数据值发生时间距离当前时刻越远，则其在建模时所占重要程度就越小，应该占有更小的比重；当数据值发生时间离当前时刻越近的数据，则在建模时的重要程度就越高，应该占有更大的比重。本文参考何霞等 [18] 用赋予残差平方和(式(11))权值w来反应建模序列值的重要程度，令：

$s={\displaystyle \underset{k=3}{\overset{n}{\sum }}w\left(k\right){\left(x^{(0)}\left(k\right)-{\hat{x}}^{(0)}\left(k\right)\right)}^2}$ (12)

其中 $w\left(k\right)=\frac{1}{{\left(z^{(1)}\left(k\right)\right)}^2}$ 。

将式(2)改为

$x^{(0)}\left(k\right)=-a{z}^{(1)}\left(k\right)+bk+c$ (13)

代入式(11)得：

$s={\displaystyle \underset{k=3}{\overset{n}{\sum }}w\left(k\right){\left(x^{(0)}\left(k\right)+a{z}^{(1)}\left(k\right)-bk-c\right)}^2}$ (14)

令 $y=\left[\begin{array}{l}w\left(3\right)\cdot x^{(0)}\left(3\right)\\ w\left(4\right)\cdot x^{(0)}\left(4\right)\\ ⋮\\ w\left(n\right)\cdot x^{(0)}\left(n\right)\end{array}\right]$ ， $B=\left[\begin{array}{l}w\left(3\right)\cdot z^{(1)}\left(3\right)3w\left(3\right)w\left(3\right)\\ w\left(4\right)\cdot z^{(1)}\left(4\right)4w\left(4\right)w\left(4\right)\\ ⋮⋮⋮\\ w\left(n\right)\cdot z^{(1)}\left(n\right)nw\left(n\right)w\left(n\right)\end{array}\right]$ ， $d={\left[a,b,c\right]}^{\text{T}}$

使s最小的灰参数a、b和c应该满足：

$\{\begin{array}{l}\frac{\partial s}{\partial a}=2{\displaystyle \underset{k=3}{\overset{n}{\sum }}w\left(k\right)\left(x^{(0)}\left(k\right)+a{z}^{(1)}\left(k\right)-bk-c\right)\cdot w\left(k\right)\cdot z^{(1)}\left(k\right)=0}\\ \frac{\partial s}{\partial b}=2{\displaystyle \underset{k=3}{\overset{n}{\sum }}w\left(k\right)\left(x^{(0)}\left(k\right)+a{z}^{(1)}\left(k\right)-bk-c\right)\cdot \left(-k\right)\cdot w\left(k\right)=0}\\ \frac{\partial s}{\partial c}=2{\displaystyle \underset{k=3}{\overset{n}{\sum }}w\left(k\right)\left(x^{(0)}\left(k\right)+a{z}^{(1)}\left(k\right)-bk-c\right)\cdot \left(-1\right)\cdot w\left(k\right)=0}\end{array}$ (15)

根据式(15)解得 $d={\left[a,b,c\right]}^{\text{T}}=\left(\begin{array}{l}\frac{\left(CE-ZF\right)H+\left(ZE-\frac{C{E}^2}{F}\right)I}{2ZEC-\frac{C^2{E}^2}{F}-Z^{2}F}\\ \frac{\left(CE-ZF\right)G+\left(AF-C^2\right)H+\left(ZC-AE\right)I}{2ZEC-\frac{C^2{E}^2}{F}-Z^{2}F}\\ \frac{\left(ZE-\frac{C{E}^2}{F}\right)G+\left(ZC-AE\right)H+\left(\frac{A{E}^2}{F}-Z^2\right)I}{2ZEC-\frac{C^2{E}^2}{F}-Z^{2}F}\end{array}\right)$

其中：

$A={\displaystyle \underset{k=3}{\overset{N}{\sum }}{w}^2\left(k\right){\left(z^{(1)}\left(k\right)\right)}^2}$ ， $Z={\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right)z^{(1)}\left(k\right)k}$ ， $C={\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right)z^{(1)}\left(k\right)}$ ， $E={\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right)k}$ ， $F={\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right)}$ ，

$G={\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right)z^{(1)}\left(k\right)x^{(0)}\left(k\right)}$ ， $H={\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right)x^{(0)}\left(k\right)}$ ， $I={\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right)x^{(0)}(\; k\; )}$

$\begin{array}{c}{B}^{\text{T}}B={\left(\begin{array}{l}w\left(3\right)z^{(1)}\left(3\right)3w\left(3\right)w\left(3\right)\\ w\left(4\right)z^{(1)}\left(4\right)4w\left(4\right)w\left(4\right)\\ ⋮⋮⋮\\ w\left(n\right)z^{(1)}\left(n\right)nw\left(n\right)w\left(n\right)\end{array}\right)}^{\text{T}}\left(\begin{array}{l}w\left(3\right)z^{(1)}\left(3\right)3w\left(3\right)w\left(3\right)\\ w\left(4\right)z^{(1)}\left(4\right)4w\left(4\right)w\left(4\right)\\ ⋮⋮⋮\\ w\left(n\right)z^{(1)}\left(n\right)nw\left(n\right)w\left(n\right)\end{array}\right)\\ =\left(\begin{array}{l}{\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right){\left(z^{(1)}\left(k\right)\right)}^2{\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right)z^{(1)}\left(k\right)k{\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right)z^{(1)}\left(k\right)}}}\\ {\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right)z^{(1)}\left(k\right)k}{\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right)k^2}{\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right)k}\\ ⋮⋮⋮\\ {\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right)z^{(1)}\left(k\right)}{\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right)k}{\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right)}\end{array}\right)=\left(\begin{array}{l}AZC\\ Z\frac{E^2}{F}E\\ CEF\end{array}\right)\end{array}$

$B^{\text{T}}Y={\left(\begin{array}{l}w\left(3\right)z^{(1)}\left(3\right)3w\left(3\right)w\left(3\right)\\ w\left(4\right)z^{(1)}\left(4\right)4w\left(4\right)w\left(4\right)\\ ⋮⋮⋮\\ w\left(n\right)z^{(1)}\left(n\right)nw\left(n\right)w\left(n\right)\end{array}\right)}^{\text{T}}\left(\begin{array}{l}w\left(3\right)x^{(0)}\left(3\right)\\ w\left(4\right)x^{(0)}\left(4\right)\\ ⋮\\ w\left(n\right)x^{(0)}\left(n\right)\end{array}\right)=\left(\begin{array}{l}{\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right)z^{(1)}\left(k\right)x^{(0)}\left(k\right)}\\ {\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right)x^{(0)}\left(k\right)k}\\ {\displaystyle \underset{k=3}{\overset{n}{\sum }}{w}^2\left(k\right)x^{(0)}\left(k\right)}\end{array}\right)=\left(\begin{array}{l}G\\ H\\ I\end{array}\right)$

$\left|B^{\text{T}}B\right|=2ZEC-\frac{C^2{E}^2}{F}-Z^{2}F=\theta$

由克拉默法则求得：

${\left(B^{\text{T}}B\right)}^{-1}=\left(\begin{array}{l}0\frac{CE-ZF}{\theta }\frac{ZE-\frac{C{E}^2}{F}}{\theta }\\ \frac{CE-ZF}{\theta }\frac{AF-C^2}{\theta }\frac{ZC-AE}{\theta }\\ \frac{ZE-\frac{C{E}^2}{F}}{\theta }\frac{ZC-AE}{\theta }\frac{\frac{A{E}^2}{F}-Z^2}{\theta }\end{array}\right)$

则：

$\begin{array}{c}{\left(B^{\text{T}}B\right)}^{-1}{B}^{\text{T}}Y=\left(\begin{array}{l}0\frac{CE-ZF}{\theta }\frac{ZE-\frac{C{E}^2}{F}}{\theta }\\ \frac{CE-ZF}{\theta }\frac{AF-C^2}{\theta }\frac{ZC-AE}{\theta }\\ \frac{ZE-\frac{C{E}^2}{F}}{\theta }\frac{ZC-AE}{\theta }\frac{\frac{A{E}^2}{F}-Z^2}{\theta }\end{array}\right)\left(\begin{array}{l}G\\ H\\ I\end{array}\right)\\ =\left(\begin{array}{l}\frac{\left(CE-ZF\right)H+\left(ZE-\frac{C{E}^2}{F}\right)I}{\theta }\\ \frac{\left(CE-ZF\right)G+\left(AF-C^2\right)H+\left(ZC-AE\right)I}{\theta }\\ \frac{\left(ZE-\frac{C{E}^2}{F}\right)G+\left(ZC-AE\right)H+\left(\frac{A{E}^2}{F}-Z^2\right)I}{\theta }\end{array}\right)={\left(a,b,c\right)}^{\text{T}}\end{array}$

所以利用加权最小二乘法求解模型参数依然满足：

$d={\left[a,b,c\right]}^{\text{T}}={\left(B^{\text{T}}B\right)}^{-1}{B}^{T}Y$ (16)

3. NGM(1,1)模型应用

首先，本节利用的数据实例来验证本文构建的组合优化模型的可行性，然后，通过与其他模型进行对比得出本文构建的组合优化模型的优化性，最后将本文构建的组合优化模型应用于武汉市2020年2月6日~2月15日的新增感染人数的拟合与预测，来检验本文构建模型的适用性。

3.1. 算例分析

通过文献 [16] 中的算例分析数据，分别建立GM(1,1)模型(模型1)、DGM(1,1)模型(模型2)、文献 [14] 改进初始值模型(模型3)及本文建立背景值、初始值和参数估计方式组合优化模型(模型4)，并且计算模型的相对误差(δ)、平均相误差(avg_δ)、方差比(C)和小误差概率(P)来比较模型的拟合性能和精度。

Table 1. Raw data

表1. 原始数据

根据表1中的原始数据，计算灰色预测模型模拟值、相对误差和平均相对误差如表2所示。其中相对误差(δ)和平均相误差(avg_δ)的计算公式分别为：

$\delta \left(k\right)=\left|\frac{x^{(0)}\left(k\right)-{\hat{x}}^{(0)}\left(k\right)}{x^{(0)}\left(k\right)}\right|\times 100\%$ (17)

$avg\_\delta =\frac{{\displaystyle \underset{k=1}{\overset{n}{\sum }}\delta \left(k\right)}}{n-1}$ (18)

通过寻优得到最优初始值为 ${\xi }_1=21.5204$ ， ${\xi }_1=44.3504$ ，代入式(10)，通过寻优得到的最优的三参数背景值加权权数为 $\alpha =0.2365$ ， $\beta =0.7635$ ， $\gamma =0.1744\times {10}^{-10}$ ，最后通过加权最小二乘法(式(16))计算得到灰参数 $a=-\text{0}\text{.3904}$ ， $b=-\text{8}\text{.8525}$ ， $c=\text{30}\text{.7150}$ ，从而通过计算得参数 $a_1=1.4302$ ， $a_2=7.0521\times {10}^{-12}$ ， $a_3=-9.7532$ ， $a_4=33.8402$ ，根据式(10)可知，表1数据序列的NGM(1,1)模型为：

$x^{(1)}\left(k\right)=\text{1}\text{.4302}{x}^{(1)}\left(k-1\right)+\text{7}\text{.5021}\times {\text{10}}^{-\text{12}}{x}^{(1)}\left(k-2\right)-\text{9}\text{.7532}k+\text{33}\text{.8402}$ (19)

其中 $k=3,4,\cdots ,n$ 。同时计算GM(1,1)模型(模型1)、DGM(1,1)模型(模型2)、文献 [14] 改进初始值模型(模型3)及本文建立背景值、初始值和参数估计方式组合优化模型(模型4)得模拟值、相对误差和平均相对误差来比较模型的模拟性能，如表2所示。

Table 2. Average error and average relative error of simulated values

表2. 模拟值的平均误差和平均相对误差

为了进一步比较模型的模拟精度，计算出四个模型的方差比(C)和小误差概率(P)，首先计算出四个模型得出的模拟值与原始序列(表1)的残差值，残差计算公式：

${\epsilon }_i\left(k\right)=x^{(0)}\left(k\right)-{\hat{x}}^{(0)}\left(k\right)$ (20)

计算得到残差序列如表3所示：

Table 3. Residual

表3. 残差

通过表3计算得到的残差序列值，进一步计算原始序列的均值( ${\bar{x}}^{(0)}$ )和标准差( $s_1$ )以及三个模型与原始序列的残差序列均值( ${\bar{\epsilon }}_i$ )与标准差( $s_{1i}$ )、方差比(C)和小误差概率(P)，计算公式如下：

${\bar{x}}^{(0)}=\frac{1}{n}{\displaystyle \underset{k=1}{\overset{n}{\sum }}{x}^{(0)}\left(k\right)}$ (21)

$s_1=\sqrt{\frac{{\displaystyle \underset{k=1}{\overset{n}{\sum }}\left(x^{(0)}\left(k\right)-{\bar{x}}^{(0)}\right)}}{n}}$ (22)

${\bar{\epsilon }}_i=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{4}{\sum }}{\epsilon }_i}$ (23)

$s_{1i}=\sqrt{\frac{{\displaystyle \underset{i=1}{\overset{4}{\sum }}{\epsilon }_{1i}}}{n}}$ (24)

$C_i=\frac{s_1}{s_{1i}}$ (25)

$P_i=p\left(\left|{\epsilon }_i\left(t\right)-{\bar{\epsilon }}_i\right|<0.6745s_1\right)$ (26)

式(20)~(21)计算表1原始数据的均值为27.91，标准差为5.76，通过式(22)~(26)计算四个模型的平均相对残差(q)、方差比(C)和小误差概率(P)如表4所示。小误差概率和方差比对于模型预测精度的评价准则如下：

Table 4. Accuracy evaluation criteria

表4. 精度评价准则

Table 5. Mean relative error, varianceratio, and small error probability test

表5. 平均相对残差、方差比和小误差概率检验

通过表2、表5和表4的评价准则比较4个模型相关指标得出每个模型的性能，从表2可以看出本文所构建初始值、背景值和参数估计方法组合优化的模型(模型4)的平均相对误差为3.021%，都小于GM(1,1)模型(模型1)、将作为初始值的改进GM(1,1)模型(模型2)和文献 [14] 改进初始值模型(模型3)，通过表4的预测精度评价准则以及表5计算得出的指标值比较4个模型的方差比(C)和小误差概率(P)，得出本文所构建模型(模型4)的模拟精度等级、将作为初始值的改进GM(1,1)模型(模型2)和文献 [14] 改进初始值模型(模型3)的模拟精度等级最高，且本文所构建模型(模型4)的模拟精度最高，GM(1,1)模型(模型1)的模拟精度等级为2最差。通过该算例得出本文将初始值通过文献 [16] 将初始值设为未知变量 ${\xi }_1$ 和 ${\xi }_2$ ，以平均相对误差最小为目标利用粒子群优化算法求最优初始值的取值，将三参数的紧邻均值背景值改为加权的三参数背景值，利用粒子群优化算法寻求最优的加权权数，最后将参数估计方法最小二乘法改为加权最小二乘法求模型的灰参数，构成NGM(1,1)模型初始值的选取、背景值的构造和参数估计方法的组合优化一定程度上使得模型精度提高，因此本文所构建组合优化模型优于其他三种模型。

3.2. 模型应用

下面通过将模型应用于武汉市新型冠状病毒肺炎新增感染人数的预测来分析本文的组合优化模型的效果，将武汉市2020年2月6日~2020年2月15日(封城阶段，假设不存在人口流动)作为原始数据，比较传统的GM(1,1)模型、文献 [16] 模型、文献 [14] 模型和本文组合优化模型平均相对误差来比较模型的预测精度，如表6所示。通过寻优得到最优初始值为 ${\xi }_1=1686$ ， ${\xi }_2=3671$ ，代入式(10)，通过寻优得到的最优的三参数背景值加权权数为 $\alpha =0.4910$ ， $\beta =\text{6}\text{.1000}\times {\text{10}}^{-\text{14}}$ ， $\gamma =0.5090$ ，最后通过加权最小二乘法(式(16))计算得到灰参数 $a=-\text{0}\text{.8661}$ ， $b=-\text{1}\text{.3818}\times {\text{10}}^{\text{3}}$ ， $c=\text{2816}\text{.5495}$ ，从而通过计算得参数 $a_1=\text{1}\text{.7399}$ ， $a_2=\text{0}\text{.7670}$ ， $a_3=-\text{2}\text{.4042}\times {\text{10}}^{\text{3}}$ ， $a_4=\text{4}\text{.9004}\times {\text{10}}^{\text{3}}$ ，根据式(10)可知，表1数据序列的NGM(1,1)模型为：

$x^{(1)}\left(k\right)=\text{1}\text{.7399}{x}^{(1)}\left(k-1\right)+\text{0}\text{.7670}{x}^{(1)}\left(k-2\right)-\text{2}\text{.4042}\times {\text{10}}^{\text{3}}k+\text{4}\text{.9004}\times {\text{10}}^{\text{3}}$ (27)

其中 $k=3,4,\cdots ,n$ ，通过式(27)计算出本文组合优化模型的 $x^{(1)}\left(k\right)$ 的模拟值，之后根据 $x^{(1)}\left(k\right)$ 、 $x^{(1)}\left(k-1\right)$ 与 $x^{(0)}\left(k\right)$ 的关系，计算出原始数据序列的模拟值，根据模拟值与原始数据本文组合优化模型的相对误差和平均相对误差。同时计算GM(1,1)模型、文献 [16] 模型、文献 [14] 模型的模拟值、相对误差和平均相对误差。

从表6可以看出本文组合优化模型的平均相对误差小于GM(1,1)模型，且远小于文献 [16] 和文献 [14] 模型的平均相对误差，且本文基于的文献 [16] 的三参数紧邻均值背景值NGM(1,1)模型平均相对误差最大，这也说明本文在此基础上对初始值、背景值和参数估计的组合优化的有效性。为了更加直观的比较GM(1,1)模型、文献 [16] 、文献 [14] 和本文组合优化模型对于武汉市2020年2月6日~2020年2月15日的新增感染肺炎人数的模拟情况，本文绘制GM(1,1)模型、文献 [16] 、文献 [14] 和本文组合优化模型的相对误差折线图，如图1所示：

Table 6. Average error and average relative error

表6. 平均误差和平均相对误差

Figure 1. Relative error line chart

图1. 相对误差折线图

从图1可以更加清楚地看出本文构建的组合优化模型优于GM(1,1)模型，且远优于文献 [16] 和文献 [14] 模型。这一结果表明本文构建的组合优化模型具有合理性和有效性。

4. 结论

本文构建对于近似非齐次指数数据模拟效果好的NGM(1,1)预测模型，考虑到灰色模型的背景值容易受到极端值的影响，所以为了减少原始数据中极端值对于灰色预测模型模拟性能的影响，达到提高灰色预测模型模拟性能稳定性的目的，本文参考文献 [16] 构建NGM(1,1)模型的三参数紧邻均值背景值，基于此，依次对三参数紧邻均值背景值的NGM(1,1)模型的初始值、背景值和参数估计方式进行优化。首先，以平均相对误差最小为目标利用粒子群优化算法求最优初始值的取值，然后将三参数的紧邻均值背景值改为加权的三参数背景值，利用粒子群优化算法寻求最优的加权权数，最后将参数估计方法最小二乘法改为加权最小二乘法求模型的灰参数，实现NGM(1,1)模型初始值、背景值和参数估计方式的协同优化。通过算例分析，验证了本文构建的初始值、背景值和参数估计方式组合优化模型模拟性能优于经典GM(1,1)模型、DGM(1,1)模型和文献 [14] 的初始值优化模型。最后，将本文的组合优化模型应用于武汉市2020年2月6日~2020年2月15日新增感染人数的模拟，模拟值的平均相对误差为0.274低于其余三个灰色预测模型的平均相对误差，再次验证了本文构建的组合优化模型的有效性和合理性。

基金项目

基金项目：国家自然科学基金青年项目(71701066)；上海理工大学管理学院2021年度科研启动项目(KYQD202101)资助课题。

参考文献

NOTES

^*通讯作者。