# -*- 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()

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

相關文章