灰色预测

数学建模
数学建模
发布日期:   2023-09-01
更新日期:   2024-04-15
文章字数:   4.8k
阅读时长:   24 分
阅读次数:  

灰色预测

灰色预测模型（Grey Prediction Model）是一种用于处理少量数据、难以建立准确数学模型的情况下进行预测的方法。这种模型常用于时间序列数据的预测，特别是在数据量较小、趋势难以准确刻画的情况下。

灰色预测模型的核心思想是通过对数据序列进行处理，将原始的不规则数据序列转化为灰色数据（累加生成数据）序列，然后基于灰色数据序列进行近似指数规律的预测。该方法的主要步骤包括建立灰色微分方程、求解参数、检验模型拟合度和进行预测等。

  • 优点是不需要很多的数据，一般只需要4个数据就够，能解决历史数据少、序列的完整性及可靠性低的问题；能利用微分方程来充分挖掘系统的本质，精度高；能将无规律的原始数据进行生成得到规律性较强的生成序列，运算简便，易于检验，具有不考虑分布规律，不考虑变化趋势。
  • 缺点是只适用于中短期的预测，只适合指数增长的预测。

（累加AGO与累减IAGO）在处理一些实际问题中，往往会遇到随机干扰，导致一些数据具有很大的波动性，为处理这些问题提出了累加和累减的概念。设原始数据序列 x(0)=[x1(0),x2(0),⋯ ,xn(0)]\mathbf{x}^{(0)}=[x^{(0)}_1,x^{(0)}_2,\cdots,x^{(0)}_n] ，则一次累加生成序列为x(1)=[x1(1),x2(1),⋯ ,xn(1)]\mathbf{x}^{(1)}=[x^{(1)}_1,x^{(1)}_2,\cdots,x^{(1)}_n]，其中 xk(1)=∑i=1kxi(0)x^{(1)}_k = \sum_{i=1}^{k} x^{(0)}_i。相应地，也有rr次累加生成序列。

（级数比）若对于所有 k=2,3,⋯ ,nk=2,3,\cdots,n，级比 λk=xk−1xk\lambda_k=\frac{x_{k-1}}{x_k} 落在区间 (e−2n+1,e2n+1)(e^{-\frac{2}{n+1}},e^{\frac{2}{n+1}}) 内，则称数据列满足指数形式增长。

累加的主要目的是把非负的波动数列转化成具有一定规律性（例如，指数形式单调增加）的数列。如果实际问题中出现负数（如温度数列），累加生成就不一定是好的处理办法，因为会出现正负抵消现象，这个时候会削弱原始数据的规律性。因此，此时应首先化为非负数列。具体做法是数列中每个数据同时减去原始数列中最小的元素值，得到非负数列后再进行累加运算。当然，在进行误差计算或预测时，应进行相应的逆运算。

（关联系数）选取参考数列 x(0)=[x1(0),x2(0),⋯ ,xn(0)]\mathbf{x}^{(0)}=[x^{(0)}_1,x^{(0)}_2,\cdots,x^{(0)}_n] 。假设有 mm个比较数列 y(i)y^{(i)}，则称

ξi(k)=min⁡smint∣xt(0)−yt(s)∣+ρmax⁡smaxt∣xt(0)−yt(s)∣∣xt(0)−yk(i)∣+ρmax⁡smaxt∣xt(0)−yt(s)∣\xi _{i}\left ( k \right ) = \frac{\min_{s}min_{t} \left | x^{(0)}_t-y^{(s)}_t \right |+\rho\max_{s}max_{t} \left | x^{(0)}_t - y^{(s)}_t \right | }{\left | x^{(0)}_t -y^{(i)}_k \right | + \rho \max_{s}max_{t} \left | x^{(0)}_t-y^{(s)}_t \right | }

为比较数列 y(i)\mathbf{y}^{(i)} 对参考数列 x(0)\mathbf{x}^{(0)} 在 kk 时刻的关联系数， 其中 ρ∈[0, 1]\rho \in \left [ 0,\: 1 \right ] 为分辨系数， min⁡smint∣xt(0)−yt(s)∣ ;  max⁡smaxt∣xt(0)−yt(s)∣\min_{s}min_{t} \left | x^{(0)}_t-y^{(s)}_t \right | \, ; \: \, \max_{s}max_{t} \left | x^{(0)}_t-y^{(s)}_t \right | 分别为两级最小差、两级最大差。一般来讲，分辨系数 ρ\rho 越大，分辨率越大； ρ\rho 越小，分辨率越小。

（关联度）称 ri=1n∑k=1nξi(k)r_i=\frac{1}{n}\sum_{k=1}^{n}\xi _{i}\left ( k \right ) 为序列 y(i)y^{(i)} 对 x(0)x^{(0)} 的关联度。

import numpy as np

def calculate_association(x0, xi, rho):
    min_diff = np.min(np.abs(np.subtract.outer(x0, xi)))
    max_diff = np.max(np.abs(np.subtract.outer(x0, xi)))
    
    association = (min_diff + rho * max_diff) / (np.abs(x0 - xi) + rho * max_diff)
    
    return association

# Example usage
x0 = np.array([1, 2, 3, 4, 5])
xi = np.array([2, 3, 4, 5, 7])
rho = 0.5

result = calculate_association(x0, xi, rho)
print(result)
print(sum(result))

GM(1,1)预测模型

数据序列在累加后呈现出指数形式的单调递增规律，联想到微分方程 y′=ayy'=ay 具有指数形式的解 y=eaxy=e^{ax}，由此提出一阶灰色方程模型GM(1,1)模型，其中的第一个1表示1阶微分方程，第二个1表示只含1个变量的灰色模型。

  1. 对数据序列 x(0)\mathbf{x}^{(0)} 进行累加生成

    x(1)=[x1(1),x2(1),⋯ ,xn(1)]=[x1(0),x1(0)+x2(0),⋯ ,∑i=1nxi(0)],\mathbf{x}^{(1)}=[x^{(1)}_1,x^{(1)}_2,\cdots,x^{(1)}_n]=[x^{(0)}_1,x^{(0)}_1+x^{(0)}_2,\cdots,\sum_{i=1}^{n} x^{(0)}_i],

  2. 对 x(1)\mathbf{x}^{(1)} 做均值生成

    z(1)=[z(2),z(3),⋯ ,z(n)],\mathbf{z}^{(1)}=[z^{(2)},z^{(3)},\cdots,z^{(n)}],

    其中 zk(1)=0.5xk(1)+0.5xk−1(1),k=2,3,⋯ ,nz^{(1)}_k=0.5x^{(1)}_k+0.5x^{(1)}_{k-1},k=2,3,\cdots,n。

  3. 数据的检验与处理

    为了保证数据建模的可行性，需要对已知数据序列作必要的检验处理。计算参考序列的级比。如果所有的级比都落在区间范围内，则该参考序列可以作为GM(1,1)的数据进行灰色预测。否则，继续对参考序列做必要的变化处理，使其落在区间内容。例如，选取适当的正数，做平移处理

    yk(0)=xk(0)+C,k=1,2,⋯ ,n.y^{(0)}_k=x^{(0)}_k+C,\quad k=1,2,\cdots,n.

    使得序列y\mathbf{y}满足级比。

  4. 建立微分方程模型

    dxt(1)dt+axt(1)=b\frac{dx^{(1)}_t}{dt}+ax^{(1)}_t=b

    为了计算系数 a,ba,b，在区间 k−1<t<kk-1<t<k 上，令 xt(1)=zk(1)x^{(1)}_t = z^{(1)}_k，构造 n−1n-1 组数据，

    dxt(1)dt=xk(1)−xk−1(1)=xk(0)\frac{dx^{(1)}_t}{dt}=x^{(1)}_k-x^{(1)}_{k-1}=x^{(0)}_k

    模型变为离散形式，

    xk(0)+azk(1)=bx^{(0)}_k+az^{(1)}_k=b

  5. 利用最小二乘法求解 a^,b^\hat{a},\hat{b}

    [a^,b^]T=(ATA)−1ATY[\hat{a},\hat{b}]^T=(\mathbf{A}^T\mathbf{A})^{-1}\mathbf{A}^T\mathbf{Y}

    其中 Y=[x2(0),x3(0),⋯ ,xn(0)]\mathbf{Y}=[x^{(0)}_2,x^{(0)}_3,\cdots,x^{(0)}_n]。

  6. 求解微分方程

    x^k+1(1)=(x1(0)−b^a^)e−a^k+b^a^,k=0,1,2,⋯\hat{x}^{(1)}_{k+1} =\left(x^{(0)}_{1}-\frac{\hat{b}}{\hat{a}}\right)e^{-\hat{a}k}+\frac{\hat{b}}{\hat{a}},\quad k=0,1,2,\cdots

    且 x^1(0)=x^1(1),x^k+1(0)=x^k+1(1)−x^k(1),k=1,2,⋯\hat{x}^{(0)}_1=\hat{x}^{(1)}_1,\hat{x}^{(0)}_{k+1}=\hat{x}^{(1)}_{k+1}-\hat{x}^{(1)}_{k},k=1,2,\cdots。

  7. 误差检验

    • 绝对误差分析（残差）
    • 相对误差分析（一般小于0.2认为达到一般要求，小于0.1认为达到较高的要求）
    • 级比偏差值分析
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sympy as sp

def GM_11(x0):
    # 灰色预测模型
    x1 = np.cumsum(x0)  # 累加生成序列
    z1 = (x1[:-1] + x1[1:]) / 2.0  # 紧邻均值生成序列
    B = np.vstack([-z1, np.ones_like(z1)]).T
    Yn = x0[1:].reshape((-1, 1))
    a, b = np.linalg.lstsq(B, Yn, rcond=None)[0]
    
    t = sp.symbols('t')
    x = sp.Function('x')
    eq = x(t).diff(t)+a*x(t)-b
    xt=sp.dsolve(eq,ics={x(0):x0[0]})
    print(xt[0].rhs)    
    xt = sp.lambdify(t, xt[0].rhs, 'numpy')

    t = np.arange(0, len(x0))
    x_predict = xt(t)
    
    #x_predict = (x0[0] - b / a) * np.exp(-a * np.arange(len(x0))) + b / a
    return x_predict

# 示例数据
data = np.array([25723, 30379, 34473, 38485, 40514, 42400, 48337])

n=len(data)
jibi = data[1:] / data[:-1]  # 计算级比
bd = [jibi.min(), jibi.max()]  # 级比范围
bd2 = [np.exp(-2/(n+1)),np.exp(2/(n+1))] # 级比范围

#判断 bd是否在bd2区间内
if bd[0] >= bd2[0] and bd[1] <= bd2[1]:
    # 灰色预测
    predicted_data = GM_11(data)
    
    predicted_data[1:] = predicted_data[1:]-predicted_data[:-1]
    
    print("最大相对误差：%d%%" % np.max(np.abs(data-predicted_data)/data*100))
    print("平均相对误差：%d%%" % np.mean(np.abs(data-predicted_data)/data*100))
    # 绘制原始数据和预测结果
    plt.plot(data, marker='o', label='Original Data')
    plt.plot(predicted_data, marker='s', linestyle='dashed', label='Predicted Data')
    plt.title('Grey Prediction Model')
    plt.xlabel('Time')
    plt.ylabel('Value')
    plt.legend()
    plt.show()
else:
    print('级比判别法判断级比序列为不稳定序列，不适用灰色预测模型进行预测')

GM(1,N)预测模型

设系统有 NN 个指标变量，对应的参考序列分别为 xi(0)=[xi1(0),xi2(0),⋯ ,xin(0)],i=1,2,⋯ ,Nx^{(0)}_i=[x^{(0)}_{i1},x^{(0)}_{i2},\cdots,x^{(0)}_{in}],i=1,2,\cdots,N，依次做累加生成。

  1. 建立微分方程

    dxit(1)dt=a1x1t(1)+a2x2t(1)+⋯++aNxNt(1)\frac{dx^{(1)}_{it}}{dt}=a_1x^{(1)}_{1t}+a_2x^{(1)}_{2t}+\cdots++a_Nx^{(1)}_{Nt}

  2. 化为离散形式

    xit(0)=a1x1t(1)+a2x2t(1)+⋯++aNxNt(1)x^{(0)}_{it}=a_1x^{(1)}_{1t}+a_2x^{(1)}_{2t}+\cdots++a_Nx^{(1)}_{Nt}

    其中系数 aka_k 是否为0与指标变量之间的相关性有关。

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

x10 = np.array([7.123, 7.994, 8.272, 7.960, 7.147, 7.092, 6.858,5.804,6.433,6.354,6.254,5.197,5.654])
x20 = np.array([0.796, 0.832, 0.917, 0.976, 1.075, 1.121,1.281,1.350,1.410,1.432,1.507,1.559,1.611])
x30 = np.array([13.108,12.334,12.216,12.201,12.132,11.990,11.926, 10.478,9.176,11.368,12.764,11.143,10.737])
x40 = np.array([27.475,27.473,27.490,27.500,27.510,27.542,27.536,27.550,27.567,27.584,27.600,27.602,27.630])

n=len(x10)
x11 =np.cumsum(x10) # 累加生成序列
x21 =np.cumsum(x20) # 累加生成序列
x31 =np.cumsum(x30) # 累加生成序列
x41 =np.cumsum(x40) # 累加生成序列

z1 = (x11[:-1] + x11[1:]) / 2.0  # 紧邻均值生成序列
z2 = (x21[:-1] + x21[1:]) / 2.0  # 紧邻均值生成序列
z3 = (x31[:-1] + x31[1:]) / 2.0  # 紧邻均值生成序列
z4 = (x41[:-1] + x41[1:]) / 2.0  # 紧邻均值生成序列

B1=np.c_[z1,np.ones((n-1,1))]
u1=np.linalg.pinv(B1).dot(x10[1:]); print(u1)
B2=np.c_[z1,z2]    # x2与x1有关
u2=np.linalg.pinv(B2).dot(x20[1:]); print(u2)
B3=np.c_[z3,np.ones((n-1,1))] 
u3=np.linalg.pinv(B3).dot(x30[1:]); print(u3)
B4=np.c_[z1,z3,z4]  # x4是x1,x2,x3的综合反应
u4=np.linalg.pinv(B4).dot(x40[1:]); print(u4)

def Pfun(x,t):
    x1, x2, x3, x4 = x
    return np.array([u1[0]*x1+u1[1], u2[0]*x1+u2[1]*x2, u3[0]*x3+u3[1], u4[0]*x1+u4[1]*x3+u4[2]*x4])

t=np.arange(0, n)

X0=np.array([x10[0],x20[0],x30[0],x40[0]])
s1=odeint(Pfun, X0, t)
s2=np.diff(s1,axis=0)
xh=np.vstack([X0,s2])

# Plot the original data
plt.figure(figsize=(10, 6))
plt.plot(t, x10, label='x1')
plt.plot(t, x20, label='x2')
plt.plot(t, x30, label='x3')
plt.plot(t, x40, label='x4')

# Plot the simulated data
plt.plot(t, xh[:, 0], 'o--', label='Simulated x1')
plt.plot(t, xh[:, 1], 'o--', label='Simulated x2')
plt.plot(t, xh[:, 2], 'o--', label='Simulated x3')
plt.plot(t, xh[:, 3], 'o--', label='Simulated x4')

plt.title('Simulation Results')
plt.xlabel('Time')
plt.ylabel('Values')
plt.legend()
plt.show()

注：GM(1,1)模型适用于具有较强指数规律的序列，只能描述单调的变化过程，对于非单调的摆动发展序列或有饱和的S形序列，可以考虑建立GM(2,1)，DGM和Verhulst等模型。

GM(2,1)

适用于非单调摆动发展序列。更多地用于处理不确定性和变动性较大的系统，进行短期预测。对原始序列做一次累加和一次累减。

  1. 建立微分方程

    d2xt(1)d2t+a1dxt(1)dt+a2xt(1)=b\frac{d^2x^{(1)}_t}{d^2t}+a_1\frac{dx^{(1)}_t}{dt}+a_2x^{(1)}_t = b

  2. 转化为离散形式

    xk(−1)+a1xk(0)+a2zk(1)=bx^{(-1)}_k+a_1x^{(0)}_k+a_2z^{(1)}_k=b 

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sympy as sp


def GM_11(x0):
    # 灰色预测模型
    x1 = np.cumsum(x0)  # 累加生成序列
    z1 = (x1[:-1] + x1[1:]) / 2.0  # 紧邻均值生成序列
    B = np.vstack([-x0[1:], -z1, np.ones_like(z1)]).T
    Yn = np.diff(x0)
    a1, a2, b = np.linalg.lstsq(B, Yn, rcond=None)[0]
    print(a1, a2, b)
    t = sp.symbols('t')
    x = sp.Function('x')
    eq = x(t).diff(t, 2) + a1 * x(t).diff(t) + a2 * x(t) - b
    xt = sp.dsolve(eq, ics={x(0): x0[0], x(5): x1[-1]})
    print(xt.rhs)
    xt = sp.lambdify(t, xt.rhs, 'numpy')

    t = np.arange(0, len(x0))
    x_predict = xt(t)

    # x_predict = (x0[0] - b / a) * np.exp(-a * np.arange(len(x0))) + b / a
    return x_predict

# 示例数据
data = np.array([41, 49, 61, 78, 96, 104])

n = len(data)
jibi = data[1:] / data[:-1]  # 计算级比
bd = [jibi.min(), jibi.max()]  # 级比范围
bd2 = [np.exp(-2 / (n + 1)), np.exp(2 / (n + 1))]  # 级比范围

# 判断 bd是否在bd2区间内
if bd[0] >= bd2[0] and bd[1] <= bd2[1]:
    # 灰色预测
    predicted_data = GM_11(data)

    predicted_data[1:] = predicted_data[1:] - predicted_data[:-1]

    print("最大相对误差：%d%%" % np.max(np.abs(data - predicted_data) / data * 100))
    print("平均相对误差：%d%%" % np.mean(np.abs(data - predicted_data) / data * 100))
    # 绘制原始数据和预测结果
    plt.plot(data, marker='o', label='Original Data')
    plt.plot(predicted_data, marker='s', linestyle='dashed', label='Predicted Data')
    plt.title('Grey Prediction Model')
    plt.xlabel('Time')
    plt.ylabel('Value')
    plt.legend()
    plt.show()
else:
    print('级比判别法判断级比序列为不稳定序列，不适用灰色预测模型进行预测')

DGM(2,1)模型

适用于非单调摆动发展序列。对原始序列做一次累加和一次累减。

  1. 建立微分方程

    d2xt(1)d2t+adxt(1)dt=b\frac{d^2x^{(1)}_t}{d^2t}+a\frac{dx^{(1)}_t}{dt}= b

  2. 转化为离散形式

    xk(−1)+axk(0)=bx^{(-1)}_k+ax^{(0)}_k=b 

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sympy as sp


def GM_11(x0):
    # 灰色预测模型
    x1 = np.cumsum(x0)  # 累加生成序列
    z1 = (x1[:-1] + x1[1:]) / 2.0  # 紧邻均值生成序列
    B = np.vstack([-x0[1:], np.ones_like(z1)]).T
    Yn = np.diff(x0)
    a, b = np.linalg.lstsq(B, Yn, rcond=None)[0]
    print(a, b)
    t = sp.symbols('t')
    x = sp.Function('x')
    eq = x(t).diff(t, 2) + a * x(t).diff(t) - b
    xt = sp.dsolve(eq, ics={x(0): x0[0], x(5): x1[-1]})
    print(xt.rhs)
    xt = sp.lambdify(t, xt.rhs, 'numpy')

    t = np.arange(0, len(x0))
    x_predict = xt(t)

    # x_predict = (x0[0] - b / a) * np.exp(-a * np.arange(len(x0))) + b / a
    return x_predict

# 示例数据
data = np.array([2.874,3.278,3.39,3.679,3.77,3.8])

n = len(data)
jibi = data[1:] / data[:-1]  # 计算级比
bd = [jibi.min(), jibi.max()]  # 级比范围
bd2 = [np.exp(-2 / (n + 1)), np.exp(2 / (n + 1))]  # 级比范围

# 判断 bd是否在bd2区间内
if bd[0] >= bd2[0] and bd[1] <= bd2[1]:
    # 灰色预测
    predicted_data = GM_11(data)

    predicted_data[1:] = predicted_data[1:] - predicted_data[:-1]

    print("最大相对误差：%d%%" % np.max(np.abs(data - predicted_data) / data * 100))
    print("平均相对误差：%d%%" % np.mean(np.abs(data - predicted_data) / data * 100))
    # 绘制原始数据和预测结果
    plt.plot(data, marker='o', label='Original Data')
    plt.plot(predicted_data, marker='s', linestyle='dashed', label='Predicted Data')
    plt.title('Grey Prediction Model')
    plt.xlabel('Time')
    plt.ylabel('Value')
    plt.legend()
    plt.show()
else:
    print('级比判别法判断级比序列为不稳定序列，不适用灰色预测模型进行预测')

Verhulst模型

该模型主要用来描述具有饱和状态的过程，即 S 型过程，常用于人口预测，生物生长，繁殖预测及产品经济寿命预测等。

  1. 建立微分方程

    dxt(1)dt+axt(1)=b(xt(1))2\frac{dx^{(1)}_t}{dt}+ax^{(1)}_t=b\left(x^{(1)}_t\right)^2

  2. 转为离散形式

    xk(0)+azk(1)=b(zk(1))2x^{(0)}_k+az^{(1)}_k=b\left(z^{(1)}_k\right)^2

灰色波形预测

当原始数据频繁波动且摆动幅度较大时，可以考虑根据原始数据的波形预测未来行为数据发展变化的波形。

  1. 绘制原始序列的 kk 段折线图

    {yK=xk+(t−k)[xk+1−xk]}\{y_K = x_k+(t-k)[x_{k+1}-x_k]\}

  2. 找等高线

找等高线极点 σmax=max1≤i≤nxi,σmin=min1≤i≤nxi\sigma_{max}=max_{1 \leq i \leq n} x_i,\sigma_{min}=min_{1 \leq i \leq n} x_i。∀β∈[σmin,σmax]\forall \beta \in [\sigma_{min},\sigma_{max}] 建立等高线方程求等高点。一般等距离平均分成 pp 段。

  1. 对等高点序列（这里指的是日期序列，而不是数值序列）建立GM(1,1)模型并求解
  2. 剔除无效点

观察上一步中得到的等高预测序列中是否存在相同值，若存在即意味着该时刻对应两个不同预测值，该点无效，应当剔除。如果不存在无效点，因而只需对等高预测序列得到的时刻按从小到大排序，最后将每个时刻与其所在序列的等高线数值对应，作出预测波形。

灰色波形预测的缺点是，只能在设定的点处变化，不会超出最大最小区间。

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt

data = np.array(
    [2276.67, 2150.76, 1929.05, 2216.81, 2092.22, 1994.67, 1895.82, 1719.81, 1760.61, 1859.11, 2017.47, 1897.88,
     1965.41, 2079.12, 1976.82, 1863.8, 1820.81, 1924.01, 1928.87, 1985.02, 2107.75, 2260.82, 2209.86, 2206.57, 2198.11,
     2139.03])

Q_dict = {}

step = 9
thresholds = np.linspace(min(data), max(data), step)

for threshold in thresholds:
    q_values = []
    for i in range(0, len(data) - 1):
        if (threshold >= data[i] and threshold <= data[i + 1]) or (threshold <= data[i] and threshold >= data[i + 1]):
            q = i + 1 + (threshold - data[i]) / (data[i + 1] - data[i])
            q_values.append(q)
    if len(q_values) >= 4:
        Q_dict[threshold] = np.array(q_values)


def GM_11(x0, stop_threshold_diff=50):
    # 灰色预测模型
    x1 = np.cumsum(x0)  # 累加生成序列
    z1 = (x1[:-1] + x1[1:]) / 2.0  # 紧邻均值生成序列
    B = np.vstack([-z1, np.ones_like(z1)]).T
    Yn = x0[1:].reshape((-1, 1))
    a, b = np.linalg.lstsq(B, Yn, rcond=None)[0]

    t = sp.symbols('t')
    x = sp.Function('x')
    eq = x(t).diff(t) + a * x(t) - b
    xt = sp.dsolve(eq, ics={x(0): x0[0]})
    print(xt[0].rhs)
    xt = sp.lambdify(t, xt[0].rhs, 'numpy')

    t_val = len(x0) - 1
    x_predict = np.array([])

    while True:
        x_val = xt(t_val)

        if len(x_predict) > 0 and x_val > x_predict[-1] + stop_threshold_diff:
            break

        x_predict = np.append(x_predict, x_val)
        t_val += 1

    return x_predict


# 初始化拼接数组
all_x_values = np.array([])
all_y_values = np.array([])

for threshold, q_values in Q_dict.items():
    # 灰色预测
    predicted_data = GM_11(q_values, stop_threshold_diff=50)

    # 计算差分
    predicted_data = predicted_data[1:] - predicted_data[:-1]

    # 将预测数据添加到拼接数组中
    all_x_values = np.append(all_x_values, predicted_data)
    all_y_values = np.append(all_y_values, np.full_like(predicted_data, threshold))

# 对数组按照 x 值进行排序
selected_indices = all_x_values > len(data)
all_x_values = all_x_values[selected_indices]
all_y_values = all_y_values[selected_indices]

sorted_indices = np.argsort(all_x_values)
all_x_values = all_x_values[sorted_indices]
all_y_values = all_y_values[sorted_indices]

all_x_values = np.append(len(data), all_x_values)
all_y_values = np.append(data[-1], all_y_values)

# 绘制折线图
plt.plot(data, marker='o', label='Original Data')
plt.plot(all_x_values, all_y_values, marker='s', label='Predicted Data')
plt.title('Grey Prediction Model')
plt.xlabel('Time')
plt.ylabel('Values')
plt.legend()
plt.show()
文章作者: JiJunhao
文章链接: http://jijunhao.github.io/2023/09/01/article20230901/
版权声明: 本博客所有文章除特別声明外，均采用 CC BY 4.0 许可协议。转载请注明来源 JiJunhao !
数学建模
赏
  目录