前言

灰色预测是一种对含有不确定因素的系统进行预测的方法，有助于解决样本数量少、数据信息少和信息不完整等不确 定的问题。灰色预测 ＧＭ（１，１）模型可利用随机原始的离散非负数据列，通过一次累加的方式得到新的离散数据列，建立相应的微分方程模型，用一阶线性微分方程的解来逼近时间累加后形成的新时间序列所呈现的规律。 它通过一个序列即可进行计算，是一阶线性的动态模型。

一、灰色模型的建模步骤及原理

1.灰色模型的基本概念

基本定义：灰色系统是介于黑色模型及白色模型之间的,既包含未知信息又包含己知信息同时系统内各因素具有不确定性关系。
一般表达式：GM(a,b)，用a 阶微分方程去对 b个变量进行求解。
基本求解原理：通过对原始数据数列进行处理,生成相关性较强的数据数列，同时建立相关的微分方程进而预测其未来的发展趋势。

2.灰色模型的建模步骤

把原始数据加工成相关数列；
利用生成的相关数列进行求解，得到相应的微分方程参数；
利用得到的参数进行求解预测值；
对预测值进行相关检验。

二、模型建立

1.确定原始数据

首先确定一组原始数据如下:

其中 n 为未知数的个数。

2.累加数列和邻均值等权数列

对原始数据进行累加得到新数据：

即每一个对应的数值为原始数据序列对应位置的值以及之前的数值的累加和，同时生成一个新数列的邻均值等权数列如下：

3.建立关于t 的白化形式一阶一元微分方程 GM(1,1)

对于一组单调递增的数据，累加后的数据可以看作是强烈单调增的数据，可以近似认为符合指数增长。对于指数曲线，它一定是某个一阶线性常系数微分方程，于是可以写出如下表达式

其中 a，b是我们所要求解的系数，分别为发展系数和灰色作用量，a 的有效区间是（-2，2）。

三、模型求解

1.利用矩阵求解参数(最小二乘法)

根据最小二乘法求解得到：

2.对一阶微分方程中进行求解

由白化方程可以得到：


将y和C代回ln y = -at+C，最终得到该模型的时间响应式，为：

四、模型检验

之后从3项指标检验模型的预测精度，分别是平均相对误差、精度和方差比，计算公式如下：

均方差比表明预测值与实际值误差的大小，其值越小则误差越小。

此外，还需计算小概率误差，小误差概率是指即残差与残差均值的差值的绝对值小于0.6745倍原始数据均方差的概率，其值大，则表示精度高。计算方式如下:

根据计算结果对比以下GM（1,1）模型预测精度等级表得到检验结果。

五、数据检验方法

两种检验一般在建模前进行，二选一即可，也可以都进行。

1.光滑比检验

判断数据是否为准光滑序列，采用光滑比分析原始数据光滑性，数据变化越平稳，光滑比越小，计算方法如下：

若数据序列满足以下两个条件即可判定为准光滑序列，可用于 GM(1,1)预测模型。

2.级比检验

若级别检验后发现不通过，则需要添加常数项，使其落入区间内，即有

总结

微分方程模型作为灰色系统理论预测用的模型，其主要凭借以下几点：
(1)灰色系统理论把原始数列在固定范围变化的灰色量，认定随机过程是在一定时区一定范围内变化的灰色过程。
(2) 一阶微分方程有着指数增长形式的解值，无规律的原始数列经过累加生成得到累加数列后，此数列也具备指数增加的规律。灰色模型是生成的数据数列所建模型。
(3) 灰色系统理论可建立 GM(1，n)模型来解决高阶系统建模问题。一阶微分方程所组成的灰色模型称为 GM 模型群。
(4) GM模型采用累加生成做生成处置，为了获得预测的数据，由 GM 模型所得的数据对其进行逆生成处理，即累减生成还原处理。
(5) 在灰色系统理论中，为了提高预测精度，可以对数据进行合适的取舍、对灰数采用不同的生成方式和调节残差的级别来实现。

代码实现

构建数列进行光滑比检验

import numpy as np
X0 = np.array(['数据'])
X1 = X0.cumsum()
rho = [X0[i] / X1[i - 1] for i in range(1,len(X0))]
rho_ratio = [rho[i + 1] / rho[i] for i in range(len(rho) - 1)]
print("rho:",rho)
print("rho_ratio:",rho_ratio)
flag = True
for i in range(1,len(rho) - 1):
    if rho[i] > 0.5 or rho[i + 1] / rho[i] >= 1:
        flag = False
if rho[-1] > 0.5:
    flag = False
if flag:
    print("数据通过光滑校验")
else:
    print("该数据未通过光滑校验")

构建数列进行级比检验

import math as mt
X0 = ['数据']
for i in range(len(X0)-1):
    l = X0[i]/X0[i+1]
    if l <= mt.exp(-2/(len(X0)+1)) or l >= mt.exp(2/(len(X0)+1)):
        break
    else:
        pass
if i == len(X0)-2 and l > mt.exp(-2/(len(X0)+1)) and l < mt.exp(2/(len(X0)+1)):
    print('级比检验通过')
else:
    print('级比检验不通过')

级比检验不通过处理

import math as mt
X0 = [88898,97242,97990,93363,88093,103852,105475,114183,133762,162214,173977,200973]
j = 1
while True:
    YO = [k+j for k in X0]
    j += 1
    for m in range(len(YO) - 1):
        l = YO[m] / YO[m + 1]
        if l > mt.exp(-2 / (len(X0) + 1)) and l < mt.exp(2 / (len(X0) + 1)):
            b = True
        else:
            b = False
            break
    if b == True:
        print("新的原始数列为：",YO)
        c = j -1
        print("c的值为：",c)
        break
    else:
        continue

构建模型并求解

import numpy as np
import math as mt

X0 = ['数据']
# 累加数列
X1 = [X0[0]]
add = X0[0] + X0[1]
X1.append(add)
i = 2
while i < len(X0):
    add = add + X0[i]
    X1.append(add)
    i += 1
    
# 紧邻均值序列
Z = []
j = 1
while j < len(X1):
    num = (X1[j] + X1[j - 1]) / 2
    Z.append(num)
    j = j + 1

# 最小二乘法计算
Y = []
x_i = 0
while x_i < len(X0) - 1:
    x_i += 1
    Y.append(X0[x_i])
Y = np.mat(Y)
Y = Y.reshape(-1,1)
B = []
b = 0
while b < len(Z):
    B.append(-Z[b])
    b += 1
B = np.mat(B)
B = B.reshape(-1,1)
c = np.ones((len(B),1))
B = np.hstack((B,c))
print("B",B)

# 求出参数
alpha = np.linalg.inv(B.T.dot(B)).dot(B.T).dot(Y)
a = alpha[0,0]
b = alpha[1,0]
print('alpha',alpha)
print("a=",a)
print("b=",b)

# 生成预测模型
GM = []
GM.append(X0[0])
did = b/a
k = 1
while k < len(X0):
    GM.append((X0[0] - did) * mt.exp(-a * k) + did)
    k += 1

# 做差得到预测序列
G = []
G.append(X0[0])
g = 1
while g < len(X0):
    G.append(round(GM[g] - GM[g - 1]))
    g += 1
print("预测数列为：",G)

绘图

import matplotlib.pyplot  as plt
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False 

X0 = ['原始数据']
G = ['预测数据']
G = [g-37315 for g in G]
r = range(len(X0))
t = list(r)
plt.plot(t,X0,color='r',linestyle="--",label='true')
plt.plot(t,G,color='b',linestyle="--",label="predict")
plt.legend()
plt.show()

模型检验

import numpy as np
import math as mt
X0 = ['原始数据']
G = ['预测数据']
X0 = np.array(X0)
G = np.array(G)
e = X0 - G  #残差
q = e / X0  # 相对误差
w = 0
for q_i in q:
    w += q_i
w = w/len(X0)
print('精度为{}%'.format(round((1-w)*100,2)))

s0 = np.var(X0)
s1 = np.var(e)
S0 = mt.sqrt(s0)
S1 = mt.sqrt(s1)
C = S1 / S0
print('方差比为：',C)

p = 0
for s in range(len(e)):
    if (abs(e[s]-np.mean(e)) < 0.6745 * S1):
        p = p + 1
P = p / len(e)
print('小概率误差为：',P)