灰色预测(Python)

灰色系统

        用颜色深浅来反映信息量的多少，黑色表示信息量太少，白色表示这个系统是清楚的、信息量充足。处于黑白之间的，即信息不完全的系统，称灰色系统。

灰色预测

        对既含有已知信息又含有不确定信息的系统进行预测，就是对再一定范围内变化的、与时间有关的灰色过程的预测。

        适用情况

•         数据是以年份度量的非负数据（如果是月份或者季度数据一般要用时间序列模型）
•         数据能经过准指数规律的检验
•         数据的期数较短且和其他数据之间关联性不强

        GM（1，1）模型

                只含有一个变量的一阶微分方程模型

        步骤

•                 根据原始离散非负数据列，通过累加等方式削弱随机性、获得有规律的离散数据列
•                 建立相应的微分方程模型，得到离散点处的解
•                 再通过累减求得的原始数据的估计值，从而对原始数据预测

 Python代码

import numpy as np
import math

# 历史数据
load_history = [724.57, 746.62, 778.27, 800.8, 827.75, 871.1, 912.37, 954.28, 995.01, 1037.2]
n = len(load_history)
load_array = np.array(load_history)

# 累加生成
cumulative_load = [sum(load_history[0:i + 1]) for i in range(n)]
cumulative_array = np.array(cumulative_load)

# 计算数据矩阵B和数据向量Y
B_matrix = np.zeros([n - 1, 2])
Y_vector = np.zeros([n - 1, 1])
for i in range(0, n - 1):
    B_matrix[i][0] = -0.5 * (cumulative_array[i] + cumulative_array[i + 1])
    B_matrix[i][1] = 1
    Y_vector[i][0] = load_array[i + 1]

# 计算GM(1,1)微分方程的参数a和u
A_matrix = np.linalg.inv(B_matrix.T.dot(B_matrix)).dot(B_matrix.T).dot(Y_vector)
a = A_matrix[0][0]
u = A_matrix[1][0]

# 建立灰色预测模型
predicted_load = np.zeros(n)
predicted_load[0] = load_array[0]
for i in range(1, n):
    predicted_load[i] = (load_array[0] - u / a) * (1 - math.exp(a)) * math.exp(-a * (i))

# 模型精度的后验差检验
residual_mean = 0
for i in range(0, n):
    residual_mean += (load_array[i] - predicted_load[i])
residual_mean /= n

# 求历史数据平均值
average_load = np.mean(load_array)

# 求历史数据方差
historical_variance = np.var(load_array)

# 求残差方差
residual_variance = 0
for i in range(0, n):
    residual_variance += ((load_array[i] - predicted_load[i]) - residual_mean) ** 2
residual_variance /= n

# 求后验差比值
C_ratio = residual_variance / historical_variance

# 求小误差概率
count = 0
for i in range(0, n):
    if abs((load_array[i] - predicted_load[i]) - residual_mean) < 0.6754 * math.sqrt(historical_variance):
        count += 1
P_probability = count / n

if C_ratio < 0.35 and P_probability > 0.95:
    # 预测精度为一级
    years_to_predict = 10
    print('往后%d年为：' % years_to_predict)
    future_load = np.zeros(years_to_predict)
    for i in range(0, years_to_predict):
        future_load[i] = (load_array[0] - u / a) * (1 - math.exp(a)) * math.exp(-a * (i + n))
    print(future_load)
else:
    print('灰色预测法不适用')