本篇探究如何把灰色预测模型（主要是GM(1, N)）用好。
灰色预测模型（Grey Prediction Model）是一种用于处理小样本和不确定性数据的预测方法。其核心思想是利用数据序列的内部规律，通过生成灰色数列来进行预测。GM(1,N)是其中的多变量模型，可以处理多个因素影响下的预测问题。
基本GM(1, N)模型
假如现在已有DataFrame格式原始数据df_org，其columns中的第一个为目标变量，其余的为解释变量。
Python
import numpy as np
import pandas as pd
from numpy.linalg import inv
from sklearn.preprocessing import StandardScaler
def cumulative_series(series):
return np.cumsum(series)
def predict_future(x0, a, b, k):
return (x0[0] - sum(b) / a) * np.exp(-a * k) + sum(b) / a
# 0 标准化
scaler = StandardScaler()
df = pd.DataFrame(scaler.fit_transform(df_org), columns=df_org.columns)
# 1 累加生成序列
df_cum = df.apply(cumulative_series)
Y = df_cum.iloc[1:, 0].values.reshape(-1, 1)
# 2 生成背景值（均值生成序列）
background_value = -0.5 * (df_cum.iloc[1:, 0].values + df_cum.iloc[:-1, 0].values)
B = np.column_stack([background_value, df_org.iloc[1:, 1:]])
# 3 求解模型参数
B_inv = inv(B.T @ B) @ B.T
coeff = B_inv @ Y
print(coeff.flatten()) # 如果要输出参数
# 4 建立预测模型
a = coeff[0, 0] # 参数a
b = coeff[1:, 0] # 参数b
cnt_list = []
for year in range(5): # 比如预测未来5年
future_cnt = predict_future(df_org.iloc[0].values, a, b, year + 1)
transform_item = [future_cnt] + [0 for _ in range(len(coeff) - 1)]
future_cnt_org = scaler.inverse_transform([transform_item])[0, 0]
cnt_list.append(future_cnt_org)
print(f'未来第{year + 1}年预测值:', future_cnt_org) # 如果要输出结果
可能的问题？
指数爆炸
原始模型是难以规避指数爆炸问题的。为了解决这一问题，可选方式包括：
- 添加基础值：同列解释变量增加一个基础值，使不同值之间的差异更加缓和。
- 添加衰减量：对进入predict_future()的代表年份的k值进行衰减，如增加一个负指数或对数系数项。
- 添加衰减拐点：使预测值从某个年份开始出现衰减或衰减量增大。
解释变量相关性
有时纳入模型的解释变量本身具有强相关性，此时可以考虑对其进行合成，即：对所有解释变量进行标准化 -> 同类解释变量加和合并 -> 再次进行标准化。
修改版GM(1, N)模型
沿用基础版本的df_org，假设其columns=[x, y1, y2, y3, y4, y5]，且其中y2和y3有强相关性，y4需要增加基础值。
Python
import numpy as np
import pandas as pd
from numpy.linalg import inv
from sklearn.preprocessing import StandardScaler
def cumulative_series(series):
return np.cumsum(series)
def predict_future(x0, a, b, k):
def get_k_prime(k): # 改动3 & 改动4的应用
if k + 1 <= inflection_key:
return k
else:
return k * np.exp(-alpha * (np.log(k + 1))) # 这里只是简单的例子，若要确保预测曲线平滑则可能需要更多的工作
k_prime = get_k_prime(k)
return (x0[0] - sum(b) / a) * np.exp(-a * k_prime) + sum(b) / a
# 0 标准化
df_org['y4'] += 100 # 【改动1】 解释变量添加基础值
scaler = StandardScaler()
df = pd.DataFrame(scaler.fit_transform(df_org), columns=df_org.columns)
df['y23'] = df['y2'] + df['y3'] # 【改动2】 同类解释变量合并
df_inter = df.drop(['y2', 'y3'], axis=1)
df = pd.DataFrame(scaler.fit_transform(df_inter), columns=df_inter.columns)
# 1 累加生成序列
df_cum = df.apply(cumulative_series)
Y = df_cum.iloc[1:, 0].values.reshape(-1, 1)
# 2 生成背景值（均值生成序列）
background_value = -0.5 * (df_cum.iloc[1:, 0].values + df_cum.iloc[:-1, 0].values)
B = np.column_stack([background_value, df_org.iloc[1:, 1:]])
# 3 求解模型参数
B_inv = inv(B.T @ B) @ B.T
coeff = B_inv @ Y
print(coeff.flatten()) # 如果要输出参数
# 4 建立预测模型
a = coeff[0, 0] # 参数a
b = coeff[1:, 0] # 参数b
alpha = 0.4 # 【改动3】 添加衰减量
inflection_key = 3 # 【改动4】 添加衰减拐点
cnt_list = []
for year in range(5): # 比如预测未来5年
future_cnt = predict_future(df_org.iloc[0].values, a, b, year + 1)
transform_item = [future_cnt] + [0 for _ in range(len(coeff) - 1)]
future_cnt_org = scaler.inverse_transform([transform_item])[0, 0]
cnt_list.append(future_cnt_org)
print(f'未来第{year + 1}年预测值:', future_cnt_org) # 如果要输出结果
