数学建模算法：灰色预测模型GM(1,1)及Python代码

灰色预测模型GM(1,1)

灰色预测模型\(GM(1,1)\)是在数学建模比赛中常用的预测值方法，常用于中短期符合指数规律的预测。

其数学表达与原理分析参考文章尾部网页与文献资料。

  • 预处理

灰色模型要求数据前后级比落入范围 \(\displaystyle \Theta\left(e^{-\frac{2}{n+1}},e^{\frac{2}{n+2}}\right)\) ，因此做线性平移预处理使得元数据满足要求。

该预处理保证了数据一定落入级比范围内。

线性平移：将数据平移至不小于1，检查级比，若不满足要求则将数据向上平移一个最小值直到满足要求。

级比范围

可以推断出，级比的上下界在给定数据点数越多的情况下，越趋于1。

  • 代码

经过整理，以下附上Python代码：

import numpy as np
import matplotlib.pyplot as plt

# 线性平移预处理，确保数据级比在可容覆盖范围
def greyModelPreprocess(dataVec):
    "Set linear-bias c for dataVec"
    import numpy as np
    from scipy import io, integrate, linalg, signal
    from scipy.sparse.linalg import eigs
    from scipy.integrate import odeint

    c = 0
    x0 = np.array(dataVec, float)
    n = x0.shape[0]
    L = np.exp(-2/(n+1))
    R = np.exp(2/(n+2))
    xmax = x0.max()
    xmin = x0.min()
    if (xmin < 1):
        x0 += (1-xmin)
        c += (1-xmin)
    xmax = x0.max()
    xmin = x0.min()
    lambda_ = x0[0:-1] / x0[1:]  # 计算级比
    lambda_max = lambda_.max()
    lambda_min = lambda_.min()
    while (lambda_max > R or lambda_min < L):
        x0 += xmin
        c += xmin
        xmax = x0.max()
        xmin = x0.min()
        lambda_ = x0[0:-1] / x0[1:]
        lambda_max = lambda_.max()
        lambda_min = lambda_.min()
    return c

# 灰色预测模型
def greyModel(dataVec, predictLen):
    "Grey Model for exponential prediction"
    # dataVec = [1, 2, 3, 4, 5, 6]
    # predictLen = 5
    import numpy as np
    from scipy import io, integrate, linalg, signal
    from scipy.sparse.linalg import eigs
    from scipy.integrate import odeint

    x0 = np.array(dataVec, float)
    n = x0.shape[0]
    x1 = np.cumsum(x0)
    B = np.array([-0.5 * (x1[0:-1] + x1[1:]), np.ones(n-1)]).T
    Y = x0[1:]
    u = linalg.lstsq(B, Y)[0]

    def diffEqu(y, t, a, b):
        return np.array(-a * y + b)

    t = np.arange(n + predictLen)
    sol = odeint(diffEqu, x0[0], t, args=(u[0], u[1]))
    sol = sol.squeeze()
    res = np.hstack((x0[0], np.diff(sol)))
    return res

# 输入数据
x = np.array([-18, 0.34, 4.68, 8.49, 29.84, 50.21, 77.65, 109.36])
c = greyModelPreprocess(x)
x_hat = greyModel(x+c, 5)-c

# 画图
t1 = range(x.size)
t2 = range(x_hat.size)
plt.plot(t1, x, color='r', linestyle="-", marker='*', label='True')
plt.plot(t2, x_hat, color='b', linestyle="--", marker='.', label="Predict")
plt.legend(loc='upper right')
plt.xlabel('xlabel')
plt.ylabel('ylabel')
plt.title('Prediction by Grey Model (GM(1,1))')
plt.show()

预测值与真实值图像

灰色预测模型GM(2,1)

根据资料，\(GM(2,1)\)适用于非单调的摆动发展序列，或有饱和的S型序列。但是从图像上观察，数据预测由于为指数类型，变化过于夸张、预测趋势也有偏离的状况，可能实用性和普适性并不如\(GM(1,1)\)。

# 灰色预测模型GM(2,1)
def greyModel2(dataVec, predictLen):
    "Grey Model for exponential prediction"
    # dataVec = [1, 2, 3, 4, 5, 6]
    # predictLen = 5
    import numpy as np
    import sympy as sy
    from scipy import io, integrate, linalg, signal

    x0 = np.array(dataVec, float)
    n = x0.shape[0]
    a_x0 = np.diff(x0) # 1次差分序列
    x1 = np.cumsum(x0) # 1次累加序列
    z = 0.5 * (x1[0:-1] + x1[1:]) # 均值生成序列
    B = np.array([-x0[1:], -z, np.ones(n-1)]).T
    u = linalg.lstsq(B, a_x0)[0]

    def diffEqu(x, f, a1, a2, b):
        return sy.diff(f(x), x, 2) + a1*sy.diff(f(x), x) + a2*f(x) - b # f''(x)+a1*f'(x)+a2*f(x)=b 二阶常系数齐次微分方程

    t = np.arange(n + predictLen)
    x = sy.symbols('x')  # 约定变量
    f = sy.Function('f')  # 约定函数
    eq = sy.dsolve(diffEqu(x, f, u[0], u[1], u[2]), f(x), ics={f(t[0]): x1[0], f(t[n-1]): x1[-1]})
    f = sy.lambdify(x, eq.args[1], 'numpy')
    sol = f(t)
    res = np.hstack((x0[0], np.diff(sol)))
    return res

误差分析部分：可就绝对误差、相对误差、级比、残差做数据分析，以下示例为最小二乘法线性回归分析。

def regressionAnalysis(xVec, yVec):
    import numpy as np
    from scipy import linalg

    x = np.array([xVec, np.ones(xVec.size)]).T
    u = linalg.lstsq(x, yVec)
    return u

res = regressionAnalysis(x_hat[0:x.size], x)
print(res[0]) # 回归系数a, b
print(res[1]) # residuals 残差平方和

参考

posted @ 2022-03-03 20:11  Chiron-zy  阅读(10151)  评论(0)  编辑  收藏  举报