灰色預測改進—三角殘差擬合_python

# -*- coding: utf-8 -*-

# @Time : 2022/3/30 14:44

# @Author : Orange

# @File : gm_trigonometric.py

from decimal import *

import math

class GM_trig():

     def __init__(self, X0, X0_hat, L):

         '''

         :param X0: 原始序列

         :param X0_hat: 採用 GM(1,1) 後獲得的序列

         :param L: 使用者自定義的迴圈週期分量

         '''

         self.X0 = X0

         self.X0_hat = X0_hat

         self.R0 = (X0 - X0_hat)[1:]  # 殘差

         self.L = L

         self.n = len(self.X0)

         self.B = None

     def train(self):

         self.B = np.array(

             [[1] * (self.n - 1), np.arange(1, self.n), [np.sin(2 * i * math.pi / self.L) for i in range(1, self.n)],

              [np.cos(2 * i * math.pi / self.L) for i in range(1, self.n)]]).T

         R_n = np.array(self.R0).reshape(self.n - 1, 1)

         b_hat = np.linalg.inv(np.matmul(self.B.T, self.B)).dot(self.B.T).dot(R_n)

         self.f_R0 = lambda k: b_hat[0][0] + b_hat[1][0] * k + b_hat[2][0] * np.sin(2 * k * math.pi / self.L) + b_hat[3][

             0] * np.cos(2 * k * math.pi / self.L)  

     def predict(self, k, X_all_0_hat):

         '''

         :param k: 給出從 0 ， k 的預測值

         :param X_all_0_hat: 所有資料（訓練 + 測試）的 GM(1,1) 預測值列表

         :return:

         '''

         R0_hat = [self.f_R0(k) for k in range(1, k)]

         R0_hat.insert(0, 0)

         X_tr_hat = X_all_0_hat + R0_hat

         return X_tr_hat

     def interval_pred(self, X_tr_hat, t_val, k):

         s = math.sqrt(sum((X_tr_hat[:self.n] - self.X0) ** 2) / (self.n - 6))

         l_bound = []

         h_bound = []

         for k in range(1, k):

             Y_k = np.array([1, k, np.sin(2 * k * math.pi / self.L), np.cos(2 * k * math.pi / self.L)]).reshape(4, 1)

             h_kk = Y_k.T.dot(np.linalg.inv(np.matmul(self.B.T, self.B))).dot(Y_k)[0][0]

             ll_k = X_tr_hat[k] - t_val * s * math.sqrt(1 + h_kk)

             hh_k = X_tr_hat[k] + t_val * s * math.sqrt(1 + h_kk)

             l_bound.append(ll_k)

             h_bound.append(hh_k)

         return l_bound, h_bound

class GM11():

     def __init__(self):

         self.f = None

     def isUsable(self, X0):

         ''' 判斷是否透過光滑檢驗 '''

         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(2, 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(" 該資料未透過光滑校驗 ")

         ''' 判斷是否透過級比檢驗 '''

         lambds = [X0[i - 1] / X0[i] for i in range(1, len(X0))]

         X_min = np.e ** (-2 / (len(X0) + 1))

         X_max = np.e ** (2 / (len(X0) + 1))

         for lambd in lambds:

             if lambd < X_min or lambd > X_max:

                 print(' 該資料未透過級比檢驗 ')

                 return

         print(' 該資料透過級比檢驗 ')

     def train(self, X0):

         X1 = X0.cumsum()

         Z = (np.array([-0.5 * (X1[k - 1] + X1[k]) for k in range(1, len(X1))])).reshape(len(X1) - 1, 1)

         # 資料矩陣 A 、 B

         A = (X0[1:]).reshape(len(Z), 1)

         B = np.hstack((Z, np.ones(len(Z)).reshape(len(Z), 1)))

         # 求灰引數

         a, u = np.linalg.inv(np.matmul(B.T, B)).dot(B.T).dot(A)

         u = Decimal(u[0])

         a = Decimal(a[0])

         print(" 灰引數 a ： ", a, " ，灰引數 u ： ", u)

         self.f = 外匯跟單gendan5.comlambda k: (Decimal(X0[0]) - u / a) * np.exp(-a * k) + u / a

     def predict(self, k):

         X1_hat = [float(self.f(k)) for k in range(k)]

         X0_hat = np.diff(X1_hat)

         X0_hat = np.hstack((X1_hat[0], X0_hat))

         return X0_hat

     def evaluate(self, X0_hat, X0):

         '''

         根據後驗差比及小誤差機率判斷預測結果

         :param X0_hat: 預測結果

         :return:

         '''

         S1 = np.std(X0, ddof=1)  # 原始資料樣本標準差

         S2 = np.std(X0 - X0_hat, ddof=1)  # 殘差資料樣本標準差

         C = S2 / S1  # 後驗差比

         Pe = np.mean(X0 - X0_hat)

         temp = np.abs((X0 - X0_hat - Pe)) < 0.6745 * S1

         p = np.count_nonzero(temp) / len(X0)  # 計算小誤差機率

         print(" 原資料樣本標準差： ", S1)

         print(" 殘差樣本標準差： ", S2)

         print(" 後驗差比： ", C)

         print(" 小誤差機率 p ： ", p)

def MAPE(y_true, y_pred):

     """ 計算 MAPE"""

     n = len(y_true)

     mape = sum(np.abs((y_true - y_pred) / y_true)) / n * 100

     return mape

if __name__ == '__main__':

     import matplotlib.pyplot as plt

     import numpy as np

     import pandas as pd

     plt.rcParams['font.sans-serif'] = ['SimHei']  # 步驟一（替換 sans-serif 字型）

     plt.rcParams['axes.unicode_minus'] = False  # 步驟二（解決座標軸負數的負號顯示問題）

     # 原始資料 X

     data = pd.read_csv("test.csv")

     X = data["val"].values

     # 訓練集

     X_train = X[:-4]

     # 測試集

     X_test = X[-4:]

     model = GM11()

     model.isUsable(X_train)  # 判斷模型可行性

     model.train(X_train)  # 訓練

     Y_pred = model.predict(len(X))  # 預測

     Y_train_pred = Y_pred[:len(X_train)]

     Y_test_pred = Y_pred[len(X_train):]

     score_test = model.evaluate(Y_test_pred, X_test)  # 評估

     print("gm(1,1)_mape:", MAPE(Y_train_pred, X_train), "%")

     model_trig = GM_trig(X_train, Y_train_pred, L=23)

     model_trig.train()

     result = model_trig.predict(len(X), Y_pred)

     X_train_pred = result[:-4]

     X_test_pred = result[-4:]

     l_bound, h_bound = model_trig.interval_pred(result, 2.179, len(X))

     # 視覺化

     plt.grid()

     plt.plot(np.arange(len(X_train)), X_train, '->')

     plt.plot(np.arange(len(X_train)), X_train_pred, '-o')

     plt.plot(np.arange(len(X_train)), Y_train_pred, '-*')

     plt.legend([' 負荷實際值 ', ' 三角殘差預測值 ', 'GM(1,1) 預測值 '])

     print("gm(1,1)_trig_mape:", MAPE(X_train_pred, X_train), "%")

     plt.title(' 訓練集 ')

     plt.show()

     # 視覺化

     plt.grid()

     plt.plot(np.arange(len(X_test)), X_test, '->')

     plt.plot(np.arange(len(X_test)), X_test_pred, '-o')

     plt.plot(np.arange(len(X_test)), Y_test_pred, '-*')

     plt.legend([' 負荷實際值 ', ' 三角殘差預測值 ', 'GM(1,1) 預測值 '])

     plt.title(' 測試集 ')

     plt.show()

     # 區間預測視覺化

     plt.figure(figsize = (10, 4))

     plt.grid(axis='y')

     plt.plot(np.arange(1, 22), h_bound,'k--.')

     plt.plot(np.arange(1, 22), l_bound,'k-->')

     plt.plot(np.arange(22), X,'r-o')

     plt.plot(np.arange(22), result,'b-d')

     plt.legend([' 上界 ', ' 下界 ', ' 真實值 ',' 預測值 '])

     plt.title(' 灰度預測——區間預測 ')

     plt.show()