python 實現 灰色預測 GM(1,1)模型 灰色系統 預測 灰色預測公式推導

 來源公式推導連接

　　https://blog.csdn.net/qq_36387683/article/details/88554434

 

關鍵詞：灰色預測 python 實現 灰色預測 GM(1,1)模型 灰色系統 預測 灰色預測公式推導

一、前言

  本文的目的是用Python和類對灰色預測進行封裝

二、原理簡述

1.灰色預測概述

  灰色預測是用灰色模型GM(1,1)來進行定量分析的，通常分為以下幾類：
    (1) 灰色時間序列預測。用等時距觀測到的反映預測對象特征的一系列數量（如產量、銷量、人口數量、存款數量、利率等）構造灰色預測模型，預測未來某一時刻的特征量，或者達到某特征量的時間。
    (2) 畸變預測（災變預測）。通過模型預測異常值出現的時刻，預測異常值什么時候出現在特定時區內。
    (3) 波形預測，或稱為拓撲預測，它是通過灰色模型預測事物未來變動的軌跡。
    (4) 系統預測，對系統行為特征指標建立一族相互關聯的灰色預測理論模型，在預測系統整體變化的同時，預測系統各個環節的變化。
  上述灰色預測方法的共同特點是：
    （1）允許少數據預測；
    （2）允許對灰因果律事件進行預測，例如：
      灰因白果律事件：在糧食生產預測中，影響糧食生產的因子很多，多到無法枚舉，故為灰因，然而糧食產量卻是具體的，故為白果。糧食預測即為灰因白果律事件預測。
      白因灰果律事件：在開發項目前景預測時，開發項目的投入是具體的，為白因，而項目的效益暫時不很清楚，為灰果。項目前景預測即為灰因白果律事件預測。
    （3）具有可檢驗性，包括：建模可行性的級比檢驗（事前檢驗），建模精度檢驗（模型檢驗），預測的滾動檢驗（預測檢驗）。

 

2.GM(1,1)模型理論

  GM(1,1)模型適合具有較強的指數規律的數列，只能描述單調的變化過程。已知元素序列數據：

做一次累加生成（1-AGO）序列：

其中

，

令為的緊鄰均值生成序列：

其中，

建立GM(1,1)的灰微分方程模型為：

其中，為發展系數，為灰色作用量。設為待估參數向量，即，則灰微分方程的最小二乘估計參數列滿足

其中

再建立灰色微分方程的白化方程（也叫影子方程）：

白化方程的解（也叫時間響應函數）為

那么相應的GM(1,1)灰色微分方程的時間響應序列為：

取，

則

再做累減還原可得

即為預測方程。

注1：原始序列數據不一定要全部使用，相應建立的模型也會不同，即和不同；

注2：原始序列數據必須要等時間間隔、不間斷。

3.算法步驟

(1) 數據的級比檢驗
  為了保證灰色預測的可行性，需要對原始序列數據進行級比檢驗。
  對原始數據列，

計算序列的級比：  

若所有的級比都落在可容覆蓋內，則可進行灰色預測；否則需要對做平移變換，，使得滿足級比要求。

(2) 建立GM(1,1)模型，計算出預測值列。

(3) 檢驗預測值：

① 相對殘差檢驗，計算

  

若 ，則認為達到一般要求，若 ，則認為達到較高要求；

② 級比偏差值檢驗

  根據前面計算出來的級比, 和發展系數, 計算相應的級比偏差：

  

若, 則認為達到一般要求，若, 則認為達到較高要求。

(4) 利用模型進行預測。

三、程序實現 python實現

　　需要安裝numpy和Torch加速等

 

# condig:utf-8
import torch as th
import numpy as np
class GM():

    def __init__(self):
        # 判斷是否可用 gpu 編程 , 大量級計算使用GPU
        self._is_gpu = False  # th.cuda.is_available()

    def fit(self,dt:list or np.ndarray):
        self._df :th.Tensor = th.from_numpy(np.array(dt,dtype=np.float32))

        if self._is_gpu:
            self._df.cuda()

        self._n:int = len(self._df)

        self._x,self._max_value = self._sigmod(self._df)

        z:th.Tensor = self._next_to_mean(th.cumsum(self._x,dim=0))

        self.coef:th.Tensor = self._coefficient(self._x, z)

        del z

        self._x0:th.Tensor = self._x[0]

        self._pre:th.Tensor = self._pred()

    # 歸一化
    def _sigmod(self,x:th.Tensor):
        _maxv:th.Tensor = th.max(x)
        return th.div(x,_maxv),_maxv

    # 計算緊鄰均值數列
    def _next_to_mean(self, x_1:th.Tensor):

        z:th.Tensor = th.zeros(self._n-1)
        if self._is_gpu:
            z.cuda()

        for i in range(1,self._n):  # 下標從0開始，取不到最大值
            z[i - 1] = 0.5 * x_1[i] + 0.5 * x_1[i - 1]

        return z

    # 計算系數 a,b
    def _coefficient(self,x:th.Tensor,z:th.Tensor):

        B:th.Tensor = th.stack((-1*z, th.ones(self._n-1)),dim=1)

        Y:th.Tensor = th.tensor(x[1:],dtype=th.float32).reshape((-1,1))

        if self._is_gpu:
            B.cuda()
            Y.cuda()

        # 返回的是a和b的向量轉置，第一個是a 第二個是b；
        return th.matmul(th.matmul(th.inverse(th.matmul(B.t(), B)), B.t()),Y)


    def _pred(self,start:int=1,end:int=0):

        les:int = self._n+end

        resut:th.Tensor = th.zeros(les)

        if self._is_gpu:
            resut.cuda()

        resut[0] = self._x0

        for i in range(start,les):
            resut[i] = (self._x0 - (self.coef[1] / self.coef[0])) * \
                            (1 - th.exp(self.coef[0])) * th.exp(-1 * self.coef[0] * (i))
        del les
        return resut

    # 計算絕對誤差
    def confidence(self):
        return round((th.sum(th.abs(th.div((self._x-self._pre),self._x)))/self._n).item(),4)

    # 預測個數，默認個數大於等於0，
    def predict(self,m:int=1,decimals:int=4):

        y_pred:th.Tensor = th.mul(self._pre,self._max_value)

        y_pred_ = th.zeros(1)

        if m<0:
            return "預測個數需大於等於0"
        elif m>0:
            y_pred_:th.Tensor = self._pred(self._n,m)[-m:].mul(self._max_value)
        else:
            if self._is_gpu:
                return list(map(lambda _: round(_, decimals), y_pred.cpu().numpy().tolist()))
            else:
                return list(map(lambda _:round(_,decimals),y_pred.numpy().tolist()))

        # cat 拼接 0 x水平拼接，1y垂直拼接
        result:th.Tensor = th.cat((y_pred,y_pred_),dim=0)

        del y_pred,y_pred_

        if self._is_gpu:
            return list(map(lambda _: round(_, decimals), result.cpu().numpy().tolist()))

        return list(map(lambda _:round(_,decimals),result.numpy().tolist()))



if __name__=="__main__":
    ls = np.arange(91,100,2)
    print(type(ls))
    # ls = list(range(91, 100, 2))
    gm = GM()
    gm.fit(ls)
    print(gm.confidence())
    print(ls)
    print(gm.predict(m=2))