1.1GM(1,1)模型

1.1.1基本概念

• 序列级比：

  已知数据序列

  $$ X=(x(1),x(2),\cdots,x(n))，x(k)\ge0,k=1,2,\cdots,n $$

  则称

  $$ \sigma(k)=\frac{x(k)}{x(k-1)},k=2,3,\cdots,n $$

  为序列X的级比。

• 一次累加生成算子—1-AGO：

  已知原始数据序列

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

  令$D$为序列算子，则

  $$ X^{(0)}D=(x^{(0)}(1)d,x^{(0)}(2)d,\cdots,x^{(0)}(n)d) $$

  其中

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

  此时称$D$为$X^{(0)}​$的一次累加生成算子，记为1-AGO。

• 紧邻均值生成序列：

  已知原始数据序列

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

  其中$x^{(0)}(k)$和$x^{(0)}(k-1)$为$X^{(0)}$的一对紧邻值，则可定义序列算子D为:

  $$ X^{(0)}D=(x^{(0)}(1)d,x^{(0)}(2)d,\cdots,x^{(0)}(n)d) $$

  其中

  $$ x^{(0)}(k)d=\alpha x^{(0)}(k)+(1-\alpha)x^{(0)}(k-1) $$

  当$\alpha=0.5$时，无偏生成算子D也称为紧邻均值生成算子,紧邻生成算子作用的序列称为紧邻均值生成序列。

• 关联度（用于模型的关联度检验）：设系统行为序列有m个，即：

$$ \begin{aligned} X^{(0)}&=(x^{(0)}(1),x^{(0)}(2),\cdots,x^{(0)}(n)),\\ X^{(1)}&=(x^{(1)}(1),x^{(1)}(2),\cdots,x^{(1)}(n)),\\ &\qquad\qquad\qquad\qquad\ \ \ \vdots\\ X^{(m)}&=(x^{(m)}(1),x^{(m)}(2),\cdots,x^{(m)}(n)) \end{aligned} $$

则称

$$ \begin{aligned} \gamma(x^{(0)}(k),x^{(i)}(k))&=\frac{min_imin_k|x^{(0)}(k)-x^{(i)}(k)|+Pmax_imax_k|x^{(0)}(k)-x^{(i)}(k)|}{|x^{(0)}(k)-x^{(i)}(k)|+Pmax_imax_k|x^{(0)}(k)-x^{(i)}(k)|}\\ \Gamma(X^{(0)},X^{(i)})&=\frac1n\sum_{k=1}^{n}{\gamma(x^{(0)}(k).x^{(i)}(k))} \end{aligned} $$

其中$i=1,2,\cdots,m；k=1,2,\cdots,n​$,$\gamma(x^{(0)}(k),x^{(i)}(k))​$表示$x^{(i)}(k)​$与$x^{(0)}(k)​$在$k​$时刻的灰色关联系数，$P\in[0,1]​$为分辨系数，$\Gamma(X^{(0)},X^{(i)})​$表示比较数列$X^{(i)}​$对参考数列$X^{(0)}​$的灰色关联度。

==灰色关联度的计算步骤：==

• [ ] 第一步：计算各序列的均值象${X^{(i)}}{'}​$

  $$ {X^{(i)}}{'}=\frac{X^{(i)}}{x^{(i)}(1)}=(x^{(i)}{'}(1),x^{(i)}{'}(2),\cdots,x^{(i)}{'}(n)),i=0,1,2,\cdots,m $$

• [ ] 第二步：计算$X^{(0)}​$与$X^{(i)}​$均值象对应分量之差的绝对值序列，记为$\Delta^{(i)}(k)​$,即

  $$ \Delta^{(i)}(k)=|x^{(0)}{'}(k)-x^{(i)}{'}(k)|,i=1,2,\cdots,n;k=1,2,\cdots,m $$

• [ ] 第三步：找到$\Delta^{(i)}(k)$的最大值与最小值，也称为两级最大差和两级最小差，分别记为

  $$ M=\max_i\max_k\Delta^{(i)}(k),L=\min_i\min_k\Delta^{(i)}(k) $$

• [ ] 第四步：计算关联系数

  $$ \gamma(x^{(0)}(k),x^{(i)}(k))=\frac{L+PM}{\Delta^{(i)}(k)+PM},P\in[0,1];i=1,2,\cdots,m;k=1,2,\cdots,n $$

• [ ] 第五步：计算关联度——即关联系数的平均值

  $$ \Gamma(X^{(0)},X^{(i)})=\frac1n\sum_{k=1}^{n}{\gamma(x^{(0)}(k).x^{(i)}(k))},i=1,2,\cdots,m;k=1,2,\cdots,n $$

1.1.2 建模过程

GM(1,1)模型有四种基本形式定义即，原始差分GM(1,1) 模型(ODGM)、均值GM(1,1)模型(EGM)、均值差分GM(1,1)模型(EDGM)、离散GM(1,1)模型(DGM)，其中均值GM(1,1)模型应用更加广泛，一般提到GM(1,1)模型往往指的是EGM(1,1)。下面对均值GM(1,1)模型建模过程进行阐述:

i)数据检验和处理——构建GM(1,1)模型的均值方程

数据检验：对已知原始数据计算所有级比$\sigma(k)$,当所有级比$\sigma(k)\in(e^{\frac{-2}{n+1}},e^{\frac{2}{n+2}})$时，则原始数据可以作为GM(1,1)模型的数据进行灰色预测。否则需要对原始数据进行平移变换使其满足条件。

数据处理：已知原始数据序列$X^{(0)}=(x^{(0)}(1),x^{(0)}(2),\cdots,x^{(0)}(n))$

• 首先对$X^{(0)}​$进行1-AGO得到序列

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

  其中

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

• 然后由$X^{(0)}$和$X^{(1)}$得到$\alpha=\frac{1}{2}$时，邻值均值生成数列为

  $$ Z^{(0)}=(z^{(1)}(2),z^{(1)}(3),\cdots,z^{(1)}(n), $$

  其中

$$ z^{(1)}(k)=\alpha[{x^{(1)}(k)+x^{(1)}(k-1)}]=\frac{1}{2}[{x^{(1)}(k)+x^{(1)}(k-1)}] $$

• 将

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

  代入上面式(1)可知：

$$ z^{(1)}(k)=\sum_{i=1}^{k-1}{x^{(0)}(i)}+\frac{1}{2}x^{(0)}(k) $$

• 由上式可知$z^{(1)}(k)​$与$x^{(0)}(k)​$存在线性关系，为了方便计算简化方程形式，我们对$z^{(1)}(k)​$与$x^{(0)}(k)​$的线性关系构建如下：

$$ x^{(0)}(k)+az^{(1)}(k)=b $$

• 上面式（3）即为GM(1,1)模型的均值形式，==GM(1,1)模型均值形式实质上是一个差分方程==，a称为发展系数，b称为灰作用量，其中上述方程中$x^{(0)}(k)和z^{(1)}(k)$为变量，接下来需要对a和b进行求解

ii)矩阵建立B、Y,求解a、b——实质是利用最小二乘法进行参数估计得出$\hat a,\hat b​$

• 根据GM(1,1)模型均值形式方程，记$\hat a=[a,b]^T​$，记

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

• 则GM(1,1)模型均值形式方程可以表示为$Y=B\hat a​$

• 利用最小二乘法对a和b进行参数估计,即运用最小二乘法估计$\hat a$,则

  $$ \hat a=(B^TB)^{-1}B^TYiii) $$

iii)求解GM(1,1)模型均值形式的时间响应序列式——类比微分方程求解方法——构造该模型的白化微分方程

• 把均值GM(1,1)模型的差分形式连续化，即把$k=2,3,\cdots,n$视为连续变量t,则

$$ \begin{aligned} x^{(0)}(k)&=x^{(1)}(k)-x^{(1)}(k-1)\\ &=\int_{k-1}^{k}d{{x^{(1)}(t)}}\\ &=\int_{k-1}^{k}\frac{d{{x^{(1)}(t)}}}{dt}{dt} \end{aligned} $$

• 同时根据积分的几何意义对$z^{(1)}(k)=\frac{1}{2}[{x^{(1)}(k)+x^{(1)}(k-1)}]​$为在点$t=k和t=k-1与曲线x^{(1)}(t)​$所围成的矩形面积近似转化为积分形式，即

$$ \begin{aligned} z^{(1)}(k)&=\frac{1}{2}[{x^{(1)}(k)+x^{(1)}(k-1)}\\ &\approx\int_{k-1}^{k}{x^{(1)}(t)}dt \end{aligned} $$

• 把上面$x^{(0)}(k)和z^{(1)}(k)​$代入均值GM(1,1)模型方程式(3)可得：

$$ x^{(0)}(k)=-az^{(1)}(k)+b $$

​ 即:

$$ \begin{aligned} &\int_{k-1}^{k}\frac{d{{x^{(1)}(t)}}}{dt}{dt}\\ &\quad=-a\int_{k-1}^{k}{x^{(1)}(t)}dt+b\\ &\quad=-\int_{k-1}^{k}a({x^{(1)}(t)}-\frac{b}{a})dt \end{aligned} $$

• 对上式移项后左右两边同时求导可得：

$$ \frac{d{{x^{(1)}(t)}}}{dt}{dt}=-ax^{(1)}(t)+b $$

• 故上面白化微分方程式(9)即为均值GM(1,1)模型的影子方程，==由常微分方程求解方法==可得式(9)的通解为：

$$ \begin{aligned} x^{(1)}(t)&=e^{-{\int adt}}[\int {be^{\int adt}}dt+c]\\ &=e^{(-at)}[\frac{b}{a}e^{(at)}+c]\\ &=\frac ba+ce^{(-at)} \end{aligned} $$

• 由上述可知初值条件$x^{(1)}(1)=x^{(0)}(1)$代入上式通解即$x^{(1)}(1)=x^{(0)}(1)=\frac ba+ce^{-a}$,由此可得$c=e^{a}[x^{(0)}(1)-\frac ba]$故白化微分方程的解为：

$$ x^{(1)}(t)=\frac ba+[x^{(0)}(1)-\frac ba]e^{-a(t-1)} $$

对于一阶微分方程形如

$$ y^{'}=-P(x)y+Q(x) $$

可以对该微分方程利用常数变异法可得该微分方程的通解

$$ y=\Big[\int Q(x)e^{\int P(x)dx}dx+C\Big]e^{-\int P(x)dx} $$

iv)根据白化微分方程（影子方程）的解——得出GM(1,1)模型时间响应序列式——GM(1,1)模型时间响应序列式即为影子方程解的离散形式

$$ \hat x^{(1)}(k)=\frac{b}{a}+\Big[x^{(0)}(1)-\frac{b}{a}\Big]e^{-a(k-1)},k=1,2,\cdots,n $$

v)累减还原——得出预测值

$$ \begin{aligned} \hat x^{(0)}(k)&=\hat x^{(1)}(k)-\hat x^{(1)}(k-1)\\ &=(1-e^a)\Big[x^{(0)}(1)-\frac ba\Big]e^{-a(k-1)}，\\ \text{其中}k&=1,2,\cdots,n \end{aligned} $$

1.1.3模型检验

1.1.3.1残差检验

残差大小检验，即对模型值和实际值的残差进行逐点检验

• 绝对残差序列$\epsilon^{(0)}$计算：首先按模型计算$\hat x^{(1)}(k+1)$,然后通过累减还原得到$\hat x^{(0)}(k)$,最后计算原始序列$x^{(0)}(k)$和模型计算值$\hat x^{(0)}(k)$的绝对残差序列:

  $$ \epsilon^{(0)}={\{\epsilon^{(0)}(k),k=1,2,\cdots,n}\},其中\epsilon^{(0)}(k)=|x^{(0)}(k)-\hat x^{(0)}(k)| $$

• 相对残差序列$\Delta^{(0)}$计算:

  $$ \begin{aligned} \Delta^{(0)}&={\{\Delta^{(0)}(k)},k=1,2,\cdots,n\}\\ &={\{\frac{\epsilon^{(0)}(k)}{x^{(0)}(k)}},k=1,2,\cdots,n\} \end{aligned} $$

• 平均相对残差$\bar \Delta$计算：

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

给定$\alpha$,当$\bar\Delta\lt\alpha$,且$\Delta\lt\alpha$成立时，称模型为残差合格。注：$\alpha常取0.01、0.05、0.1、0.2$

1.1.3.2关联度检验

$\Gamma(X^{(0)},\hat X^{(0)})$表示模型计算数列$\hat X^{(0)}$对原始数列$X^{(0)}$的灰色关联度，对于给定的$\Gamma_0(\Gamma_0常取0.9、0.8、0.7、0.6)，有\Gamma(X^{(0)},\hat X^{(0)})\gt\Gamma_0$，则称模型为关联度合格模型​。

1.1.3.3后验差检验（残差检验的基础上）

后验差检验，及对残差分布的统计特性进行验证。

• 计算原始序列平均值$\bar x^{(0)}$

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

• 计算原始序列均方差$S_1$

  $$ S_1=\sqrt{\frac{\sum_{k=1}^n(|x^{(0)}(k)-\bar x^{(0)}|)^2}{n-1}} $$

• 计算残差的均值$\bar\Delta$

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

• 计算残差的均方差$S_2$

  $$ S_2=\sqrt{\frac{\sum_{k=1}^n(|\Delta^{(0)}(k)-\bar \Delta|)^2}{n-1}} $$

• 计算方差比$C$

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

• 计算小残差概率$P$

  $$ P=P{\{|\Delta^{(0)}(k)-\bar \Delta|\lt0.6745S_1\}} $$

  若对于给定的$C_0,(C_0常取0.35、0.50、0.65、0.80)，当C\lt C_0$时，称模型为均方差比合格模型；对于给定的$P_0,(P_0常取0.95、0.80、0.70、0.60)，当P\gt P_0$时，称模型为小残差概率合格模型。

  ==若残差检验、关联度检验、后验差检验在允许范围内，则可以利用所建模型进行预测，否则需要进行残差修正——GM(1,1)残差模型。==

  1.1.3.4小结

  上述检验模型的三种方法都是通过对残差的考察来判断模型的精度,其中平均相对误差$\bar \Delta $和模拟误差都要求越小越好,关联度$\Gamma$要求越大越好，均方差比值$C$越小越好(因为$C$小说明$S_2$小，$S_1$大,即残差方差小，原始数据方差大,说明残差比较集中,摆动幅度小，原始数据比较分散,摆动幅度大,所以模拟效果好要求$S_2$与$S_1$,相比尽可能小)，以及小误差概率P越大越好,给定$\alpha, \Gamma_0, C_0, P_0$的一组取值,就确定了检验模型模拟精度的一个等级.常用的精度等级见表1,可供检验模型参考。

  精度等级	相对误差$\alpha$	关联度$\Gamma_0$	均方差比值$C_0$	小残差概率$P_0$
  一级	0.01	0.90	0.35	0.95
  二级	0.05	0.80	0.50	0.80
  三级	0.10	0.70	0.65	0.70
  四级	0.20	0.60	0.80	0.60

1.1.5模型说明

1）给定原始数据$X^{(0)}$中的数据不一定要全部用来建模，对原始数据的取舍不同，得到的模型不同，即a,b不同；

2）建模数据的取舍应保证建模序列等时距、相连、不得有跳跃出现；

3）一般建模的数据序列应当由最新的数据以及相邻数据构成，当再出现新数据时有两种处理方法：i）将新信息加入到原始序列中，重新估参；ii)去掉原始数据中最老的一个数据，再加上最新的数据，所形成的序列和原序列维数相等，再重新估参数。