数学建模 | 灰色预测原理及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 np

class 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=N
        self.n=n
        self.data=data


    def 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 predition
        pre_data=[]

        # calculating the prediction
        for 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_data


    def 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 data
        pre_rawdata=self.sequence_prediction(a,b)

        # calculating absolute residual
        abs_residual=[]
        for i in range(self.N):
            r=abs(pre_rawdata[i]-self.data[i])
            abs_residual.append(r)

        # calculating relative residual
        rel_residual=[]
        for i in range(self.N):
            rel=abs_residual[i]/abs(self.data[0])
            rel_residual.append(rel)

        # calculating average relative residual
        avg_residual=0
        for i in range(self.N):
            avg_residual=avg_residual+rel_residual[i]
        avg_residual=avg_residual/self.N
        print("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_rawdata


    def GM11(self):
        '''
        :return: n prediction if the grey prediction model can be used
                 and parameter a and parameter b
        '''

        # accumulate raw data
        cumdata=[]
        for i in range(self.N):
            s=0
            for j in range(i+1):
                s=s+self.data[j]
            cumdata.append(s)

        # calculating the nearest neighbor mean generation sequence
        Z=[]  # len(Z)=N-1
        for 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 Y
        B=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 b
        A=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 used
        if (-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 predict
            pre_data=self.prediction(a,b)
            return pre_data,a,b


if __name__=="__main__":
    data=[2.67,3.13,3.25,3.36,3.56,3.72]
    N=6
    n=1
    # create object
    grey_prediction=GM(N,n,data)
    # get prediction
    pre_data,a,b=grey_prediction.GM11()
    # get average relative residual
    avg_residual=grey_prediction.residual_test(a,b)
    print("predicted data: {}".format(pre_data[0]))

运行结果：