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

巴山学长

发布于 2021-04-22 07:12:28

3.7K00

代码可运行

文章被收录于专栏：巴山学长巴山学长

运行总次数：0

代码可运行

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

GM(1,1)编程步骤：

1.建立时间序列

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

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

4.计算邻均值等权数列

5.构造矩阵B、Y

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

7.计算模型拟合值

8.模型精度评定（后验差检验）

①计算残差

②计算标准差

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

④查表定级

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);

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);