灰色预测模型在matlab数据预测中的应用【编程算法】

5259754a7f0f3bfa5bc45f06ed425e65.png

  概述算法:灰色预测模型用于对原始数据(≥4个)做中短期预测,其中,GM(1,1)模型适用于具有较强的指数规律的序列,只能描述单调的变化过程,而GM(2,1)模型适用于非单调的摆动发展序列或具有饱和的S形序列。

GM(1,1)编程步骤:

1.建立时间序列

   ac2959ababbb32b36a846a138bb1d6ac.png

       2.检验数据是否符合要求

         102580a9e42da3d5ee0a4cd6b2d42444.png

       3.计算一次累加生成序列

          5a6bf93c6427a6219b64d7bfe2be0816.png

       4.计算邻均值等权数列

         3213a71cabdb117dc76f62710248122d.png

       5.构造矩阵B、Y

          04f01256d2d63c7d330ff5a6ce73be41.png

       6.计算发展系数a和灰作用量b

          ccd426893b6f97bb81758543dc117d54.png

       7.计算模型拟合值

         7f3ac526ce94df3e76453e4c53bcb7a6.png

       8.模型精度评定(后验差检验)

       ①计算残差

         e7c64e0275d2f0cb03ad8a26e786fc01.png

       ②计算标准差

          c010387f0b01a917f1933436d7f48eda.png

       ③计算后验差比值、小误差概率

          f54192bdf39d9d6d4f3897b54d2ad529.png

④查表定级

        b4166939e93db05ede0ad9cf88fa74ac.png

GM(2,1)编程步骤与GM(1,1)类似。

下面就一起来看看如何将优雅的数学语言转换成matlab语言吧。

GM(1,1)源代码

clear;clc;
% 建立时间序列【输入】
x0 = [15.9 15.4 18.1 21.3 20.1 22.0 22.6 21.4]';
% 需要预测几期数据【输入】,预测数据见x0_hat变量
count = 2;
% 检验数据是否符合要求
n1 = length(x0);
lmd = x0(1:end-1)./x0(2:end);
flag = min(lmd)>exp(-2/(n1+1)) & max(lmd)<exp(2/(n1+1));
if ~flag
    error('数据级比不符合要求');
end
% 计算一次累加生成序列
x1 = cumsum(x0);
% 计算邻均值等权数列
z1 = 0.5 * ( x1(1:end-1) + x1(2:end) );
% 构造矩阵B、Y
B = [-z1,ones(n1-1,1)];
Y = x0(2:n1);
% 计算发展系数a和灰作用量b
u = (B'*B)\(B'*Y);
a = u(1);
b = u(2);
% 计算模型拟合值
k = (1:n1-1+count)';
x0_hat = [x0(1);(1-exp(a))*(x0(1)-b/a)*exp(-a*k)];
disp('预测数据:')
x0_hat(n1+1:end)
% 模型精度评定
e = x0(1:n1)-x0_hat(1:n1);
s1 = std(x0);
s2 = std(e);
c = s2/s1;
p = length(find(e-mean(e)<0.6745*s1))/n1;
if p>=0.95 && c<=0.35
    disp('精度等级:1级(好)');
elseif p>=0.80 && c<=0.5
    disp('精度等级:2级(合格)');
elseif p>=0.7 && c<=0.65
    disp('精度等级:3级(勉强)');
else
    disp('精度等级:4级(不合格)');
end
% 绘图说明
plot(1:n1,x0,'ro','LineWidth',2.5);
hold on
plot(1:n1+count,x0_hat,'bo','LineWidth',2.5);
plot(1:n1+count,x0_hat,'b-','LineWidth',2.5);
hold off
legend('实测值','预测值','FontSize',18);

003aa2d48eb5a0128f13f55b552055fa.png

GM(2,1)代码

clear;clc;
% 建立时间序列【输入】
x0 = [5.6 4.2 3.3 2.5 3.1 4.4 5.8]';
n1 = length(x0);
% 需要预测几期数据【输入】,预测数据见x0_hat变量
count = 2;
% 计算一次累加生成序列
x1 = cumsum(x0);
% 计算一次累减生成序列
alpx0 = x0(2:end)-x0(1:end-1);
% 计算邻均值等权数列
z1 = 0.5 * ( x1(1:end-1) + x1(2:end) );
% 构造矩阵B、Y
B = [ -x0(2:end),-z1,ones(n1-1,1)];
Y = alpx0;
% 计算发展系数a和灰作用量b
u = (B'*B)\(B'*Y);
% 计算模型拟合值
syms x(t)
x=dsolve(diff(x,2)+u(1)*diff(x)+u(2)*x==u(3),x(0)==x1(1),x(n1-1)==x1(n1));
xt=vpa(x,2);
x1_hat=subs(x,t,0:n1-1+count);
x1_hat=(double(x1_hat))';
x0_hat = [x0(1);diff(x1_hat)];
disp('预测数据:')
x0_hat(n1+1:end)
% 模型精度评定
e = x0(1:n1)-x0_hat(1:n1);
s1 = std(x0);
s2 = std(e);
c = s2/s1;
p = length(find(e-mean(e)<0.6745*s1))/n1;
if p>=0.95 && c<=0.35
    disp('精度等级:1级(好)');
elseif p>=0.80 && c<=0.5
    disp('精度等级:2级(合格)');
elseif p>=0.7 && c<=0.65
    disp('精度等级:3级(勉强)');
else
    disp('精度等级:4级(不合格)');
end
% 绘图说明
plot(1:n1,x0,'ro','LineWidth',2.5);
hold on
plot(1:n1+count,x0_hat,'bo','LineWidth',2.5);
plot(1:n1+count,x0_hat,'b-','LineWidth',2.5);
hold off
legend('实测值','预测值','FontSize',18);

b2ee6867c8ebdfc5034b7359932cebd3.png

通过学习相关算法并将算法转变为实际的编程语言是练习编程的一种重要途径,这不仅可以提升理论认知,还能提高实践动手能力。鉴于此,matlab爱好者公众号计划推出【编程算法】系列,将逐一介绍各类算法在matlab中实现,与大家一起来在算法的海洋里畅游。

若您对算法感兴趣,并有一定的matlab编程基础,欢迎将所学算法整理成文推送给我们。

欢迎大家向matlab爱好者公众号投稿!

投稿请加群

参考资料:

[1]《数学建模算法与应用(第2版)—司守奎

封面图片:由 mohamed Hassan 在Pixabay上发布

  • 0
    点赞
  • 5
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值