分类： Python/Ruby
2022-05-06 17:49:29
# -*- 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 = lambda 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'] =外汇跟单gendan5.com 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()
