灰色预测GM(1,1)方法步骤与MATLAB代码实现

最新推荐文章于 2024-09-25 23:17:36 发布

Dlemolve

最新推荐文章于 2024-09-25 23:17:36 发布

阅读量1.2k

点赞数 13

文章标签： matlab 算法

本文链接：https://blog.csdn.net/Dlemolve/article/details/140571258

版权

1. 前言

预测在日常生活和决策过程中扮演着关键角色，它能够帮助我们做出更明智的选择和规划未来的行动。通过预测，我们能够理解可能的发展趋势和潜在结果，从而在面对不确定性时减少风险，优化资源分配，提升效率，以及为未来做出更合适的准备。预测不仅仅是一种工具，更是决策过程中的重要参考，使我们能够更有自信地面对未知的挑战和机会。

在预测模型的选择上，现阶段常见的模型有神经网络模型、回归预测模型和灰色预测模型等。其中，灰色预测模型可以对较少数量的样本进行预测，且不需要考虑数据的分布规律和变化趋势。本文采用灰色预测 GM(1,1)模型进行预测，其思想是对原始数据累加生成数据序列，并通过微分函数建立预测模型，通过累加和均值序列求出原始数据在未来的变化情况。

GM(1,1)也是本蛆毕业论文所使用的预测方法，用来分析供应链风险的变化趋势。这个方法并不复杂，不然本蛆也不会写在论文里面。学不会也妹有关系，文末为您附上了MATLAB代码，换个数字就能直接使用，彻底解放大脑。

2.基础概念

2.1 定义

GM(1,1)是由邓聚龙(1987)首次提出的一种灰色预测模型，它可以有效处理具有指数增长或衰减趋势的时间序列数据。该模型通过对原始数据序列进行累加生成一次累加序列，并利用一次指数拟合对累加序列进行模型建立，从而预测原始序列的未来发展趋势。GM(1,1)模型广泛应用于各种领域的数据分析和预测，特别适用于短期数据和变化规律不明显的情况下的预测工作。

2.2 特点

GM(1,1)具有以下特点：

（1）适合于小样本数据的预测。

（2）不需要考虑数据的分布规律和变化趋势。

（3）缺点在于预测结果可能不够精确。

3. 方法步骤

Step.1 构建累加序列

假设总共存在 𝑚 项影响供应链的风险因素指标，对每一项指标采样 𝑛 个值，建立原始不完全信息非负序列

$X^{(0)}=\left(x^{(0)}\begin{pmatrix}1\end{pmatrix},x^{(0)}\begin{pmatrix}2\end{pmatrix},\cdots,x^{(0)}\begin{pmatrix}n\end{pmatrix}\right)$ 。

对 $x^{(0)}(k)$ 进行一次累加，生成一次累加序列(1-Accumulating generation operational, 1-AGO)

$X^{(1)}=\Big(x^{(1)}\big(1\big),x^{(1)}\big(2\big),\cdots,x^{(1)}\big(n\big)\Big),\quad x^{(1)}(k)=\sum_{r=1}^kx^{(0)}\big(r\big)$ 。

其中， $k=1,2,\cdots,n$ 。同时，1-AGO序列满足 $x^{(1)}\left(1\right)=x^{(0)}\left(1\right)$ 。

由此可得多指标的 1-AGO 累加矩阵

$X=\begin{bmatrix}x_1^{(1)}\left(1\right)&x_1^{(1)}\left(2\right)&\cdots&x_1^{(1)}\left(n\right)\\x_2^{(2)}\left(2\right)&x_2^{(2)}\left(2\right)&\cdots&x_2^{(2)}\left(n\right)\\\cdots&\cdots&\cdots&\cdots\\x_n^{(m)}\left(1\right)&x_n^{(m)}\left(2\right)&\cdots&x_n^{(m)}\left(n\right)\end{bmatrix}$

Step.2 生成均值序列

根据1-AGO累加矩阵，得到各项指标均值

$z^{(1)}\big(k\big)=\frac{1}{2}\Big(x^{(1)}\big(k\big)+x^{(1)}\big(k-1\big)\Big)$ 。

由此生成均值序列

$Z^{(1)}=\left(z^{(1)}(2),z^{(1)}(3),\cdots,z^{(1)}(n)\right)$ 。

Step.3 计算模型参数

GM(1,1)模型可表示为 $x^{(0)}(k)+az^{(1)}(k)=b$ 。

其中 𝑎 与 𝑏 为GM(1,1)模型的两个被识别的参数。

根据上式，GM(1,1)模型的白化微分方程为

$\frac{dx^{(1)}(k)}{dt}+ax^{(1)}(k)=b$ 。

令 $\hat{a}=\begin{bmatrix}a,b\end{bmatrix}^T$ 为参数序列，应用最小二乘法对其求解

$\hat{a}=\begin{bmatrix}a,b\end{bmatrix}^T=\begin{pmatrix}\mathbf{B}^T\mathbf{B}\end{pmatrix}^{-1}\mathbf{B}^T\mathbf{Y}$ 。

其中 $\mathbf{B}$ 和 $\mathbf{Y}$ 分别为常数向量和累加矩阵，且

$\mathbf{B}=\begin{pmatrix}-Z_i^{(2)}&1\\-Z_i^{(3)}&1\\\vdots&\vdots\\-Z_i^{(n)}&1\end{pmatrix},\quad\mathbf{Y}=\begin{pmatrix}x_1^{(0)}\left(2\right)\\x_2^{(0)}\left(3\right)\\\vdots\\x_n^{(0)}\left(n\right)\end{pmatrix}$ 。

Step.4 求解预测值

将求解的参数值代入白化微分方程，对其求解后得到时间响应序列

$\hat{x}^{(1)}\left(k\right)=(x^{(0)}(1)-\frac ba)e^{-a(k-1)}+\frac ba$ 。

最后，将 𝑥(1)(𝑘) 还原成原始数列即得到预测函数

$\hat{x}^{(0)}\left(k\right)=\hat{x}^{(1)}\left(k\right)-\hat{x}^{(1)}\left(k-1\right)$ 。

Step.5 验证预测结果

在预测之前，对原始序列进行级比检验，判断序列是否能够使用 GM(1,1)模型进行预测。级比的计算公式为

$\vartheta(k)=\frac{x^{(0)}\left(k-1\right)}{x^{(0)}\left(k\right)}$ 。

如果级比 𝜗(𝑘) 在区间 $\left[e^{-2/\left(n+1\right)},e^{2/\left(n+2\right)}\right]$ 内，则说明序列的预测可行，否则预测结果无效。

注：为节省篇幅，该步骤只采用残差检验对模型进行验证，检验结果的精度和可靠性。

将 $x^{(0)}(k)$ 与 $\hat{x}^{(0)}(k)$ 进行对比，得到残差

$f^{(0)}(k)=x^{(0)}(k)-\hat{x}^{(0)}(k),k=2,3,...,n$ 。

则相对误差为

$R(k)=(x^{(0)}(k)-\hat{x}^{(0)}(k))/x^{(0)}(k),k=2,3,...,n$ 。

所对应的平均绝对百分比误差(Mean Absolute Percentage Error, MAPE)为

$MAPE=\frac1{n-1}\sum_{k=1}^n\frac{f^{(0)}(k)}{x^{(0)}(k)}$ 。

一般地，MAPE的精度等级参照如下。

检验结果	精度等级
MAPE < 10%	非常高(Highly accurate)
10% ≤ MAPE < 20%	较高(Good)
20% ≤ MAPE < 50%	一般(Reasonable)
MAPE ≥ 50%	较低(Inaccurate)

4. 算例分析

据传奇名著《斗破苍穹》记载，萧炎吞噬青莲地心火时难以承受其炎压，导致自身游离在生死边缘。此时药老提醒他需要将异火用斗气不断炼化，以减少异火对身体造成的伤害。药老从纳戒里掏出了盖革计数器用来测量萧炎体内未炼化的异火值，仪器显示，从凌晨1:00到4:00未炼化的异火值分别为[5.8, 5.6, 5.2, 4.9] (单位：µSv/h)。虽然数值在不断下降，但只有异火值低于4.0µSv/h时人体才会脱离危险。请通过GM(1,1)模型，预测在早上7:00时萧炎是否脱离危险。

Step.1

首先根据以上信息，建立累加序列 $X^{(0)}=(5.8000, 5.6000,5.2000,4.9000)$ 。

接着，对 $X^{(0)}$ 进行级别检验，检验结果显示结果符合要求。

例： $\vartheta(2)=\frac{x^{(0)}\left(1\right)}{x^{(0)}\left(2\right)}=0.9655$ ， $e^{-2/\left(2+1\right)}\leq0.9655 \leq e^{2/\left(2+2\right)}$ ，则当 𝑛=2 时级别检验通过。

则1-AGO序列为 $X^{(1)}=(5.8000, 11.4000, 16.6000, 21.5000)$ 。

例：16.6000=5.8000+5.6000+5.2000。

Step.2

根据1-AGO累加矩阵的指标均值，得到均值序列 $Z^{(1)}=(8.6000 , 14.0000 , 19.0500)$ 。

例：8.6000=0.5×(5.8000+11.4000)。

Step.3

设置常数向量和累加矩阵。

$\mathbf{B}=\begin{pmatrix}-8.6000&1\\-14.000&1\\-19.000&1\end{pmatrix},\quad\mathbf{Y}=\begin{pmatrix}5.6000\\5.2000\\4.9000\end{pmatrix}$

可以计算求得参数 𝑎=0.0671 ， 𝑏=6.1645 。

计算过程略。

Step.4

将参数 𝑎 和 𝑏 代入白化微分方程，对其求解后得到时间响应序列

$\hat{x}^{(1)}=(5.8000 , 11.3860 , 16.6097 , 21.4946 ,26.0626 , 30.3342 , 34.3288 )$ 。

例：11.3860=(5.8000-(0.0671÷6.1645))×EXP(-0.0671×(2-1))+(0.0671÷6.1645)。

对时间响应序列进行还原，得到预测结果

$\hat{x}^{(0)}=(5.8000 , 5.5860 , 5.2237 , 4.8848 , 4.5680 , 4.2717 , 3.9946)$ 。

例：5.5860=11.3860-5.8000。

Step.5

对预测结果进行验证。

残差序列为 $f^{(0)}=(0 , 0.0140 , -0.0237, 0.0152)$ 。

例：0.0140=5.6000-5.5860。

相对误差为 $R=(0 , 0.0025 , 0.0046 , 0.0031)$ 。

例：0.0025=0.0140÷5.6000。

最后，得到 $MAPE=0.1690$ 。

计算过程略。

5. MATLAB代码

本方法的MATLAB代码如下。

clear;
clc;
% 作者：
% 如蛆虫般遨游（知乎）
% Dlemolve（CSDN）
% 阴间偷渡客（Bilibili）
% 代码免费分享，拿去赚钱似劝架。
format short
% 【Determine the original sequence.】
X0=[5.8	 5.6  5.2  4.9];
display(X0)
% 【Test the ratio level.】
Judge_1=exp((-2)/(length(X0)+1));
Judge_2=exp(2/(length(X0)+2));
for n=(2:length(X0))
    if Judge_1<=X0(n)/X0(n-1)<=Judge_2
        disp('Qualified')
    else disp('Unqualified')
    end
end
% 【Calculate the 1-AGO(Accumulating Generation Operational) sequence.】
for n=(1:length(X0))
    if n==1
        X1(n)=X0(n);
    else X1(n)=X0(n)+X1(n-1);
    end
end
display(X1)
% 【Compute the consecutive neighbor sequence of X1.】
Gamma=0.5;
for n=(1:length(X1)-1)
    Z1(n)=Gamma*X1(n+1)+(1-Gamma)*X1(n);
end
display(Z1)
% 【Set the parameters by least square method.】
Value_Token=1;
Value=repmat(Value_Token,1,length(Z1));
B=[-Z1' Value'];
Y=X0(2:length(X0))';
display(B)
display(Y)
% Construct the GM(1,1) model by least square method: X0(n)+aZ1(n)=b,
% and then get the parameters of a and b.
ab_Token=(B'*B)^(-1)*B'*Y;
a=ab_Token(1);
b=ab_Token(2);
display(a)
display(b)
% 【Get the time response sequence.】
Expected=7;  % It can be changed according to the DM.
Xt1=[ ];
Xt1(1)=X0(1);
for n=(2:Expected)
    Xt1(n)=(X0(1)-(b/a))*exp(-a*(n-1))+(b/a);
end
display(Xt1)
% 【Obtain the estimated value of the original sequence.】
Xt0=[ ];
Xt0(1)=X0(1);
for n=(2:Expected)
    Xt0(n)=Xt1(n)-Xt1(n-1);
end
display(Xt0)
% 【Verify the prediction result.】
% Compute the residuals.
Residuals=X0-Xt0(1:length(X0));
display(Residuals)
% Compute the relative residuals.
Relative_error=abs(Residuals./X0);
display(Relative_error)
% Calculate MAPE.
MAPE=abs(sum(Relative_error)/(n-1))*100;
display(MAPE)

注：

（1）作者并非为编程专业人士，代码可能存在优化空间或纰漏，欢迎各位读者提出改进策略，非常感谢。

（2）代码免费分享，拿去赚钱死欠佳。

6. 参考文献

[1] Yao T , Wang Z . Crude oil price prediction based on LSTM network and GM (1,1) model[J]. Grey Systems Theory and Application, 2020, ahead-of-print(ahead-of-print).

[2] Liu S , Yin C , Cao D . Weapon equipment management cost prediction based on forgetting factor recursive GM (1,1) model[J]. Grey Systems Theory & Application, 2019, 10(1):38-45.

感谢您的阅读，祝您科研有成、收获丰硕。若发现文章存在错误，恳请您作出批评与指正。

本文由“如蛆虫般遨游”学术交流账号原创，如有侵权请联系删除。如需转载，请联系作者。