作者：張張張張
github地址：https://github.com/zhanghekai
【转载请注明出处，谢谢！】

文章目录

一、灰色预测模型GM(1,1)
- 1、数据检验与数据预处理
- - 1.1 构建模型前检验
  - 1.2 数据生成
  - - (1) 累加生成（AGO）
    - (2)累减生成（IAGO）
    - (3)加权邻值生成
- 2 、灰色模型GM（1，1）构建
- - 2.1 定义$x^{(1)}$的灰导数为
  - 2. 2 GM(1,1)的白化模型（也叫影子模型）
- 3、检验预测值
- - 3.1 残差检验：计算相对残差
  - 3.2 级比偏差检验
二、使用灰色模型进行预测案例

一、灰色预测模型GM(1,1)

灰色预测是用灰色模型GM(1,1)来进行定量分析的，通常分为以下几类：
(1) 灰色时间序列预测。用等时距观测到的反映预测对象特征的一系列数量（如产量、销量、人口数量、存款数量、利率等）构造灰色预测模型，预测未来某一时刻的特征量，或者达到某特征量的时间。
(2) 畸变预测（灾变预测）。通过模型预测异常值出现的时刻，预测异常值什么时候出现在特定时区内。
(3) 波形预测，或称为拓扑预测，它是通过灰色模型预测事物未来变动的轨迹。
(4) 系统预测，对系统行为特征指标建立一族相互关联的灰色预测理论模型，在预测系统整体变化的同时，预测系统各个环节的变化。

灰色模型预测流程

数据检验与数据预处理
模型构建
模型检验

1、数据检验与数据预处理

1.1 构建模型前检验

\qquad给定观测数据列x(0)={x1(0),x2(0),⋯ ,xn(0)}x^{(0)}=\{x_1^{(0)},x_2^{(0)},\cdots,x_n^{(0)}\}x(0)={x1(0),x2(0),⋯,xn(0)}一般采用x(0)x^{(0)}x(0)的级比λk(0)\lambda^{(0)}_kλk(0)的大小与所属区间来判断。
\qquad级比定义：λk(0)=xk−1(0)xk(0),k=2,3,⋯ ,n\lambda^{(0)}_k=\frac{x^{(0)}_{k-1}}{x^{(0)}_k},\quad k=2,3,\cdots,nλk(0)=xk(0)xk−1(0),k=2,3,⋯,n
\qquad若满足λk(0)∈(e−2n+1,e2n+1)\lambda^{(0)}_k\in (e^{-\frac{2}{n+1}},e^{\frac{2}{n+1}})λk(0)∈(e−n+12,en+12)，则认为x(0)x^{(0)}x(0)是可作为GM(1,1)GM(1,1)GM(1,1)建模的。如果原始序列不满足，需要对原始序列做必要的变换处理，使其落入可容覆盖内。即取适当的常数CCC，做平移变换：
y(0)=xi(0)+C,i=1,2,⋯ ,ny^{(0)}=x^{(0)}_i+C,\quad i=1,2,\cdots,ny(0)=xi(0)+C,i=1,2,⋯,n
\qquad使序列y(0)=(y1(0),y2(0),⋯ ,yn(0))y^{(0)}=(y^{(0)}_1,y^{(0)}_2,\cdots,y^{(0)}_n)y(0)=(y1(0),y2(0),⋯,yn(0))的级比λy(k)=yk−1(0)y(k)\lambda_y^{(k)}=\frac{y^{(0)}_{k-1}}{y^{(k)}}λy(k)=y(k)yk−1(0)落入上述区间范围。

1.2 数据生成

(1) 累加生成（AGO）

\qquad设原始序列为x(0)=(x1(0),x2(0),⋯ ,xn(0))x^{(0)}=(x^{(0)}_1,x^{(0)}_2,\cdots,x^{(0)}_n)x(0)=(x1(0),x2(0),⋯,xn(0))，令：xk(1)=∑i=1kxi(0),k=1,2,⋯ ,nx^{(1)}_k=\sum_{i=1}^{k}x^{(0)}_i,\quad k=1,2,\cdots ,nxk(1)=i=1∑kxi(0),k=1,2,⋯,n有：x(1)=(x1(1),x2(1),⋯ ,xn(1))x^{(1)}=(x^{(1)}_1,x^{(1)}_2,\cdots,x^{(1)}_n)x(1)=(x1(1),x2(1),⋯,xn(1))，称所得到的新数列x(1)x^{(1)}x(1)为数列x(0)x^{(0)}x(0)的1次累加生成的数列。

(2)累减生成（IAGO）

\qquad通过累加数列得到的新数列，可以通过累减生成还原出原始数列，设原始数列x(1)=(x1(1),x2(1),⋯ ,xn(1))x^{(1)}=(x^{(1)}_1,x^{(1)}_2,\cdots,x^{(1)}_n)x(1)=(x1(1),x2(1),⋯,xn(1))，令：
xk(0)=(xk(1)−xk−1(1)),k=2,3,⋯nx^{(0)}_k=(x^{(1)}_k-x^{(1)}_{k-1}),\quad k=2,3,\cdots nxk(0)=(xk(1)−xk−1(1)),k=2,3,⋯n称所得到的数列x(0)x^{(0)}x(0)为x(1)x^{(1)}x(1)的1次累减生成数列。

(3)加权邻值生成

\qquad令z(1)z^{(1)}z(1)为x(1)x^{(1)}x(1)的紧邻均值生成序列：
z(1)=(z2(1),z3(1),⋯ ,zn(1))z^{(1)}=(z^{(1)}_2,z^{(1)}_3,\cdots,z^{(1)}_n)z(1)=(z2(1),z3(1),⋯,zn(1))其中，
zk(1)=αxk(1)+(1−α)xk−1(1),k=2,3,⋯ ,nz^{(1)}_k=\alpha x^{(1)}_k+(1-\alpha)x^{(1)}_{k-1},\quad k=2,3,\cdots,nzk(1)=αxk(1)+(1−α)xk−1(1),k=2,3,⋯,nα\alphaα称为生成系数，当α=0.5\alpha = 0.5α=0.5时，则称该数列为均值生成数，也称为等权邻值生成数。

2 、灰色模型GM（1，1）构建

2.1 定义x(1)x^{(1)}x(1)的灰导数为

dk=xk(0)=xk(1)−xk−1(1)d_k=x^{(0)}_k=x^{(1)}_k-x^{(1)}_{k-1}dk=xk(0)=xk(1)−xk−1(1) 于是GM(1,1)的灰微分方程模型为：
dk+azk(1)=b或xk(0)+azk(1)=b(1)d_k+a z^{(1)}_k = b \quad或 \quad x^{(0)}_k+a z^{(1)}_k = b\qquad\qquad\qquad(1)dk+azk(1)=b或xk(0)+azk(1)=b(1)
其中，xk(0)x^{(0)}_kxk(0)称为灰导数，α\alphaα称为发展系数，zk(1)z^{(1)}_kzk(1)称为白化背景值，bbb称为灰作用量。将时刻k=2,3,⋯ ,nk=2,3,\cdots,nk=2,3,⋯,n代入上式，有：
{x2(0)+az2(1)=bx3(0)+az3(1)=b⋯xn(0)+azn(1)=b\begin{cases}x^{(0)}_2+a z^{(1)}_2=b\\ x^{(0)}_3+a z^{(1)}_3=b\\ \qquad\quad\cdots\\ x^{(0)}_n+a z^{(1)}_n=b \end{cases}⎩⎪⎪⎪⎨⎪⎪⎪⎧x2(0)+az2(1)=bx3(0)+az3(1)=b⋯xn(0)+azn(1)=b
引入矩阵向量记号：
μ=[ab],Y=[x2(0)x3(0)⋮xn(0)],B=[−z2(1)1−z3(1)1⋮−zn(1)1]\mu=\begin{bmatrix}a\\b\end{bmatrix} ,Y=\begin{bmatrix}x^{(0)}_2\\x^{(0)}_3\\ \vdots\\x^{(0)}_n\end{bmatrix},B=\begin{bmatrix}-z^{(1)}_2\quad 1\\-z^{(1)}_3\quad 1\\\vdots\\-z^{(1)}_n\quad 1 \end{bmatrix}μ=[ab],Y=⎣⎢⎢⎢⎢⎡x2(0)x3(0)⋮xn(0)⎦⎥⎥⎥⎥⎤,B=⎣⎢⎢⎢⎢⎡−z2(1)1−z3(1)1⋮−zn(1)1⎦⎥⎥⎥⎥⎤
于是GM（1，1）模型可以表示为Y=BμY=B\muY=Bμ
那么现在的问题就是求aaa和bbb的值，可以用最小二乘法求它们的估计值为：
μ=[ab]=(BTB)−1BTY\mu=\begin{bmatrix}a\\b\end{bmatrix}=(B^TB)^{-1}B^TYμ=[ab]=(BTB)−1BTY

2. 2 GM(1,1)的白化模型（也叫影子模型）

\qquad对于GM(1,1)GM(1,1)GM(1,1)的灰微分方程。如果将时刻k=2,3,⋯ ,nk=2,3,\cdots,nk=2,3,⋯,n视为连续变量ttt，则之前的x(1)x^(1)x(1)视为时间ttt函数，于是灰导数x(0)x^{(0)}x(0)变为连续函数的导数dxt(1)dt\frac{dx^{(1)}_t}{dt}dtdxt(1),白化背景值zk(1)z^{(1)}_kzk(1)对应于导数xt(1)x^{(1)}_txt(1)。于是GM（1，1）的灰微分方程对应于的白微分方程为：
dxt(1)dt+axt(1)=b(2)\frac{dx^{(1)}_t}{dt}+ax^{(1)}_t=b\qquad\qquad\qquad(2)dtdxt(1)+axt(1)=b(2)
解为：
xt(1)=(x1(0)−ba)e−a(t−1)+bax^{(1)}_t=(x^{(0)}_1-\frac{b}{a})e^{-a(t-1)}+\frac{b}{a}xt(1)=(x1(0)−ab)e−a(t−1)+ab

注1：以往的文献中并没有对这种从离散形式到连续形式的直接跳跃做出解释。

于是得到预测值：
x^k+1(1)=(x1(0)−ba)e−ak+ba,k=1,2,⋯ ,n−1\hat{x}^{(1)}_{k+1}=(x^{(0)}_1-\frac{b}{a})e^{-ak}+\frac{b}{a},\qquad k=1,2,\cdots,n-1x^k+1(1)=(x1(0)−ab)e−ak+ab,k=1,2,⋯,n−1从而得到相应的预测值：x^k+1(0)=x^k+1(1)−x^k(1),k=1,2,⋯ ,n−1\hat{x}^{(0)}_{k+1}=\hat{x}^{(1)}_{k+1}-\hat{x}^{(1)}_k,\qquad k=1,2,\cdots,n-1x^k+1(0)=x^k+1(1)−x^k(1),k=1,2,⋯,n−1

注2：原始序列数据不一定要全部使用，相应建立的模型也会不同，即a和b不同。
注3：原始序列数据必须要等时间间隔、不间断。

3、检验预测值

3.1 残差检验：计算相对残差

ξk=xk(0)−x^k(0)xk(0),k=1,2,⋯ ,n\xi_k=\frac{x^{(0)}_k-\hat{x}^{(0)}_k}{x^{(0)}_k},\qquad k=1,2,\cdots,nξk=xk(0)xk(0)−x^k(0),k=1,2,⋯,n如果对所有的∣ξk∣<0.1|\xi_k|<0.1∣ξk∣<0.1，则认为达到较高的要求；
如果对所有的∣ξk∣<0.2|\xi_k|<0.2∣ξk∣<0.2，则认为达到一般的要求。

3.2 级比偏差检验

根据前面计算出来的级比λk\lambda_kλk和发展系数aaa，计算相应的级比偏差：
ρk=1−1−0.5a1+0.5aλk\rho_k=1-\frac{1-0.5a}{1+0.5a}\lambda_kρk=1−1+0.5a1−0.5aλk若∣ρk∣<0.1|\rho_k|<0.1∣ρk∣<0.1，则认为达到较高要求；
若∣ρk∣<0.2|\rho_k|<0.2∣ρk∣<0.2，则认为达到一般要求；

二、使用灰色模型进行预测案例

北方某城市1986~1992年道路交通噪声平均声级数据见表1，为近年来交通噪声数据【dB(A)】

序号年份	eq L
1986	71.1
1987	72.4
1988	72.4
1989	72.1
1990	71.4
1991	72.0
1992	71.6

第一步：级比检验
\quad建立交通噪声平均声级数据时间序列：x(0)=(x1(0),x2(0),⋯ ,x7(0))=(71.1,72.4,72.4,72.1,71.4,72.0,71.6)x^{(0)}=(x^{(0)}_1,x^{(0)}_2,\cdots,x^{(0)}_7)=(71.1,72.4,72.4,72.1,71.4,72.0,71.6)x(0)=(x1(0),x2(0),⋯,x7(0))=(71.1,72.4,72.4,72.1,71.4,72.0,71.6)

\quad（1）求级比λk\lambda_kλk：λk=xk−1(0)xk(0)\lambda_k=\frac{x^{(0)}_{k-1}}{x^{(0)}_k}λk=xk(0)xk−1(0)
λ=(λ2,λ3,⋯ ,λ7)=(0.982,1,1.0042,1.0098,0.9917,1.0056)\lambda=(\lambda_2,\lambda_3,\cdots,\lambda_7)=(0.982,1,1.0042,1.0098,0.9917,1.0056)λ=(λ2,λ3,⋯,λ7)=(0.982,1,1.0042,1.0098,0.9917,1.0056)

\quad（2）级比判断：λk(0)∈(e−2n+1,e2n+1)\lambda^{(0)}_k\in (e^{-\frac{2}{n+1}},e^{\frac{2}{n+1}})λk(0)∈(e−n+12,en+12)
由于所有的λk∈[0.982,1.0098],k=2,3,⋯ ,7\lambda_k\in[0.982,1.0098],k=2,3,\cdots,7λk∈[0.982,1.0098],k=2,3,⋯,7则可以用x(0)x^{(0)}x(0)进行GM（1，1）建模。

第二步：GM(1,1)建模
\quad（1）对原始数据x(0)x^{(0)}x(0)作一次累加
x(1)=(71.1,143.5,215.9,288,359.4,431.4,503)x^{(1)}=(71.1,143.5,215.9,288,359.4,431.4,503)x(1)=(71.1,143.5,215.9,288,359.4,431.4,503)

\quad（2）构造数据矩阵B及数据向量Y
B=[−12(x1(1)+x2(1))1−12(x2(1)+x3(1))1⋮−12(x6(1)+x7(1))1],Y=[x2(0)x3(0)⋮x7(0)]B=\begin{bmatrix}-\frac{1}{2}(x^{(1)}_1+x^{(1)}_2)\qquad 1\\[1ex]-\frac{1}{2}(x^{(1)}_2+x^{(1)}_3)\qquad 1\\[1ex]\vdots\\[1ex]-\frac{1}{2}(x^{(1)}_6+x^{(1)}_7)\qquad 1\end{bmatrix},\quad Y=\begin{bmatrix}x^{(0)}_2\\[1ex]x^{(0)}_3\\[1ex]\vdots\\[1ex] x^{(0)}_7\end{bmatrix}B=⎣⎢⎢⎢⎢⎢⎢⎡−21(x1(1)+x2(1))1−21(x2(1)+x3(1))1⋮−21(x6(1)+x7(1))1⎦⎥⎥⎥⎥⎥⎥⎤,Y=⎣⎢⎢⎢⎢⎢⎢⎡x2(0)x3(0)⋮x7(0)⎦⎥⎥⎥⎥⎥⎥⎤

\quad（3）计算μ\muμ
μ=[ab]=(BTB)−1BTY=[0.002372.6573]\mu=\begin{bmatrix}a\\b\end{bmatrix}=(B^TB)^{-1}B^TY=\begin{bmatrix}0.0023\\72.6573\end{bmatrix}μ=[ab]=(BTB)−1BTY=[0.002372.6573]

\quad（4）建立模型
dx(1)dt+0.0023x(1)=72.6573\frac{dx^{(1)}}{dt}+0.0023x^{(1)}=72.6573dtdx(1)+0.0023x(1)=72.6573
\quad\qquad求解得
x^k+1(1)=(x1(0)−ba)e−ak+ba=−30929e−0.0023k+31000\hat{x}^{(1)}_{k+1}=(x^{(0)}_1-\frac{b}{a})e^{-ak}+\frac{b}{a}=-30929e^{-0.0023k}+31000x^k+1(1)=(x1(0)−ab)e−ak+ab=−30929e−0.0023k+31000

\quad（5）求生成数列值x^k+1(1)\hat{x}^{(1)}_{k+1}x^k+1(1)及模型还原值x^k+1(0)\hat{x}^{(0)}_{k+1}x^k+1(0),令：k=1，2，⋯，6k=1，2，\cdots ，6k=1，2，⋯，6，由上面的时间响应函数可算得x^(1)\hat{x}^{(1)}x^(1)，其中取x^1(1)=x^1(0)=x1(0)=71.1\hat{x}^{(1)}_1=\hat{x}^{(0)}_1=x^{(0)}_1=71.1x^1(1)=x^1(0)=x1(0)=71.1由x^k(0)=x^k(1)−x^k−1(1)\hat{x}^{(0)}_k=\hat{x}^{(1)}_k-\hat{x}^{(1)}_{k-1}x^k(0)=x^k(1)−x^k−1(1)，取k=2,3,⋯ ,7k=2,3,\cdots,7k=2,3,⋯,7，得
x^(0)=(x^1(0)，x^2(0)，⋯，x^7(0))=（71.1，72.4，72.2，72.1，71.9，71.7，71.6）\hat{x}^{(0)}=(\hat{x}^{(0)}_1，\hat{x}^{(0)}_2，\cdots，\hat{x}^{(0)}_7)=（71.1，72.4，72.2，72.1，71.9，71.7，71.6）x^(0)=(x^1(0)，x^2(0)，⋯，x^7(0))=（71.1，72.4，72.2，72.1，71.9，71.7，71.6）

第三步：模型检验
\quad模型的各种检验指标的计算结果如下表:

年份	原始值	模型值	残差	相对误差	级比偏差
1986	71.1	71.1
1987	72.4	72.4	-0.0057	0.01%	0.0023
1988	72.4	72.2	0.1638	0.23%	0.0203
1989	72.1	72.1	0.0329	0.05%	-0018
1990	71.4	71.9	-0.4984	0.7%	-0.0074
1991	72.0	71.7	0.2699	0.37%	0.0107
1992	71.6	71.6	0.0378	0.05%	-0.0032

\quad经验证，该模型的精度较高，可进行预测和预报。

【参考文献】

徐华锋，方志耕《优化白化方程的GM(1, 1)模型》
灰色预测模型GM(1,1) 与例题分：https://blog.csdn.net/qq547276542/article/details/77865341
【数学建模】灰色预测及Python实现：https://www.jianshu.com/p/a35ba96d852b

【数学建模常用算法】之灰色预测模型GM相关推荐

【Python数学建模常用算法代码（二）之BP神经网络】
Python数学建模常用算法代码(二) BP神经网络模型Python代码 import numpy as np import math import random import string impo ...
数学建模常用算法—灰色预测
今天数模君给大家讲解一下数学建模比赛中常用的一种预测方法:灰色预测法. 目录模型的含义灰色预测的原理实例模型的含义灰色预测模型 ( Gray Forecast Model )是通过少量的.不 ...
数学建模|预测方法：灰色预测模型
简介灰色系统理论是由华中理工大学邓聚龙教授于1982年提出并加以发展的.二十几年来,引起了不少国内外学者的关注,得到了长足的发展.目前,在我国已经成为社会.经济.科学技术在等诸多领域进行预测.决策. ...
二维动态规划降维误差一般为多少_数学建模常用算法模型
数学模型的分类按模型的数学方法分: 几何模型.图论模型.微分方程模型.概率模型.最优控制模型.规划论模型.马氏链模型等按模型的特征分: 静态模型和动态模型,确定性模型和随机模型,离散模型和连续性模 ...
数学建模常用算法汇总及python，MATLAB实现（五) —— 拟合
拟合比较重要的就是2.1和2.3 2.2可以浅看一下, 自己敲着试一试就拟合部分来说, MATLAB比python强大很多, 自带cftool工具包, 并且有很多快速的函数, 个人建议使用MATL ...
数学建模常用算法—马尔可夫预测
今天数模君带大家学习一下数学建模中的预测算法之马尔科夫预测. 目录模型的含义实例分析模型的含义马尔可夫(Markov)预测法,就是一种关于事件发生的概率预测方法.它是根据事件的目前状况来预测其 ...
数学建模常用算法汇总及python，MATLAB实现（六) —— pandas和matlab实现插值
插值 2的拉格朗日插值用的其实比较少, 可以看一下了解一下插值的原理主要看看3里的结论和4的实现代码文章目录插值 1. 定义 2.拉格朗日插值 2.1 概念 3. Rouge现象 3.1 是什么 ...
数学建模常用算法—灰色关联分析法(GRA)
解决问题灰色关联分析的基本思想是根据序列曲线几何形状的相似程度来判断其联系是否紧密.曲线越接近,相应序列之间的关联度就越大,反之就越小. 一般的抽象系统,如社会系统.经济系统.农业系统.生态系统.教 ...
机器学习数据挖掘十大经典算法数学建模常用算法
国际权威的学术组织the IEEE International Conference on Data Mining (ICDM) 2006 (香港召开)年12月评选出了数据挖掘领域的十大经典算法.不仅 ...

【数学建模常用算法】之灰色预测模型GM