常用的预测方法（如回归分析），需要较大的样本，如果样本较小，常造成较大的误差，使预测目标失效。灰度预测模型（Gray Forecast Model）是通过利用少量的、不完全的信息，建立数学模型并做出预测的一种预测方法。本文将对这种预测模型进行详细介绍。由于部分公式书写麻烦，本文直接插入了公式图片。

1. 灰度系统理论

灰度系统理论是由华中理工大学邓聚龙教授于1982年提出并发展而成的解决灰色系统的理论和方法，该方法将一般的系统论、信息论和控制论的观点与数学方法结合在一起。首先看一下灰度理论中的几个概念。

1.1 灰度系统

灰度系统是指“部分信息已知，部分未知”的“小样本”，“贫信息”的不确定性系统。系统中各因素间有不确定的关系。
作为两个极端，如果系统中信息完全确定，则成为白色系统；相反，信息完全不确定的系统为黑色系统，只能通过它与外界的联系来加以观测研究。区别白色和黑色系统的重要标志是系统中各因素之间是否具有确定关系。

1.2 灰度系统的特点

灰色系统可以充分利用已知信息寻求系统的运动规律，使用灰色数学处理不确定量，使之量化，因此，灰色系统理论能够处理贫信息或小样本的系统。

2. 灰色预测方法

灰度预测法是对灰色系统进行预测的方法。通过关联分析和数据处理来寻找系统变动的规律，生成有较强规律性的数据序列。然后用来预测两种不同的目标，一种是预测某时刻的特征量来预测未来发展趋势。另一种是固定特征量阈值，通过分析特征量出现时刻的规律，预测未来达到阈值的时间点。

2.1 灰色预测类型

(1) 数列预测。用观察到的反映预测对象特征的时间序列构造灰色预测模型，预测未来某一时刻的特征量，或达到某一特征量的时间。
(2) 灾变或异常值预测。利用灰度模型预测异常值出现的特定时间。
(3) 季节灾变或异常值预测。预测灾变或异常值发生在一年内哪个时间区间或季节内。
(4) 拓扑预测。将原始数据做曲线，在曲线上按照阈值寻找发生的所有时间点，并以该时间点构成的序列构建预测模型，建立模型预测该阈值的时间点。
(5) 系统预测。 通过对系统行为特征指标建立一组相关联的灰度预测模型，预测系统中众多变量间的相互协调关系的变化。

下面我们对具体的集中预测模型进行介绍。

3. GM(1,1)模型

3.1 数据/数列生成

我们获取到的未作处理的数据被称为原始数列，按照某种要求经过处理后生成的数列为生成数列，我们可以从中发现原始数列中没有表现出的内在规律。
常用的数据处理方法有累加生成、累减生成、均值生成。
(1) 累加生成
通过依次累加某时刻之前的数据得到新的数列，记为AGO(Accumulating Ceneration Operator)。
设x^(0)为原始序列，既有
![][1]
[1]: http://latex.codecogs.com/svg.latex?$x^{(0)}=[x{(0)}(1),x^{{(0)}(2),...,x}{(0)}(n)]$
按照下列公式生成新的累加数列：
![][2]
[2]: http://latex.codecogs.com/svg.latex?$x^{{(1)}(k)=\sum_{i=1}}{k}x^{(0)}(i);k=1,2,...,n$
则称x^(1)为x(0)的一次累加生成数列：
![][3]
[3]: http://latex.codecogs.com/svg.latex?$x^{(1)}=[x{(1)}(1),x^{{(1)}(2),...,x}{(1)}(n)]$
类似的有r次累加生成：
![][4]
[4]: http://latex.codecogs.com/svg.latex?$x^{{(r)}(k)=\sum_{i=1}}{k}x^{(r-1)}(i);k=1,2,...,n$
(2)累减生成
对原始数列的数据依次作前后两个数据相减的运算，记为IAGO(Inverse Accumulating Ceneration Operator)。累减生成可以将累加生成序列还原成非生成序列。
我们设x^(1)为原始序列，则有
![][5]
[5]: http://latex.codecogs.com/svg.latex?$x^{(1)}=[x{(1)}(1),x^{{(1)}(2),...,x}{(1)}(n)]$
如果令：
![][6]
[6]: http://latex.codecogs.com/svg.latex?$x^{(0)}(k)=x{(1)}(k)-x^{(1)}(k-1)$
则可生成一次累减序列x^(0)：
![][7]
[7]: http://latex.codecogs.com/svg.latex?$x^{(0)}=[x{(0)}(1),x^{{(0)}(2),...,x}{(0)}(n)]$
更一般地，有：
![][8]
[8]: http://latex.codecogs.com/svg.latex?$x^{(r-1)}(k)=x{(r)}(k)-x^{(r)}(k-1)$
均值生成序列分为邻均值生成和非邻均值生成两种，就暂时不在本文详细介绍了，以后再补充。

3.2模型建立

给定观测数据序列
![][9]
[9]: http://latex.codecogs.com/svg.latex?$x^{(0)}=[x{(0)}(1),x^{{(0)}(2),...,x}{(0)}(n)]$
经过一次累加后，得到序列：
![][10]
[10]: http://latex.codecogs.com/svg.latex?$x^{(1)}=[x{(1)}(1),x^{{(1)}(2),...,x}{(1)}(n)]$
设x^(1)满足一阶常微分方程
![][11]
[11]: http://latex.codecogs.com/svg.latex?${dx^{{(1)}}/{dt}+ax}{(1)}=u$
其中a是发展灰数或发展系数，u为内生控制灰数或灰作用量，是对系统的常定输入。该方程满足当初始条件为：
![][12]
[12]: http://latex.codecogs.com/svg.latex?$when:t=t_{0},x^{(1)}=x{(1)}(t_0)$
则有解为：
![][13]
[13]: http://latex.codecogs.com/svg.latex?$x^{(1)}(t)=[x{(1)}(t_0)-u/a]e^{-a(t-t_0)}+u/a$
对等间距取样的离散值,则有
![][19]
[19]: http://latex.codecogs.com/svg.latex?$when:t_0=1\\x^{{(1)}(k+1)=[x}{(1)}(1)-u/a]e^{-ak}+u/a$
由于x^{(1)对t求导涉及x}(1)的两个时刻值，因此，x^(1)(i)取前后两个时刻的平均代替更合理。即有：
![][14]
[14]: http://latex.codecogs.com/svg.latex?$1/2[x^{(r)}(i)+x{(r)}(i-1)],(i=2,3,...,N)$
写成矩阵表达的形式为：
![][15]
[15]: http://latex.codecogs.com/svg.latex?$\begin{bmatrix}x^{(0)}(2)\x{(0)}(3)\...\x^{{(0)}(N)\end{bmatrix}=\begin{bmatrix}-1/2[x}{(1)}(2)+x^{{(1)}(1)]&1\-1/2[x}{(1)}(3)+x^{{(1)}(2)]&1\...\-1/2[x}{(1)}(N)+x^{(1)}(N-1)]&1\end{bmatrix}\begin{bmatrix}a\u\end{bmatrix}$
如果用符号代表上式的向量或矩阵，则有
![][16]
[16]: http://latex.codecogs.com/svg.latex?$y=BU$
其中，
![][17]
[17]: http://latex.codecogs.com/svg.latex?$U=[a,u]^{T}\\y=[x{(0)}(2),x^{{(0)}(3),...,x}{(0)}(N)]^{{T}\\B=\begin{bmatrix}-1/2[x}{(1)}(2)+x^{{(1)}(1)]&1\-1/2[x}{(1)}(3)+x^{{(1)}(2)]&1\...\-1/2[x}{(1)}(N)+x^{(1)}(N-1)]&1\end{bmatrix}$
利用最小二乘法估计上式，可得：
![][18]
[18]: http://latex.codecogs.com/svg.latex?\hat{U}=\begin{bmatrix}\hat{a}\\hat{u}\end{bmatrix}=(B^{T}B){-1}B^{T}y
将上式得到的估计值a和u带入求离散解的公式，得到时间响应方程
![][20]
[20]: http://latex.codecogs.com/svg.latex?\hat{x^{{(1)}}(k+1)=[x}{(1)}(1)-\frac{\hat{u}}{\hat{a}}]e^{-\hat{a}k}+\frac{\hat{u}}{\hat{a}}
当k=1,2,...,N-1时，由上式计算得到的结果是拟合值，当k>=N时，结果为预测值。这是相对于一次累加序列的拟合值，用后减运算还原，就可得原始序列的拟合值和预测值。

3.3 精度检验

(1) 残差检验：预测得到的一次累加序列进行累减计算，得到原始序列估计值，再计算其与原始序列的绝对残差序列和相对残差序列。
残差：
![][21]
[21]: http://latex.codecogs.com/svg.latex?E(k)=x^{{(0)}(k)-\hat{x}{(0)}}(k),k=2,3,...N;
相对残差：
![][22]
[22]: http://latex.codecogs.com/svg.latex?e(k)=[x^{{(0)}(k)-\hat{x}{(0)}}(k)]/x^{(0)}(k),k=2,3,...N;
(2) 后验差检验：分别计算
原始序列x^(0)的均值：
![][23]
[23]: http://latex.codecogs.com/svg.latex?\bar{X}=\frac{1}{N}\sum_{k=1}^{N}x{(0)}(k)
原始序列x^(0)的方差：
![][24]
[24]: http://latex.codecogs.com/svg.latex?S_1=\sqrt{\frac{1}{N}\sum_{k=1}^{N}[x{(0)}(k)-\bar{X}]^2}
残差的均值：
![][25]
[25]: http://latex.codecogs.com/svg.latex?\bar{E}=\frac{1}{N-1}\sum_{k=2}^{N}E(k)
残差的方差：
![][26]
[26]: http://latex.codecogs.com/svg.latex?S_2=\sqrt{\frac{1}{N-1}\sum_{k=2}^{{N}[E(k)-\bar{E}]}2}
后验差比值：
![][27]
[27]: http://latex.codecogs.com/svg.latex?C=\frac{S_2}{S_1}
小误差概率：
![][28]
[28]: http://latex.codecogs.com/svg.latex?P=P{|E(k)-\bar{E}|<0.6745S_1}
对比精度等级对照表，如下表所示。
![][29]
[29]: http://latex.codecogs.com/svg.latex?\begin{array}{c|lcr}level&\text{P}&\text{C}&\text{Level}\\hline1&{\textgreater0.95}&{\textless0.35}&Good\2&{\textgreater0.80}&{\textless0.45}&Qualified\3&{\textgreater0.70}&{\textless0.5}&Barely-passing\4&\leq0.70&\geq0.65&Unqualified\end{array}
评价模型的指标C越小越好，P越大越好，一般按照Max{P的level,C的level}来确定模型的精度级别。
(3) 关联度检验
关联度分析时根据因素之间发展态势的相似或相异程度，分析事物动态关联的特征和程度。后期再作补充。

由于灰度预测模型是给予一阶常微分方程建立的，故被称为一阶一元灰色模型，记为GM(1,1)。特别注意，原始数据中如果有负数，需要先进行“数据整体提升”处理，再做一次累加，避免累加时正负抵消。

4. GM(1,1)残差模型

如果用原始时间序列建立的GM(1,1)模型检验不合格或精度不理想时，要对建立的模型进行残差修正，可以提高模型的预测精度。
我们设
![][30]
[30]: http://latex.codecogs.com/svg.latex?\varepsilon^{{(0)}=(\varepsilon}{(0)}(1),\varepsilon^{{(0)}(2),...,\varepsilon}{(0)}(n))
其中，我们称
![][31]
[31]: http://latex.codecogs.com/svg.latex?\varepsilon(k)=x^{{(0)}(k)-\hat{x}}{(1)}(k)
为X^(1)的残差序列。若存在k0满足如下条件
![][32]
[32]: http://latex.codecogs.com/svg.latex?1.\forall{k}\gek_0,\varepsilon^{(0)}(k)\textless0or*\varepsilon^{(0)}\textgreater0;\\****2.when:n-k_0\ge4;
则称
![][33]
[33]: http://latex.codecogs.com/svg.latex?(|\varepsilon^{{(0)}(k_0)|,|\varepsilon}{(0)}(k_0+1)|,...,|\varepsilon^{(0)}(n)|)
为可建模残差尾段，仍记为
![][34]
[34]: http://latex.codecogs.com/svg.latex?(\varepsilon^{{(0)}(k_0),\varepsilon}{(0)}(k_0+1),...,\varepsilon^{(0)}(n))
对其进行一次累加，得到序列
![][35]
[35]: http://latex.codecogs.com/svg.latex?(\varepsilon^{{(1)}(k_0),\varepsilon}{(1)}(k_0+1),...,\varepsilon^{(1)}(n))
利用GM(1,1)建模后的时间响应式为：
![][36]
[36]: http://latex.codecogs.com/svg.latex?\varepsilon^{{(1)}(k+1)=(\varepsilon}{(0)}(k_0)-\frac{\hat{u}}{\hat{a}})e^{[-a(k-k_0)]}+\frac{\hat{u}}{\hat{a}},{k}\ge{k_{0}}
则残差尾段的模拟序列为：
![][37]
[37]: http://latex.codecogs.com/svg.latex?\hat{\varepsilon}^{{(0)}=(\hat{\varepsilon}}{(0)}(k_0),\hat{\varepsilon}^{{(0)}(k_0+1),...,\hat{\varepsilon}}{(0)}(n))
其中
![][39]
[39]: http://latex.codecogs.com/svg.latex?\hat{\varepsilon}^{{(0)}(k+1)=-\hat{a}(\hat{\varepsilon}}{(0)}(k_0)-\frac{\hat{u}}{\hat{a}})e^{[-\hat{a}(k-k_0)]},k\ge{k_0}
如果用上式修正X^(1)，则称修正后的时间响应式为
![][38]
[38]: http://latex.codecogs.com/svg.latex?\hat{\varepsilon}^{{(1)}(k+1)=\left{\begin{array}{ll}(x}{(0)}(1)-\frac{u}{a})e^{{-ak}+\frac{u}{a},k\textless{k_0};(1)\\(x}{(0)}(1)-\frac{u}{a})e^{{-ak}+\frac{u}{a}\pm\hat{a}(\varepsilon}{(0)}(k_0)-\frac{\hat{u}}{\hat{a}})e^{[-\hat{a}(k-k_0)]},k\ge{k_0};(2)\end{array}\right.
因此为残差修正GM(1,1)模型，简称残差GM(1,1)

5. GM(n,h)模型

GM(n,h)模型是微分方程模型，可用于对描述对象做长期、连续、动态的反映。从原则上讲，某一灰色系统无论内部机制如何，只要能将系统原始数据表示为原始时间序列，并保证值非负，则可用GM模型对系统表示和预测。