灰色关联度常用于分析影响因子与被影响因子的关联，是水论文的好东西
import numpy as np
def gray_correlation(refer, data, rho=0.5):
''' refer: 参照数列 (列向量)
data: 比较数列 (以列为单位)
rho: 分辨率
return: 灰色关联度'''
# 确保参照数列为列向量
refer = refer.reshape(-1, 1)
# 数列间的绝对距离
distance = np.abs(data - refer)
dis_min = distance.min()
dis_max = rho * distance.max()
# 关联系数: 比较数列中每个值与参照数列的关联性
cor_param = (dis_min + dis_max) / (distance + dis_max)
# 关联度: 关联系数按列求平均
degree = cor_param.mean(axis=0)
return degree
求解示例：
x0 = np.array([9, 9, 9, 9, 9, 9, 9])
xk = np.array([[8, 9, 8, 7, 5, 2, 9],
[7, 8, 7, 5, 7, 3, 8],
[9, 7, 9, 6, 6, 4, 7],
[6, 8, 8, 8, 4, 3, 6],
[8, 6, 6, 9, 8, 3, 8]])
# 比较数列 xk 应以列为单位, 故转置
degree = gray_correlation(x0, xk.T)
print(degree)
# OUTPUT: [0.71313131 0.61424774 0.68020215 0.5986346 0.68266821]
灰色预测模型
import logging
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
warnings.filterwarnings('ignore')
logging.basicConfig(format='%(message)s', level=logging.INFO)
LOGGER = logging.getLogger(__name__)
三种灰色数列的计算
# 灰积分: 累加数列
integrade = np.cumsum(seq)
# 灰微分: 差分数列
diff = np.diff(seq)
def gray_poly(seq, alpha=0.5):
''' return: 灰多项式 (加权邻值生成数列)'''
part_1 = alpha * seq[1:]
part_2 = (1 - alpha) * seq[:-1]
return part_1 + part_2
模型建立步骤：
- 记时间序列 x 的长度为 n，计算级比序列：
- 可容覆盖区间：，若 均在该区间内，则级比检验成功
- 灰导数 diff：时间序列 -> 累加数列 -> 差分数列
- 白化背景值 white：时间序列 -> 累加数列 -> 灰多项式
- 求解灰微分方程：diff + a · white = b，得到 a 和 b
- 白化模型预测函数：
- 相对残差：，均小于 0.1 则认为达到较高要求，均小于 0.2 则认为达到一般要求
- 级比偏差：，均小于 0.1 则认为达到较高要求，均小于 0.2 则认为达到一般要求
灰色预测模型一步到位代码如下：
def gray_model(seq, alpha=0.5, decimals=4, show=False):
''' 灰色预测模型建立
seq: 原始序列
alpha: 灰多项式权值
decimals: 数值精度
show: 绘制真实值和预测值的对比图
return: 灰色预测模型, 模型信息'''
seq = seq.reshape(-1)
length = len(seq)
index = np.arange(1, length + 1, 1)
# 创建数据存储表单
model_info = pd.DataFrame(index=index, columns=['真实值', '模型值', '相对误差', '级比偏差'])
model_info['真实值'] = seq
# 级比: x(k-1) / x(k)
mag_ratio = seq[:-1] / seq[1:]
mr_right = round(mag_ratio.max(), decimals)
mr_left = round(mag_ratio.min(), decimals)
# 级比检验的区间边界
left = round(np.exp(- 2 / (length + 1)), decimals)
right = round(np.exp(2 / (length + 1)), decimals)
# 级比检验: 检验数列级比是否都落在可容覆盖区间
LOGGER.info('n'.join(['级比检验:', f'MR ∈ [{mr_left}, {mr_right}]',
f'Border ∈ [{left}, {right}]', '']))
assert mr_left >= left and mr_right <= right, '序列未通过级比检验, 可通过平移变换改善'
# 灰积分: 累加数列
integrade = np.cumsum(seq)
# 灰微分: 差分数列
diff = np.diff(integrade)
# 灰多项式: 加权邻值生成数列
white = alpha * integrade[1:] + (1 - alpha) * integrade[:-1]
# 求解灰微分方程: diff + a * white = b
hidden = np.stack([-white, np.ones(white.shape)], axis=0)
# shape: [2, n] × [n, 2] × [2, n] × [n, ] -> [2, ]
a, b = (np.linalg.inv(hidden @ hidden.T) @ hidden @ diff)
LOGGER.info('n'.join([f'发展灰度: {a}', f'内生控制灰度: {b}']))
model_info['模型值'][1] = seq[0]
c2 = b / a
c1 = seq[0] - c2
def inte_pred(time):
''' 白化模型预测函数
time: 序列首个值对应的时间为 1'''
return c1 * (np.exp(- a * (time - 1)) - np.exp(- a * (time - 2)))
LOGGER.info(
f'预测方程: x(t) = {round(c1, decimals)}·[e^({-round(a, decimals)}(t-1)) - e^({-round(a, decimals)}(t-2))]n')
# 得到白化模型预测值
origin = seq[1:]
prediction = inte_pred(index[1:])
model_info['模型值'][1:] = np.round(prediction, decimals)
def check(error, label):
model_info[label][1:] = error
# 相对残差 / 级比偏差: < 0.1 达到较高要求, < 0.2 达到一般要求
if np.all(error < 0.1):
model_info[label][1] = 'Excelent'
elif np.all(error < 0.2):
model_info[label][1] = 'Common'
else:
model_info[label][1] = False
# 计算得到误差信息
rela_error = np.round(np.abs((origin - prediction) / origin), decimals)
check(rela_error, '相对误差')
radio_error = np.round(np.abs(1 - (1 - 0.5 * a) / (1 + 0.5 * a) * mag_ratio), decimals)
check(radio_error, '级比偏差')
LOGGER.info(model_info)
if show:
# 绘制真实值和预测值的对比图
plt.plot(index, model_info['真实值'], c='deepskyblue', label='true')
plt.scatter(index, model_info['模型值'], c='orange', label='pred', marker='p')
plt.legend()
plt.show()
return inte_pred, model_info
求解示例：
data = np.array([71.1, 72.4, 72.4, 72.1, 71.4, 72.0, 71.6])
# 只给模型提供前 5 个值, 利用后 2 个值验证模型精度
pred, model = gray_model(data[:5], decimals=4, alpha=0.5, show=True)
print(f'真实值: {data[-2]}, 模型值: {pred(6)}')
print(f'真实值: {data[-1]}, 模型值: {pred(7)}')
