Matlab-灰色预测模型

Matlab随笔 2021-08-25

919

灰色系统理论是一种研究少数据、贫信息不确定性问题的新方法，灰色预测模型是灰色系统理论的一个重要研究领域。其中，GM系列模型是基本模型，尤其是GM(1,1)模型，应用十分广泛。传统GM(1,1)模型建模对象主要是具有指数或近似指数增长规律的序列，对应非单调增长数据，模拟精度较差。

对于非单调的摆动发展序列或有饱和的S形序列，可以考虑建立GM(2,1)。

一.GM(1,1)

1.理论

1.建立时间序列

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

当采用GM(1,1)模型的各种形式进行模拟精度达不到要求时，可以考虑对残差序列建立GM(1,1)模型，对原来的模型进行修正，以提高模拟精度。

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

4.计算邻均值等权数列

5.构造矩阵B、Y

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

7.计算模型拟合值

8.模型精度评定

①方法1：残差检验

②方法2：后验差检验

③方法3：小误差概率检验

9.查表定级

2.代码：

方法1：（利用公式）

clear;clc;
% 建立时间序列
x0 = [89677 99215 109655 120333 135823 159878 182321 209407 246619]';
% 需要预测几期数据
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);
p = length(find(e-mean(e)<0.6745*s1))/n1;
if p>=0.95
    disp('精度等级:1级(好)');
elseif p>=0.80
    disp('精度等级:2级(合格)');
elseif p>=0.7
    disp('精度等级:3级(勉强)');
else
    disp('精度等级:4级(不合格)');
end
% 绘图说明
hold on
p1=plot(1:n1,x0,'ro','LineWidth',2.5);
p2=plot(1:n1+count,x0_hat,'bo','LineWidth',1.5);
p3=plot(1:n1+count,x0_hat,'b-','LineWidth',1.5);
legend([p1,p2],'实测值','预测值','FontSize',18);
title('GM(1,1)');
hold off

方法2：（利用微分方程）

clear;clc;
% 建立时间序列
x0 = [89677 99215 109655 120333 135823 159878 182321 209407 246619]';
% 需要预测几期数据
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);
% 计算u=[发展系数a,灰作用量b]
u = (B'*B)\(B'*Y);
% 计算模型拟合值
syms x(t)
x=dsolve(diff(x)+u(1)*x == u(2),x(0)==x1(1));
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);
p = length(find(e-mean(e)<0.6745*s1))/n1;
if p>=0.95
    disp('精度等级:1级(好)');
elseif p>=0.80
    disp('精度等级:2级(合格)');
elseif p>=0.7
    disp('精度等级:3级(勉强)');
else
    disp('精度等级:4级(不合格)');
end
% 绘图说明
hold on
p1=plot(1:n1,x0,'ro','LineWidth',2.5);
p2=plot(1:n1+count,x0_hat,'bo','LineWidth',1.5);
p3=plot(1:n1+count,x0_hat,'b-','LineWidth',1.5);
legend([p1,p2],'实测值','预测值','FontSize',18);
title('GM(1,1)');
hold off

二.GM(2,1)

1.理论

1.建立时间序列

2.计算一次累加生成序列、一次累减生成序列

3.计算邻均值等权数列

4.构造矩阵B、Y

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

6.求解白化方程

7.计算模型拟合值

8.模型精度评定同GM(1,1)

2.代码

clear;clc;
% 建立时间序列
x=-10:8;
x0=(x.^2)';
n1 = length(x0);
% 需要预测几期数据
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);
p = length(find(e-mean(e)<0.6745*s1))/n1;
if p>=0.95
    disp('精度等级:1级(好)');
elseif p>=0.80
    disp('精度等级:2级(合格)');
elseif p>=0.7
    disp('精度等级:3级(勉强)');
else
    disp('精度等级:4级(不合格)');
end
% 绘图说明
hold on
p1=plot(1:n1,x0,'ro','LineWidth',2.5);
p2=plot(1:n1+count,x0_hat,'bo','LineWidth',1.5);
p3=plot(1:n1+count,x0_hat,'b-','LineWidth',1.5);
legend([p1,p2],'实测值','预测值','FontSize',18);
title('GM(2,1)');
hold off