灰色关联度

灰色关联度常用于分析影响因子与被影响因子的关联，是水论文的好东西

如果数据的量纲不统一的话，需要先进行归一化处理

import numpy as npdef 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 warningsimport matplotlib.pyplot as plt
import numpy as np
import pandas as pdwarnings.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，计算级比序列： $m_t = \frac{x_{t+1}}{x_{t}}, t\in [1, n-1]$
可容覆盖区间： $(e^{-\frac{2}{n+1}}, e^{\frac{2}{n+1}})$ ，若 $m_t$ 均在该区间内，则级比检验成功
灰导数 diff：时间序列 -> 累加数列 -> 差分数列
白化背景值 white：时间序列 -> 累加数列 -> 灰多项式
求解灰微分方程：diff + a · white = b，得到 a 和 b
白化模型预测函数： $\hat{x_t}=(x_1-\frac{b}{a}) \cdot (e^{-a \cdot (t-1)}-e^{-a \cdot (t-2)}), t \in [1, \infty ]$
相对残差： $e_t=\frac{|\hat{x_t}-x_t|}{x_t}$ ，均小于 0.1 则认为达到较高要求，均小于 0.2 则认为达到一般要求
级比偏差： $e_k=|1-\frac{1-0.5a}{1+0.5a} m_t|$ ，均小于 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 = bhidden = 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}']))c2 = b / ac1 = seq[0] - c2def inte_pred(time):''' 白化模型预测函数time: 序列首个值对应的时间为 1'''assert np.all(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')# 得到白化模型预测值pred = inte_pred(index)model_info['模型值'] = np.round(pred, decimals)# 计算得到误差信息model_info['相对误差'] = np.round(np.abs((seq - pred) / seq), decimals)model_info['级比偏差'] = np.concatenate([np.zeros(1), np.round(np.abs(1 - (1 - 0.5 * a) / (1 + 0.5 * a) * mag_ratio), decimals)])# 相对残差 / 级比偏差: < 0.1 达到较高要求, < 0.2 达到一般要求for label in ['相对误差', '级比偏差']:error = model_info[label]if np.all(error < 0.1):message = 'Excelent'elif np.all(error < 0.2):message = 'Common'else:raise AssertionError(f'{label}: Fail {error.max()}')LOGGER.info(f'{label}: {message}')# 输出模型信息LOGGER.info('')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

求解示例：