灰色预测模型

灰色预测模型

应用

灰色预测模型（Gray Forecast Model）是通过少量的、不完全的信息，建立数学模型并做出预测的一种预测方法。是处理小样本（4个就可以）预测问题的有效工具，而对于小样本预测问题回归和神经网络的效果都不太理想。

灰色系统

我们称信息完全未确定的系统为黑色系统，称信息完全确定的系统为白色系统，灰色系统就是这介于这之间，一部分信息是已知的，另一部分信息是未知的，系统内各因素间有不确定的关系。

特点

• 用灰色数学处理不确定量，使之量化。
• 充分利用已知信息寻求系统的运动规律。
• 灰色系统理论能处理贫信息系统。

灰色生成数列

灰色系统理论认为，尽管客观表象复杂，但总是有整体功能的，因此必然蕴含某种内在规律。关键在于如何选择适当的方式去挖掘和利用它。灰色系统时通过对原始数据的整理来寻求其变化规律的，这是一种就数据寻求数据的现实规律的途径，也就是灰色序列的生产。一切灰色序列都能通过某种生成弱化其随机性，显现其规律性。数据生成的常用方式有累加生成、累减生成和加权累加生成。常用的是累加生成。
设原始数据为$x^{(0)}=(x^0(1),x^0(2),\dots \dots ,x^0(n))$

1.累加生成

$x^1(1)=x^0(1)$
$x^1(2)=x^0(1)+x^0(2)$
$x^1(3)=x^0(1)+x^0(2)+x^0(3)$
$\dots$
$x^1(n)=x^0(1)+x^0(2)+\dots \dots +x^0(n)$
累加的数据为$x^{(1)}=(x^1(1),x^1(2),\dots \dots ,x^1(n))$
例如有一组数据的的折线如下

不能看出数据有什么的规律，但经过累加生成后的结果如下

看起来就是一个递增的规律（这是20期某大坝变形位移的数据，顺水流方向的形变一定都是正的）

2.加权临值生成

$z^0(2)=a{x}^0(2)+(1-a)x^0(1)$
$z^0(3)=a{x}^0(3)+(1-a)x^0(2)$
$z^0(4)=a{x}^0(4)+(1-a)x^0(3)$
$\dots$
$z^1(n)=a{x}^0(n)+(1-a)x^0(n-1)$
由此得到的数列称为邻值生成数，权$a$也称为生成系数。 特别地，当生成系数$\alpha =0.5$时，则称该数列为均值生成数，也称为等权邻值生成数。

灰色模型GM(1,1)

GM代表grey model（灰色模型），GM(1,1)是一阶微分方程模型。

1.数据检验

使用GM(1,1)建模需要对数据进行检验，首先计算数列的级比
$\lambda (k)=\frac{x^0(k-1)}{x^0(k)},其中k=2,3,...,n$
如果所有的级比都落在可容覆盖区间$X=(e^{\frac{-2}{n+1}},e^{\frac{2}{n+1}})$内，则数列$x^{(}0)$可以建立GM(1,1)模型进行灰色预测。否则就需要对数据做适当的变换处理，如平移等。

2.构建灰色模型

定义$x^{(1)}$的灰导数为
$d(k)=x^0(k)=x^1(k)-x^1(k-1)$
令$z^1(k)$为数列$x^1$的邻值生成数列，即
$z^1(k)=a{x}^1(k)+(1-a)x^1(k-1)$
于是定义GM(1,1)的灰微分方程模型为
$d(k)+a{z}^1(k)=b$
其中，$a$称为发展系数，$z^1(k)$称为白化背景值，$b$称为灰作用量。接下来我们得到如下方程组
$x^0(2)+a{z}^1(2)=b$
$x^0(3)+a{z}^1(3)=b$
$...$
$x^0(n)+a{z}^1(n)=b$
按照矩阵的方法列出
$$u=\left[\begin{array}{c}a\\ b\end{array}\right],Y=\left[\begin{array}{c}{x}^0(2)\\ x^0(3)\\ \dots \\ x^0(n)\end{array}\right],B=\left[\begin{array}{cc}-z^1(2)& 1\\ -z^1(3)& 1\\ \dots \dots & \\ -z^1(n)& 1\end{array}\right]$$
则GM(1,1)就可以表示为$Y=Bu$，接下来就是求$a$和$b$的值，可以使用线性回归或者$(B^{T}B{)}^{-1}{B}^{T}Y$（正规方程）等按照最小二乘的原理来求出$a$和$b$的值（这里在列方程时$a$可以随便写一个数，比如$\frac{1}{2}$，这样$a$和$(1-a)$就一样的可以提出来方便写）。

3.预测

相应的白化模型为
$\frac{d{x}^1(t)}{dt}+a{x}^1(t)=b$
由此得到$x^1(t)的解为$
$x^1(t)=(x^0(1)-\frac{b}{a})e^{-a(t-1)}+\frac{b}{a}$
令$t+1=t$有
$x^1(t+1)=(x^0(1)-\frac{b}{a})e^{-a}+\frac{b}{a}，其中k=1,2,3\dots ,n-1$
这就是我们的预测值。

4.检验

灰色模型的精度检验一般有三种方法灰色模型的精度检验一般有三种方法，相对误差大小检验法，关联度检验法和后验差检验法。常用的为后验差检验法。

1. 将预测的${\hat{x}}^1$使用累减生成得到${\hat{x}}^0$
   ${\hat{x}}^0(k)={\hat{x}}^1(k)-{\hat{x}}^1(k-1),其中k=2,3,...,n$
2. 计算残差
   $e(k)=x^0(k)-{\hat{x}}^0(k),其中k=1,2,\dots ,n$
3. 计算原始序列$x^0$的方差$S_1$和残差$e$的方差$S_2$
   $S_1=\frac{1}{n}{\sum }_{k=1}^n(x^0(k)-\bar{x}{)}^2$
   $S_2=\frac{1}{n}{\sum }_{k=1}^n(e(k)-\bar{e}{)}^2$
4. 计算后验差比
   $C=\frac{S_2}{S_1}$
5. 查表观察效果

模型精度等级	均方差比值C
1 级(好)	C<=0.35
2 级（合格）	C<=0.5&c>0.35
3 级（勉强）	C<=0.65&c>0.5
4 级（不合格）	C>0.65

python代码

这里我们使用的网上一个例子，数据是几年的噪声数据，如下

以下是python代码

# 灰色模型预测

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

dir = 'C:/Users/Administrator/Desktop/'

data = pd.read_excel(dir+'test.xls', sheetname='Sheet1')
data = np.array(data['db'])
lens = len(data)  # 数据量

# 数据检验
## 计算级比
lambds = []
for i in range(1, lens):
    lambds.append(data[i-1]/data[i])
## 计算区间
X_min = np.e**(-2/(lens+1))
X_max = np.e**(2/(lens+1))
## 检验
is_ok = True
for lambd in lambds:
    if (lambd < X_min or lambd > X_max):
        is_ok = False
if (is_ok == False):
    print('该数据未通过检验')
else:
    print('该数据通过检验')

# 构建灰色模型GM(1,1)
## 累加数列
data_1 = []
data_1.append(data[0])
for i in range(1, lens):
    data_1.append(data_1[i-1]+data[i])
## 灰导数及临值生成数列
ds = []
zs = []
for i in range(1, lens):
    ds.append(data[i])
    zs.append(-1/2*(data_1[i-1]+data_1[i]))
## 求a、b
B = np.array(zs).reshape(lens-1,1)
one = np.ones(lens-1)
B = np.c_[B, one]  # 加上一列1
Y = np.array(ds).reshape(lens-1,1)
a, b = np.dot(np.dot(np.linalg.inv(np.dot(B.T, B)), B.T), Y)
print('a='+str(a))
print('b='+str(b))

# 预测
def func(k):
    c = b/a
    return (data[0]-c)*(np.e**(-a*k))+c
data_1_hat = []  # 累加预测值
data_0_hat = []  # 原始预测值
data_1_hat.append(func(0))
data_0_hat.append(data_1_hat[0])
for i in range(1, lens+5):  # 多预测5次
    data_1_hat.append(func(i))
    data_0_hat.append(data_1_hat[i]-data_1_hat[i-1])
print('预测值为：')
for i in data_0_hat:
    print(i)

# 模型检验
## 预测结果方差
data_h = np.array(data_0_hat[0:7]).T
sum_h = data_h.sum()
mean_h = sum_h/lens
S1 = np.sum((data_h-mean_h)**2)/lens
## 残差方差
e = data - data_h
sum_e = e.sum()
mean_e = sum_e/lens
S2 = np.sum((e-mean_e)**2)/lens
## 后验差比
C = S2/S1
## 结果
if (C <= 0.35):
    print('1级，效果好')
elif (C <= 0.5 and C >= 0.35):
    print('2级，效果合格')
elif (C <= 0.65 and C >= 0.5):
    print('3级，效果勉强')
else:
    print('4级，效果不合格')

# 画图
plt.figure(figsize=(9, 4), dpi=100)
x1 = np.linspace(1, 7, 7)
x2 = np.linspace(1, 12, 12)
plt.subplot(121)
plt.title('x^0')
plt.plot(x2, data_0_hat, 'r--', marker='*')
plt.scatter(x1, data, marker='^')
plt.subplot(122)
plt.title('x^1')
plt.plot(x2, data_1_hat, 'r--', marker='*')
plt.scatter(x1, data_1, marker='^')
plt.show()

输出结果如下

• 该数据通过检验
• a=[0.00234379]
  b=[72.6572696]
• 预测值为：
  [71.1]
  [72.40574144]
  [72.23623656]
  [72.0671285]
  [71.89841633]
  [71.73009912]
  [71.56217595]
  [71.39464589]
  [71.22750803]
  [71.06076145]
  [70.89440522]
  [70.72843845]
• 1级，效果好
•