Python 灰色关联度 & 灰色预测模型

灰色关联度

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

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

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

模型建立步骤：

1. 记时间序列 x 的长度为 n，计算级比序列：
2. 可容覆盖区间：，若  均在该区间内，则级比检验成功
3. 灰导数 diff：时间序列 -> 累加数列 -> 差分数列
4. 白化背景值 white：时间序列 -> 累加数列 -> 灰多项式
5. 求解灰微分方程：diff + a · white = b，得到 a 和 b
6. 白化模型预测函数：
7. 相对残差：，均小于 0.1 则认为达到较高要求，均小于 0.2 则认为达到一般要求
8. 级比偏差：，均小于 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}']))
    c2 = b / a
    c1 = seq[0] - c2

    def 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

求解示例：