你已经了解了GM(1,1)模型,它是用来处理和预测只有一个因素(变量)随时间变化情况的模型,而且它是基于一阶微分方程的。这个模型假设数据之间的增长率是恒定的,但是在实际情况中,很多时候数据的增长率并不是一成不变的。
这时候,灰色模型中的GM(2,1)就派上用场了,GM(2,1)模型仍然是针对单一变量的,只不过它使用的是二阶微分方程。那个“2”表示的是模型预测时考虑了数据变化率的变化率,也就是数据增长速度的速度。这样就能让模型更加精细地捕捉数据随时间变化的趋势,尤其是当数据变化不是那么规则,或者有点起起伏伏时。
简单来说,如果你觉得GM(1,1)模型对某些数据的预测不够精确,或者数据的走势有较大的波动,你可能会想尝试用GM(2,1)模型,它通过考虑更复杂的变化趋势来提供可能更精确的预测结果。因此,GM(2,1)模型适用于那些数据变化比较复杂,一阶模型难以准确捕捉其变化趋势的情况。
以下是GM(2,1)模型的构建和求解过程。
原始数据序列:
首先,我们有一个原始数据序列,它是时间序列中的观察值
其中 表示第 个时刻的观测值。
一次累加生成操作(I-AGO)序列 :
通过对原始数据序列进行累加,我们得到一次累加序列,这有助于展现数据的增长趋势和结构:
得到
生成一次累加逆生成操作(I-IAGO)序列 ,是 I-AGO 序列的逆操作,它的目的是从累加序列中重构原始数据序列。基本上,逆生成操作是对累加序列 进行差分,以恢复原始的数据值。
其中,
生成 的紧邻均值生成序列 ,紧邻均值生成序列 是通过计算 序列中相邻两个数的平均值得到的,它的目的是用来构建模型方程:
建立GM(2,1)模型:
这个方程中的 , , 和 是我们需要通过数据来估计的参数。
建立微分方程模型:基于一次累加序列,我们构建一个包含未知参数的二阶微分方程,即GM(2,1)的白化方程。
其中,,,和 是模型参数,需要通过数据拟合来确定。
求解模型参数:为了估计模型参数,我们构建矩阵 和向量 :
参数通过最小二乘估计得到:
模型预测:有了模型参数后,我们可以通过解上述微分方程,使用这些参数对未来的数据点进行预测。
还原预测结果:将预测的累加序列 通过逆累加生成操作转换回原始数据序列 ,得到预测的原始数据值,以便于与实际观察值进行比较。
这个过程的核心是通过累加操作平滑数据,建立和求解微分方程来预测数据的趋势。GM(2,1)模型特别适合预测少量数据和具有一定发展趋势的时间序列。
1import numpy as np
2from sympy import symbols, Function, Eq, diff, dsolve
3from numpy.linalg import lstsq
4
5t = symbols('t') # 定义符号变量 t
6x = Function('x')(t) # 定义函数 x(t)
7
8# 数据初始化
9x0 = np.array([41,49,61,78,96,104]) # 输入原始序列
10n = len(x0) # 原始序列元素个数
11x1 = np.cumsum(x0) # 原始序列的累积和
12a_x0 = np.diff(x0) # 对原始序列进行一阶差分
13a_x0 = np.concatenate([[0], a_x0]) # 在差分结果前面补0,保持序列长度一致
14
15z = np.zeros(n) # 初始化z序列
16for i in range(1, n): # 对每一个元素:
17 z[i] = 0.5 * (x1[i] + x1[i - 1]) # 计算相邻元素的平均值
18
19# 构建数组 B 和 Y,用于之后的方程求解
20B = np.column_stack((-x0[1:], -z[1:], np.ones(n - 1))) # 构建 B,其中包含负的 x0 序列元素(从第二个元素开始)、负的 z 序列元素(从第二个元素开始)以及单位列向量
21Y = a_x0[1:] # 构建 Y,由一阶差分序列 a_x0 (从第二个元素开始)
22
23# 求解方程,获得参数 u
24u, _, _, _ = lstsq(B, Y, rcond=None) # 运用 Numpy 的最小二乘法 lstsq 求解线性方程得到 u
25
26# 定义微分方程
27a1, a2, b, = u[0], u[1], u[2] # 从参数 u 提取所需参数 a1,a2 和 b
28diff_eq = Eq(diff(x, t, t) + a1 * diff(x, t) + a2 * x, b) # 定义微分方程
29
30# 设定初始条件
31ics = {x.subs(t, 0): x1[0], x.subs(t, 5): x1[-1]} # 设定初始条件为 x(0) = x1[0],x(5) = x1[5] (注意这里索引是以 0 开始的)
32
33# 解微分方程
34solution = dsolve(diff_eq, ics=ics) # 解微分方程
35
36# 运用求解出的微分方程进行预测
37yuce = [solution.rhs.subs(t, val).evalf() for val in range(n)] # 对给出的数据序列进行预测
38x0_hat = [yuce[0]] + np.diff(yuce).tolist() # 对预测结果的序列计算一阶差分
39
40# 计算预测误差和相对误差
41epsilon = x0 - x0_hat # 计算预测值和实际值之间的差,即预测误差
42delta = np.abs(epsilon / x0) # 计算绝对比例误差,也就是相对误差
1import numpy as np
2from sympy import symbols, Function, Eq, diff, dsolve
3from numpy.linalg import lstsq
4
5t = symbols('t') # 定义符号变量 t
6x = Function('x')(t) # 定义函数 x(t)
7
8# 数据初始化
9x0 = np.array([41,49,61,78,96,104]) # 输入原始序列
10n = len(x0) # 原始序列元素个数
11x1 = np.cumsum(x0) # 原始序列的累积和
12a_x0 = np.diff(x0) # 对原始序列进行一阶差分
13a_x0 = np.concatenate([[0], a_x0]) # 在差分结果前面补0,保持序列长度一致
14
15z = np.zeros(n) # 初始化z序列
16for i in range(1, n): # 对每一个元素:
17 z[i] = 0.5 * (x1[i] + x1[i - 1]) # 计算相邻元素的平均值
18
19# 构建数组 B 和 Y,用于之后的方程求解
20B = np.column_stack((-x0[1:], -z[1:], np.ones(n - 1))) # 构建 B,其中包含负的 x0 序列元素(从第二个元素开始)、负的 z 序列元素(从第二个元素开始)以及单位列向量
21Y = a_x0[1:] # 构建 Y,由一阶差分序列 a_x0 (从第二个元素开始)
22
23# 求解方程,获得参数 u
24u, _, _, _ = lstsq(B, Y, rcond=None) # 运用 Numpy 的最小二乘法 lstsq 求解线性方程得到 u
25
26# 定义微分方程
27a1, a2, b, = u[0], u[1], u[2] # 从参数 u 提取所需参数 a1,a2 和 b
28diff_eq = Eq(diff(x, t, t) + a1 * diff(x, t) + a2 * x, b) # 定义微分方程
29
30# 设定初始条件
31ics = {x.subs(t, 0): x1[0], x.subs(t, 5): x1[-1]} # 设定初始条件为 x(0) = x1[0],x(5) = x1[5] (注意这里索引是以 0 开始的)
32
33# 解微分方程
34solution = dsolve(diff_eq, ics=ics) # 解微分方程
35
36# 运用求解出的微分方程进行预测
37yuce = [solution.rhs.subs(t, val).evalf() for val in range(n)] # 对给出的数据序列进行预测
38x0_hat = [yuce[0]] + np.diff(yuce).tolist() # 对预测结果的序列计算一阶差分
39
40# 计算预测误差和相对误差
41epsilon = x0 - x0_hat # 计算预测值和实际值之间的差,即预测误差
42delta = np.abs(epsilon / x0) # 计算绝对比例误差,也就是相对误差
这个代码提供了一个基于 GM(2,1) 模型对时间序列数据进行预测的完整工作流程,包含了建模、求解、预测和可视化的各个步骤:
以上仅为模型的参考代码,如果你在实践过程中预测的值相比于原始数据变化幅度过大,说明原始数据可能并不适合该模型。
在实践中,GM(2,1) 模型通常要求数据具有指数趋势,并且估计的参数必须能够合理地反映数据的动态变化。
应用领域:
适用问题:
优点:
缺点: