import numpy as np
import pandas as pd
from decimal import *
import matplotlib.pyplot as plt
def Grade_ratio_test(X0):
    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('該數(shù)據(jù)未通過級(jí)比檢驗(yàn)')
            return False
    print('該數(shù)據(jù)通過級(jí)比檢驗(yàn)')
    return True
def model_train(X0_train):
    #AGO生成序列X1
    X1 = X0_train.cumsum()
    Z= (np.array([-0.5 * (X1[k - 1] + X1[k]) for k in range(1, len(X1))])).reshape(len(X1) - 1, 1)
    # 數(shù)據(jù)矩陣A、B
    A = (X0_train[1:]).reshape(len(Z), 1)
    B = np.hstack((Z, np.ones(len(Z)).reshape(len(Z), 1)))
    # 求灰參數(shù)
    a, u = np.linalg.inv(np.matmul(B.T, B)).dot(B.T).dot(A)
    u = Decimal(u[0])
    a = Decimal(a[0])
    print("灰參數(shù)a：", a, "，灰參數(shù)u：", u)
    return u,a
def model_predict(u,a,k,X0):
    predict_function =lambda k: (Decimal(X0[0]) - u / a) * np.exp(-a * k) + u / a 
    X1_hat = [float(predict_function(k)) for k in range(k)]
    X0_hat = np.diff(X1_hat)
    X0_hat = np.hstack((X1_hat[0], X0_hat))
    return X0_hat
'''
根據(jù)后驗(yàn)差比及小誤差概率判斷預(yù)測結(jié)果
:param X0_hat: 預(yù)測結(jié)果
:return:
'''
def result_evaluate(X0_hat,X0):
    S1 = np.std(X0, ddof=1)  # 原始數(shù)據(jù)樣本標(biāo)準(zhǔn)差
    S2 = np.std(X0 - X0_hat, ddof=1)  # 殘差數(shù)據(jù)樣本標(biāo)準(zhǔn)差
    C = S2 / S1  # 后驗(yàn)差比
    Pe = np.mean(X0 - X0_hat)
    temp = np.abs((X0 - X0_hat - Pe)) < 0.6745 * S1    
    p = np.count_nonzero(temp) / len(X0)  # 計(jì)算小誤差概率
    print("原數(shù)據(jù)樣本標(biāo)準(zhǔn)差：", S1)
    print("殘差樣本標(biāo)準(zhǔn)差：", S2)
    print("后驗(yàn)差比：", C)
    print("小誤差概率p：", p)
if __name__ == '__main__':
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 步驟一（替換sans-serif字體）
    plt.rcParams['axes.unicode_minus'] = False  # 步驟二（解決坐標(biāo)軸負(fù)數(shù)的負(fù)號(hào)顯示問題）
        # 原始數(shù)據(jù)X
    data = pd.read_excel('./siwei_day_traffic.xlsx')
    X=data[data['week_day']=='周五'].jam_num[:5].astype(float).values
    print(X)
    # 訓(xùn)練集
    X_train = X[:int(len(X) * 0.7)]
    # 測試集
    X_test = X[int(len(X) * 0.7):]
    Grade_ratio_test(X_train)  # 判斷模型可行性
    a,u=model_train(X_train)  # 訓(xùn)練
    Y_pred = model_predict(a,u,len(X),X)  # 預(yù)測
    Y_train_pred = Y_pred[:len(X_train)]
    Y_test_pred = Y_pred[len(X_train):]
    score_test = result_evaluate(Y_test_pred, X_test)  # 評(píng)估
    # 可視化
    plt.grid()
    plt.plot(np.arange(len(X_train)), X_train, '->')
    plt.plot(np.arange(len(X_train)), Y_train_pred, '-o')
    plt.legend(['負(fù)荷實(shí)際值', '灰色預(yù)測模型預(yù)測值'])
    plt.title('訓(xùn)練集')
    plt.show()
    plt.grid()
    plt.plot(np.arange(len(X_test)), X_test, '->')
    plt.plot(np.arange(len(X_test)), Y_test_pred, '-o')
    plt.legend(['負(fù)荷實(shí)際值', '灰色預(yù)測模型預(yù)測值'])
    plt.title('測試集')
    plt.show()