灰色预测及其MATLAB实现——MATLAB在数学建模中的应用

灰色预测及其MATLAB实现

灰色预测是一种常规的预测手段，具有操作简便，所需数据量少等优点，一般只需有4个数据就可以进行预测。

灰色预测是基于灰色系统理论的预测方法。灰色系统由我国著名学者邓聚龙教授在1982年提出，是相对于“白色模型”——完全信息透明的模型，和“黑色模型”——对信息一无所知的模型的模型概念。利用灰色系统解决的问题主要是具有不确定性的问题：

1. 信息具有模糊性，无法用数学方程精确刻画。
2. 机理具有不确定性。
3. 信息贫瘠的不确定性。

灰色系统基本理论

灰色关联度矩阵

灰色关联度分析了向量与向量之间以及矩阵与矩阵之间的关联度，而向量与向量之间的关联度可以看作矩阵与矩阵之间的关联度的一种特殊形式。

计算关联度，一定是计算某个待比较的数列与参照物（亦即参考数列）之间的相关程度。

假设有一组参照数列如下：
$$x_j=(x_j(1),x_j(2),x_j(3),\cdots ,x_j(k),\cdots ,x_j(n))j=1,2,3,\cdots ,s$$
假设有一组待比较数组如下：
$$x_i=(x_i(1),x_i(2),x_i(3),\cdots ,x_i(k),\cdots ,x_i(n))i=1,2,3,\cdots ,t$$
则定义关联系数如下：
$${\xi }_{ji}(k)=\frac{\underset{i}{\min }\underset{k}{\min }|x_j(k)-x_i(k)|+\rho \cdot \underset{i}{\max }\underset{k}{\max }|x_j(k)-x_i(k)|}{|x_j(k)-x_i(k)|+\rho \cdot \underset{i}{\max }\underset{k}{\max }|x_j(k)-x_i(k)|}$$
关于该式的说明如下：

1. 变量${\xi }_{ji}(k)$表示的是第$i$个数列与第$j$个参考数列第$k$个样本之间的关联系数。
2. $\underset{i}{\min }\underset{k}{\min }|x_j(k)-x_i(k)|$ 和$\underset{i}{\max }\underset{k}{\max }|x_j(k)-x_i(k)|$表示参考数列矩阵与比较数列矩阵作差后的最小值和最大值，目的是为了保证${\xi }_{ji}(k)$的值在[0,1]区间内，同时上下对称的结构可以消除量纲不同和数量级悬殊的问题。
3. $|x_j(k)-x_i(k)|$即为汉明距离（“Hamming distance”），汉明距离的倒数被称为反倒数距离，灰色关联度的本质就是通过反倒数距离的大小来判定关联程度，倒数越大，表示两条曲线之间距离越小，其曲线程度越相似。
4. $\rho$的取值约定俗成为[0,1]，但实际上$\rho$的取值范围为$(0,+\infty )$。但不管$\rho$如何取值，其只改变${\xi }_{ji}(k)$的绝对大小，而不改变关联度的相对强弱。

由于${\xi }_{ji}(k)$只能反映处点与点之间的相关性，相关性信息分散，不方便刻画数列之间的相关性，需要把${\xi }_{ji}(k)$整合起来，定义：
$$r_{ji}=\frac{\stackrel{n}{\underset{k=1}{\Sigma }}{\xi }_{ji}(k)}{n}$$
变量$r_{ji}$即为相关度，结合实际背景，有正面作用的即为正相关，反之则为负相关；$r_{ji}$大于0.7称为强相关，小于0.3称为弱相关。

将$x_i$与$x_j$之间的相关度写为矩阵形式：
$$R=\left[\begin{array}{cccc}{r}_{11}& r_{12}& \cdots & r_{1t}\\ r_{21}& r_{22}& \cdots & r_{2t}\\ ⋮& ⋮& \ddots & ⋮\\ r_{s1}& r_{s2}& \cdots & r_{st}\end{array}\right]$$
则通过观察某一列数值明显大于其他列数值，则称该列为优势子因素；假如某一行数值明显大于其他行数值，则称该行为优势母因素，优势母因素易受子因素的驱动影响。

经典灰色模型GM(1,1)

灰色模型是基于常微分方程的基础上发展起来的。灰色模型记号是$GM(M,N)$，前面的$N$表示变量个数，后面的$M$表示常微分方程的阶数，一般情况下$GM(1,1)$最为常见。

下面是$GM(1,1)$模型组建过程：

1. 数据累加：降低数据噪声，加强数学规律的显露。

   设有原始数列:$x^{(0)}=(x^{(0)}(1),x^{(0)}(2),\cdots ,x^{(0)}(n))$

   现对$x^{(0)}$进行累加得到新数列$x^{(1)}$:
   $$\begin{array}{rl}{x}^{(1)}& =(x^{(1)}(1),x^{(1)}(2),\cdots ,x^{(1)}(n))\\ & =(x^{(0)}(1),x^{(0)}(1)+x^{(0)}(2),x^{(0)}(1)+x^{(0)}(2)+x^{(0)}(3)\cdots )\end{array}$$

2. 计算数列$x^{(1)}$的灰导数方程

   假设一个连续函数$C=f(t)$，取横坐标$t_{k-2},t_{k-1},t_k,t_{k+1},\cdots$,各相邻点间距为$h$，$h$可以是相等的，也可以不相等。当$h\to 0$时，有:
   $$(\frac{dC}{dt}{)}_k=\frac{C_k-C_{k-1}}{h}$$
   由于$h$是相邻数值之间的距离也就是差值，则上式变为：
   $$(\frac{dC}{dt}{)}_k=\frac{C_k-C_{k-1}}{h}=\frac{x^{(1)}(k)-x^{(1)(k-1)}}{k-(k-1)}=x^{(1)}(k)-x^{(1)}(k-1)=x^{(0)}(k)=d(k)$$
   这一步需要$x^{(1)}$数列的值变化比较平滑，如果出现突变的数据序列，需要将突变提取出来后建立灰色模型。

3. 定义数列$x^{(1)}$的紧邻均值
   $$z^{(1)}(k)=\frac{x^{(1)(k)+x^{(1)}(k-1)}}{2}$$
   则称$z^{(1)}(k)=(z^{(1)}(2),z^{(1)}(3),z^{(1)}(4),\cdots ,z^{(1)}(n))$为$x^{(1)}$的紧邻均值数列，紧邻均值数列用于防止原始数据自身存在突变的奇异数据，平滑数列的阶跃性。

4. 定义GM(1,1)灰微分方程

   定义GM(1,1)灰微分方程为一阶线性方程
   $$d(k)+a{z}^{(1)}(k)=b$$
   可以得到:
   $$x^{(0)}(k)+a{z}^{(1)}(k)=b$$
   该式中，$a,b$是待辨识的参数，$a$被称为“发展系数”，$b$被称为“灰作用量”。将$k=2,3,\cdots ,n$代入，用最小二乘法可以得到$a,b$的值。

5. 白化GM(1,1)模型

   原来的方程$x^{(0)}(k)+a{z}^{(1)}(k)=b$是一个离散变量的方程，将其白化为对应的连续方程，这样就可以化解成我们最终需要的灰色模型形式。
   对$x^{(0)}(k)+a{z}^{(1)}(k)=b$白化后可以得到：
   $$\frac{d{x}^{(1)}}{dk}+a{x}^{(1)}=b$$
   对该微分方程求解易得：
   $$x^{(1)}(k+1)=\frac{b}{a}+(x^{(0)}(1)-\frac{b}{a})e^{-ak}$$
   最后累减还原获得最终序列

6. 模型检验

   这里介绍了三种常用的检验模型精度的方式:

   1. 相对残差$Q$检验

      原始数列为 $x^{(0)}=(x^{(0)}(1),x^{(0)}(2),\cdots ,x^{(0)}(n))$

      对应灰色模型预测得出的序列:
      $$x^{(0{)}^{\prime }}=(x^{(0{)}^{\prime }}(1),x^{(0{)}^{\prime }}(2),\cdots ,x^{(0{)}^{\prime }}(n))$$
      残差序列为：
      ϵ(0)=(ϵ1,ϵ2,ϵ3,⋯ ,ϵn)={x(0)−x(0)′∣x(0)(1)−x(0)′(1),x(0)(2)−x(0)′(2),⋯ ,x(0)(n)−x(0)′(n)} \begin{aligned} \epsilon^{(0)}&=(\epsilon_1,\epsilon_2,\epsilon_3,\cdots,\epsilon_n)\\ &=\{x^{(0)}-x^{(0)'}|x^{(0)}(1)-x^{(0)'}(1),x^{(0)}(2)-x^{(0)'}(2),\cdots,x^{(0)}(n)-x^{(0)'}(n)\} \end{aligned} ϵ(0)​=(ϵ1​,ϵ2​,ϵ3​,⋯,ϵn​)={x(0)−x(0)′∣x(0)(1)−x(0)′(1),x(0)(2)−x(0)′(2),⋯,x(0)(n)−x(0)′(n)}​
      相对误差序列为：
      $$\begin{array}{rl}\Delta & =({\Delta }_1,{\Delta }_2,{\Delta }_3,\cdots ,{\Delta }_n)\\ & =(\frac{{ϵ}_1}{x^{(0)}(1)},\frac{{ϵ}_2}{x^{(0)}(2)},\frac{{ϵ}_3}{x^{(0)}(3)},\cdots ,\frac{{ϵ}_n}{x^{(0)}(n)})\end{array}$$
      所以相对误差$Q=\frac{1}{n}\stackrel{n}{\underset{k=1}{\Sigma }{\Delta }_k}$,$Q$越小越好，表明模型越精确，根据$Q$的大小来划分模型的精度等级。

   2. 方差比$C$的检验

      $x^{(0)}$是原始数列，$x^{(0{)}^{\prime }}$是模型预测数列，${ϵ}^{(0)}$是残差序列。
      则各自的序列和方差为:
      $$\{\begin{array}{rl}\overline{x}& =\frac{1}{n}\stackrel{n}{\underset{k=1}{\Sigma }}{x}^{(0)}(k)\\ S_1^2& =\frac{1}{n}\stackrel{n}{\underset{k=1}{\Sigma }}(x^{(0)}-\overline{x}{)}^2\\ \overline{ϵ}& =\frac{1}{n}\stackrel{n}{\underset{k=1}{\Sigma }}{ϵ}^{(0)}\\ S_2^2& =\frac{1}{n}\stackrel{n}{\underset{k=1}{\Sigma }}({ϵ}^{(0)}-\overline{ϵ})\end{array}$$
      取$C=\frac{S_2}{S_1}$,$C$越小，可信度越高。

   3. 小概率误差$P$的检验

      p=P{∣ϵ−x‾∣0.6745S1} p=P\{|\epsilon-\overline{x}|0.6745S_1\} p=P{∣ϵ−x∣0.6745S1​}
      这里的置信水平取0.5，所以置信区间半长取0.6745，其意义是误差摆动$|{ϵ}_k-\overline{ϵ}|$位于指定数值范围内的概率，此指标可以反映误差摆动幅度落在指定区间的数量。因此$p$越大，表示吻合精度越好。

   灰色模型精度检验对照表

   等级	相对误差$Q$	方差比$C$	小误差概率$p$
   一级	$<0.01$	$<0.35$	$>0.95$
   二级	$<0.05$	$<0.50$	$<0.80$
   三级	$<0.10$	$<0.65$	$<0.70$
   四级	$>0.20$	$>0.80$	$<0.60$

灰色系统的程序设计

灰色关联度矩阵的程序设计

	1979	1980	1981	1982	1983
固定资产投资	308.58	310	295	346	367
工业投资	195.4	189.9	187.2	205	222.7
农业投资	24.6	21	12.2	15.1	14.57
科技投资	20	25.6	23.3	29.2	30
交通投资	18.98	19	22.3	23.5	27.655
国民收入	170	174	197	216.4	235.8
工业收入	57.55	70.74	76.8	80.7	89.85
农业收入	88.56	70	85.38	99.83	103.4
商业收入	11.19	13.28	16.82	18.9	22.8
交通收入	4.03	4.26	4.34	5.06	5.78
建筑业收入	13.7	15.6	13.77	11.98	13.95

function [R]=GrayConnect(X,Y)
[xa,xb]=size(X);
[ya,yb]=size(Y);
if (xb==yb)
else
    return ;
end
R=zeros(ya,xa);
q=0.5;
for i = 1:ya
    k=zeros(xa,xb);
    for j=1:xa
        k(j,:)=abs(X(j,:)-Y(i,:));
    end
    temp1=min(min(k));
    temp2=q*max(max(k));
    for j=1:xa
        sum=0;
        for t=1:xb
        sum=sum+(temp1+temp2)/(abs(X(j,t)-Y(i,t))+temp2);
        end
        R(i,j)=sum/xb;
    end
end
end

将数据导入MATLAB,

A=xlsread('Grey.xlsx','B1:F11');

将数据标准化,

for i=1:11
    A(i,:)=A(i,:)/A(i,1);
end

将母因素和子因素分割为X和Y,并使用函数GrayConnect,

X=A(1:5,:);
Y=A(6:11,:);
[R]=GrayConnect(X,Y)

得到结果如下：
$$R=\left[\begin{array}{ccccc}0.8016& 0.7611& 0.5570& 0.8102& 0.9355\\ 0.6887& 0.6658& 0.5287& 0.8854& 0.8004\\ 0.8910& 0.8581& 0.5786& 0.5773& 0.6749\\ 0.6776& 0.6634& 0.5675& 0.7800& 0.7307\\ 0.8113& 0.7742& 0.5648& 0.8038& 0.9205\\ 0.7432& 0.7663& 0.5616& 0.6065& 0.6319\end{array}\right]$$

GM(1,1)的程序设计

银行有各种理财产品，客户可以根据自己的资金实力和投资偏好来选择。有一款理财产品，有十天的“10天犹豫期”，在这10天里如果对自己购得的理财产品不放心或者不满意，可以自由退买，不需要支付手续费。有如下八天顾客保有量：
$$\left[\begin{array}{ccccccc}92.810& 97.660& 98.800& 99.281& 99.537& 99.817& 100.000\end{array}\right]$$

function [XR]=GM1_1(x0)
[a,b]=size(x0);
x1=zeros(a,b);
z1=zeros(a,b-1);
for i =1:b %建立数列x1
    for j = 1:i
        x1(1,i)=x1(1,i)+x0(1,j);
    end
end
for i=2:b
    z1(1,i-1)=(x1(1,i)+x1(1,i-1))/2;
end
B=[-1.*z1',ones(b-1,1)];
Y=x0(1,2:b)';
U=(B.'*B)^(-1)*B.'*Y;
p=U(1,1);
q=U(2,1);
m=q/p;
XR=zeros(a,b);
for i =1:b+1
    XR(1,i)=m+(x0(1,1)-m)*exp(-p*(i-1));
end
for i=1:b
    XR(1,b+2-i)=XR(1,b+2-i)-XR(1,b+1-i);
end
end

根据上式可以预测之后的结果:
$$XR=\left[\begin{array}{cccccccc}92.8100& 98.1158& 98.5400& 98.9660& 99.3938& 99.8235& 100.2550& 100.6885\end{array}\right]$$