数学建模 | 灰色预测原理及python实现

一、灰色预测的原理

二、灰色预测的应用及python实现

一、灰色预测的原理

灰色预测是以灰色模型为基础，灰色模型GM(n,h)是微分方程模型，可用于描述对象做长期、连续、动态的反应。其中，n代表微分方程式的阶数，h代表微分方程式的变化数目。在诸多的灰色模型中，以灰色系统中单序列一阶线性微分方程模型GM(1,1)最为常见。

下面说明GM(1,1)模式：

设有原始数据列 $X^{(0)}=\left\{x^{(0)}(1),x^{(0)}(2),...,x^{(0)}(n)\right\}$ ，n为数据个数，则可以根据以下步骤来建立GM(1,1)模型：

步骤一：

与原来统计累加以便减弱随机数据序列的波动性与随机性，从而得到新数据序列：

$X^{(1)}=\left\{x^{(1)}(1),x^{(1)}(2),...,x^{(1)}(n)\right\}$

其中， $x^{(1)}(k)$ 中各数据表示前几项的相加，即：

$x^{(1)}(k)=\sum_{i=1}^{k}x^{(0)}(i),k=1,2,...n$

步骤二：

令 $Z^{(1)}$ 为 $X^{(1)}$ 的紧邻均值生成序列：

$Z^{(1)}=\left\{z^{(1)}(2),z^{(1)}(3),...,z^{(1)}(n)\right\}$

其中， $z^{(1)}(k)=0.5x^{(1)}(k)+0.5x^{(1)}(k-1)$

步骤三：

建立GM(1,1) 的灰微分方程为：

$x^{(0)}(k)+az^{(1)}(k)=b$

灰微分方程的白化方程为：

$\tfrac{dx^{(1)}(t)}{dt}+ax^{(1)}(t)=b$

其中，a为发展系数，b为内生控制系数。

步骤四：

模型求解：构造矩阵B和向量Y：

$B=\begin{bmatrix} -z^{(1)}(2) & 1\\ -z^{(1)}(3) &1 \\ ... & ...\\-z^{(1)}(n) &1 \end{bmatrix}$ ， $Y=\begin{bmatrix} x^{(0)}(2)\\ x^{(0)}(3) \\... \\x^{(0)}(n) \end{bmatrix}$

步骤五：

设 $\hat{a}$ 为待估参数向量，即：

$\hat{a}=[ab]^{T}=(B^{T}B)^{-1}B^{T}Y$

这是利用正规方程得到的闭式解。

步骤六：

求解白化方程，可得到灰色预测的离散时间响应函数：

$\hat{x}^{(1)}(t)=(x^{(1)}(0)-\frac{b}{a})e^{-ak}+\frac{b}{a}$

那么相应的时间相应序列为：

$\hat{x}^{(1)}(k+1)=(x^{(1)}(0)-\frac{b}{a})e^{-ak}+\frac{b}{a}$

取 $x^{(1)}(0)=x^{(0)}(1)$ ，则所得的累加预测值为：

$\hat{x}^{(1)}(k+1)=(x^{(0)}(1)-\frac{b}{a})e^{-ak}+\frac{b}{a},k=1,2,...,n-1$

将预测值还原为：

$\hat{x}^{(0)}(k+1)=\hat{x}^{(1)}(k+1)-\hat{x}^{(1)}(k)=(x^{(0)}(1)-\frac{b}{a})(e^{-ak}-e^{-a(k-1)})=(x^{(0)}(1)-\frac{b}{a})(1-e^{a})e^{-ak}$

二、灰色预测的应用及python实现

某公司根据2015-2020年的产品的销售额，试构建GM(1,1)预测模型，并预测2021年的产品销售额。

原始数据为： $X^{(0)}=(2.67,3.13,3.25,3.36,3.56,3.72)$

Python实现代码：

import numpy as npclass GM:def __init__(self,N,n,data):''':param N: the number of data:param n: the number of data that is needs to be predicted:param data: time series'''self.N=Nself.n=nself.data=datadef prediction(self,a,b):''':param a: parameter a of the grey prediction model:param b: parameter b of the grey prediction model:return: a list of prediction'''# a list to save n preditionpre_data=[]# calculating the predictionfor i in range(self.n):pre=(self.data[0]-b/a)*(1-np.exp(a))*np.exp(-a*(self.N+i))pre_data.append(pre[0])return pre_datadef residual_test(self,a,b):''':param a: parameter a of the grey prediction model:param b: parameter b of the grey prediction model'''# prediction od raw datapre_rawdata=self.sequence_prediction(a,b)# calculating absolute residualabs_residual=[]for i in range(self.N):r=abs(pre_rawdata[i]-self.data[i])abs_residual.append(r)# calculating relative residualrel_residual=[]for i in range(self.N):rel=abs_residual[i]/abs(self.data[0])rel_residual.append(rel)# calculating average relative residualavg_residual=0for i in range(self.N):avg_residual=avg_residual+rel_residual[i]avg_residual=avg_residual/self.Nprint("average relative residual: {}".format(avg_residual[0]))if avg_residual<0.01:print("model accuracy: excellent(Level I)")elif avg_residual<0.05:print("model accuracy: qualified(LevelⅡ)")elif avg_residual<0.10:print("model accuracy: barely qualified(Level Ⅲ)")else:print("model accuracy: unqualified(Level Ⅳ)")def sequence_prediction(self,a,b):''':param a: parameter a of the grey prediction model:param b: parameter b of the grey prediction model:return:  prediction of raw data'''pre_rawdata=[]pre_rawdata.append(self.data[0])for i in range(1,self.N):pre_raw=(self.data[0]-b/a)*(1-np.exp(a))*np.exp(-a*(i))pre_rawdata.append(pre_raw)return pre_rawdatadef GM11(self):''':return: n prediction if the grey prediction model can be usedand parameter a and parameter b'''# accumulate raw datacumdata=[]for i in range(self.N):s=0for j in range(i+1):s=s+self.data[j]cumdata.append(s)# calculating the nearest neighbor mean generation sequenceZ=[]  # len(Z)=N-1for i in range(1,self.N):z=0.5*cumdata[i]+0.5*cumdata[i-1]Z.append(z)# construct data matrix B and data vector YB=np.array([[-Z[i],1] for i in range(self.N-1)])Y=np.array([[self.data[i]] for i in range(1,self.N)])# calculating parameter a and parameter bA=np.dot(np.dot(np.linalg.inv(np.dot(B.T,B)),B.T),Y)a=A[0]b=A[1]# determine whether the grey prediction model can be usedif (-1)*a>1:print("the grey prediction model can not be used")else:#print("a={}".format(a[0]))#print("b={}".format(b[0]))# derive a prediction model and predictpre_data=self.prediction(a,b)return pre_data,a,bif __name__=="__main__":data=[2.67,3.13,3.25,3.36,3.56,3.72]N=6n=1# create objectgrey_prediction=GM(N,n,data)# get predictionpre_data,a,b=grey_prediction.GM11()# get average relative residualavg_residual=grey_prediction.residual_test(a,b)print("predicted data: {}".format(pre_data[0]))

运行结果：