美赛数学建模之灰色预测

一、灰色预测概述
1.原理：通过鉴别系统因素之间发展趋势的相异程度，并对原始数据进行生成处理来寻找系统变动的规律，生成有较强规律性的数据序列，然后建立相应的微分方程模型，从而预测事物的未来发展趋势。灰色预测的数据是通过生成数据的模型所得到的预测值的逆处理结果。

2.特点：在求解过程中需要进行数据预测时，利用灰色模型预测的结果比较稳定，不仅适用于大数据量的预测，在数据量较少时（数据多于3个即可）预测结果依然准确。

二、单序列一阶线性微分方程模型GM(1,1)介绍
设有原始数据列 $x^{(0)}=(x^{(0)}(1),x^{(0)}(2),\cdots ,x^{(0)}(n))$ ，n为数据个数。
如果根据 $x^{(0)}$ 数据列建立GM(1,1)来实现预测功能，则基本步骤如下：
1.原始数据累加以便弱化随机序列的波动性和随机性，得到新数据序列：$$x^{(1)}=(x^{(1)}(1),x^{(1)}(2),\cdots ,x^{(1)}(n))$$ 其中，$x^{(1)}(t)$ 中各数据表示对应前几项数据的累加。 $$x^{(1)}(t)=\sum _{k=1}^t{x}^{(0)}(k)，t=1,2,\cdots ,n$$
2.对 $x^{(1)}(t)$ 建立 $x^{(1)}(t)$ 的一阶线性微分方程：$$\frac{d{x}^{(1)}}{dt}+a{x}^{(1)}=u$$ 其中，a , u 为待定系数，分别称为发展系数和灰色作用量，a 的有效区间是 $(-2,2)$,并记 a , u 构成的矩阵 $\hat{a}=\left(\begin{array}{c}a\\ u\end{array}\right)$ 。只要求出参数 a , u ，就能求出 $x^{(1)}(t)$ ，进而求出 $x^{(0)}(t)$ 的未来预测值。

3.对累加生成数据做均值生成 B 与常数项向量 Y ，即
$B=\left[\begin{array}{c}0.5(x^{(1)}(1)+x^{(1)}(2))\\ 0.5(x^{(1)}(2)+x^{(1)}(3))\\ \cdots \\ 0.5(x^{(1)}(n-1)+x^{(1)}(n))\end{array}\right]$ ， $Y_n=(x^{(0)}(2),x^{(0)}(3),\cdots ,x^{(0)}(n){)}^T$

4.用最小二乘法求解灰参数 $\hat{a}$ , 则$$\hat{a}=\left(\begin{array}{c}a\\ u\end{array}\right)=(B^{T}B{)}^{-1}{B}^T{Y}_n$$
5.将灰参数 $\hat{a}$ 代入 $\frac{d{x}^{(1)}}{dt}+a{x}^{(1)}=u$ ，并对其进行求解，得$${\hat{x}}^{(1)}(t+1)=(x^{(0)}(1)-\frac{u}{a})e^{-at}+\frac{u}{a}$$ 由于 $\hat{a}$ 是通过最小二乘法求出的近似值，所以 ${\hat{x}}^{(1)}(t+1)$ 是一个近似表达式，为了与原序列 $x^{(1)}(t+1)$ 区分开来，故记为 ${\hat{x}}^{(1)}(t+1)$ 。

6.对函数表达式 ${\hat{x}}^{(1)}(t+1)$ 及 ${\hat{x}}^{(1)}(t)$ 进行离散，并将二者做差以便还原 $x^{(0)}$ 原序列，得到近似数据序列 ${\hat{x}}^{(0)}(t+1)$ 如下：$${\hat{x}}^{(0)}(t+1)={\hat{x}}^{(1)}(t+1)-{\hat{x}}^{(1)}(t)$$
7.对建立的灰色模型进行检验，步骤如下：
(1)计算 $x^{(0)}$ 与 ${\hat{x}}^{(0)}(t)$ 之间的残差 $e^{(0)}(t)$ 和相对误差 $q(x)$ ：
$$e^{(0)}(t)=x^{(0)}-{\hat{x}}^{(0)}(t)$$ $$q(x)=\frac{e^{(0)}(t)}{x^{(0)}(t)}$$
(2)求原始数据 $x^{(0)}$ 的均值以及方差 $s_1$ 。
(3)求 $e^{(0)}(t)$ 的平均值 $\overline{q}$ 以及残差的方差 $s_2$ 。
(4)计算方差比 $C=\frac{s_1}{s_2}$ 。
(5)求小误差概率 P = P { ∣ e ( t ) ∣ < 0.6745 s 1 } P = P\{ |e(t)|<0.6745s_1 \} P=P{∣e(t)∣<0.6745s1​} 。
(6)灰色模型精度检验如下表所列。

等级	相对误差q	方差比C	小概率误差P
I级	<0.01	<0.35	>0.95
II级	<0.05	<0.50	<0.80
III级	<0.10	<0.65	<0.70
IV级	>0.20	>0.80	<0.60

8.利用模型进行预测：$${\hat{x}}^{(0)}=\left[\begin{array}{c}{\hat{x}}^{(0)}(1),{\hat{x}}^{(0)}(2),\cdots ,{\hat{x}}^{(0)}(n),{\hat{x}}^{(0)}(n+1),\cdots ,{\hat{x}}^{(0)}(n+m)\end{array}\right]$$
其中，${\hat{x}}^{(0)}(1),{\hat{x}}^{(0)}(2),\cdots ,{\hat{x}}^{(0)}(n)$ 为原数列的模拟；
而 ${\hat{x}}^{(0)}(n+1),\cdots ,{\hat{x}}^{(0)}(n+m)$ 为未来数列的预测。

三、基于MATLAB R2014a的灰色预测实现方法
1.典型程序结构(可完全按照预测模型的求解步骤)
(1)对原始数据进行累加。
(2)构造累加矩阵 B 与常数向量。
(3)求解灰参数。
(4)将参数带入预测模型进行数据预测

2.实例
已知某公司1999-2008年的利润为(单位：元/年)：[89677，99215，109655，120333，135823，159878，182321，209407，246619，300670]
具体的MATLAB程序如下：

clc
syms a b;                           % 创建符号变量(用于符号运算的变量)a ,b
c = [a b]';                         % 创建一个2 * 1的符号变量矩阵c
start = 1999;                       % 第一个数据的横坐标(起始年)
finish = 2018;                      % 想要预测到的数据的横坐标(目标年)
% % 数据读入
A = [89677,99215,109655,120333,135823,159878,182321,209407,246619,300670];
B = cumsum(A);                      % 原始数据累加
n = length(A);                      % 原始数据个数
% % 生成累加矩阵
for i = 1 : (n - 1)
    C(i) = (B(i) + B(i + 1)) / 2;
end
% % 生成常数向量
D = A; D(1) = [];                   
D = D';
% % 最小二乘法计算灰参数
E = [-C; ones(1,n-1)];
c = (E * E') \ E * D;
c = c';
a = c(1); b = c(2);                 % 将算得得灰参数传入变量a , b中
% % 预测后续数据
F = []; F(1) = A(1);
for i = 2 : (finish - start + 1)
    F(i) = (A(1) - b / a) / exp(a * (i - 1)) + b / a;
end
G = []; G(1) = A(1);
for i = 2 : (finish - start + 1)
    G(i) = F(i) - F(i - 1);         % 得到预测出来的数据
end
% % 设定时间序列
t1 = start : (start + n - 1);       % 原始数据的时间序列
t2 = start : finish;                % 预测数据的时间序列
% % 预测图输出
plot(t1,A,'o',t2,G);                % 原始数据与预测数据的比较

得到的预测图输出如下：

程序说明：

• 接口已经在程序开头拉出，因此在实际使用时可以直接套用该程序，只需要把原数据和时间序列变量 start 和 finish 替换就可以了。
• 模型误差检验可以灵活处理，这里给出的是预测数据与原始数据的比较图，实际运用时也可以对预测数据进行其他方式的精度校验。