【数模百科】一文讲清楚灰色预测模型GM（2，1）的原理（附python代码）

本文主要是介绍【数模百科】一文讲清楚灰色预测模型GM（2，1）的原理（附python代码），希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

本篇文章摘录自GM(2,1) - 数模百科，如果想了解更多有关灰色预测模型的信息，请移步灰色预测模型 - 数模百科

在学习GM（2，1）之前，强烈建议阅读这篇文章，这里包含了数据分析与检验，和灰色预测模型的基础知识

【数模百科】一篇文章讲清楚灰色预测模型GM（1，1）附python代码-CSDN博客

你已经了解了GM(1,1)模型，它是用来处理和预测只有一个因素（变量）随时间变化情况的模型，而且它是基于一阶微分方程的。这个模型假设数据之间的增长率是恒定的，但是在实际情况中，很多时候数据的增长率并不是一成不变的。

这时候，灰色模型中的GM(2,1)就派上用场了，GM(2,1)模型仍然是针对单一变量的，只不过它使用的是二阶微分方程。那个“2”表示的是模型预测时考虑了数据变化率的变化率，也就是数据增长速度的速度。这样就能让模型更加精细地捕捉数据随时间变化的趋势，尤其是当数据变化不是那么规则，或者有点起起伏伏时。

简单来说，如果你觉得GM(1,1)模型对某些数据的预测不够精确，或者数据的走势有较大的波动，你可能会想尝试用GM(2,1)模型，它通过考虑更复杂的变化趋势来提供可能更精确的预测结果。因此，GM(2,1)模型适用于那些数据变化比较复杂，一阶模型难以准确捕捉其变化趋势的情况。

定义与详解

以下是GM(2,1)模型的构建和求解过程。

原始数据序列：

首先，我们有一个原始数据序列，它是时间序列中的观察值
$x^{(0)}=\{x^{(0)}(1),x^{(0)}(2),\ldots,x^{(0)}(n)\},$
其中 $x^{(0)}(t)$ 表示第 t 个时刻的观测值。
一次累加生成操作（I-AGO）序列 $x^{(1)}$ ：

通过对原始数据序列进行累加，我们得到一次累加序列，这有助于展现数据的增长趋势和结构：
$x^{(1)}(k) = \sum_{i=1}^{k} x^{(0)}(i), \quad k = 1, 2, ..., n .$
得到
$x^{(1)}=\{x^{(1)}(1),x^{(1)}(2),\ldots,x^{(1)}(n)\}.$
生成一次累加逆生成操作（I-IAGO）序列 $\alpha^{(1)}x^{(0)}$ ，是 I-AGO 序列的逆操作，它的目的是从累加序列中重构原始数据序列。基本上，逆生成操作是对累加序列 $x^{(1)}$ 进行差分，以恢复原始的数据值。
$\alpha^{(1)}x^{(0)}=(\alpha^{(1)}x^{(0)}(2),\dots,\alpha^{(1)}x^{(0)}(n)).$
其中，
$\alpha^{(1)}x^{(0)}(k)=x^{(0)}(k)-x^{(0)}(k-1),k=2,3,\dots,n$
生成 $x^{(1)}$ 的紧邻均值生成序列 $z^{(1)}$ ，紧邻均值生成序列 $z^{(1)}$ 是通过计算 $x^{(1)}$ 序列中相邻两个数的平均值得到的，它的目的是用来构建模型方程：
$z^{(1)}=(z^{(1)}(2),z^{(1)}(3),\dots,z^{(1)}(n)).$
建立GM(2,1)模型：
$\alpha^{(1)}x^{(0)}(k)+a_1x^{(0)}(k)+a_2z^{(1)}(k)=b.$
这个方程中的 $a_1$ , $a_2$ , 和 b 是我们需要通过数据来估计的参数。
建立微分方程模型：基于一次累加序列，我们构建一个包含未知参数的二阶微分方程，即GM(2,1)的白化方程。
$\frac{d^2 x^{(1)}}{dt^2} + a_1 \frac{dx^{(1)}}{dt} + a_2 x^{(1)} = b .$
其中， $a_1$ ， $a_2$ ，和 b 是模型参数，需要通过数据拟合来确定。
求解模型参数：为了估计模型参数，我们构建矩阵 B 和向量 Y ：
$B = \begin{bmatrix} -x^{(0)}(2) & -z^{(1)}(2) & 1 \\\ -x^{(0)}(3) & -z^{(1)}(3) & 1 \\\ \vdots & \vdots & \vdots \\\ -x^{(0)}(n) & -z^{(1)}(n) & 1 \end{bmatrix}, \quad Y = \begin{bmatrix} \alpha^{(1)} x^{(0)}(2) \\\ \alpha^{(1)} x^{(0)}(3) \\\ \vdots \\\ \alpha^{(1)} x^{(0)}(n) \end{bmatrix} .$
参数通过最小二乘估计得到：
$\hat{a} = (B^T B)^{-1} B^T Y .$
模型预测：有了模型参数后，我们可以通过解上述微分方程，使用这些参数对未来的数据点进行预测。
还原预测结果：将预测的累加序列 $x^{(1)}$ 通过逆累加生成操作转换回原始数据序列 $x^{(0)}$ ，得到预测的原始数据值，以便于与实际观察值进行比较。
$\hat{x}^{(0)}(k) = \hat{x}^{(1)}(k) - \hat{x}^{(1)}(k-1), \quad k=2,3,\dots,n .$

这个过程的核心是通过累加操作平滑数据，建立和求解微分方程来预测数据的趋势。GM(2,1)模型特别适合预测少量数据和具有一定发展趋势的时间序列。

代码

import numpy as np
from sympy import symbols, Function, Eq, diff, dsolve
from numpy.linalg import lstsqt = symbols('t')  # 定义符号变量 t
x = Function('x')(t)  # 定义函数 x(t)# 数据初始化
x0 = np.array([41,49,61,78,96,104]) # 输入原始序列
n = len(x0)  # 原始序列元素个数
x1 = np.cumsum(x0)  # 原始序列的累积和
a_x0 = np.diff(x0)  # 对原始序列进行一阶差分
a_x0 = np.concatenate([[0], a_x0])  # 在差分结果前面补0，保持序列长度一致z = np.zeros(n)  # 初始化z序列
for i in range(1, n):  # 对每一个元素：z[i] = 0.5 * (x1[i] + x1[i - 1])  # 计算相邻元素的平均值# 构建数组 B 和 Y，用于之后的方程求解
B = np.column_stack((-x0[1:], -z[1:], np.ones(n - 1)))  # 构建 B，其中包含负的 x0 序列元素（从第二个元素开始）、负的 z 序列元素（从第二个元素开始）以及单位列向量
Y = a_x0[1:]  # 构建 Y，由一阶差分序列 a_x0 （从第二个元素开始）# 求解方程，获得参数 u
u, _, _, _ = lstsq(B, Y, rcond=None)  # 运用 Numpy 的最小二乘法 lstsq 求解线性方程得到 u# 定义微分方程
a1, a2, b, = u[0], u[1], u[2]  # 从参数 u 提取所需参数 a1，a2 和 b
diff_eq = Eq(diff(x, t, t) + a1 * diff(x, t) + a2 * x, b)  # 定义微分方程# 设定初始条件
ics = {x.subs(t, 0): x1[0], x.subs(t, 5): x1[-1]}  # 设定初始条件为 x(0) = x1[0]，x(5) = x1[5] （注意这里索引是以 0 开始的）# 解微分方程
solution = dsolve(diff_eq, ics=ics)  # 解微分方程# 运用求解出的微分方程进行预测
yuce = [solution.rhs.subs(t, val).evalf() for val in range(n)]  # 对给出的数据序列进行预测
x0_hat = [yuce[0]] + np.diff(yuce).tolist()  # 对预测结果的序列计算一阶差分# 计算预测误差和相对误差
epsilon = x0 - x0_hat  # 计算预测值和实际值之间的差，即预测误差
delta = np.abs(epsilon / x0)  # 计算绝对比例误差，也就是相对误差

这个代码提供了一个基于 GM(2,1) 模型对时间序列数据进行预测的完整工作流程，包含了建模、求解、预测和可视化的各个步骤：

以上仅为模型的参考代码，如果你在实践过程中预测的值相比于原始数据变化幅度过大，说明原始数据可能并不适合该模型。

在实践中，GM(2,1) 模型通常要求数据具有指数趋势，并且估计的参数必须能够合理地反映数据的动态变化。