灰色预测GM(1,1)-Matlab实现

灰色预测是通过少量的，不完全的信息，建立数学模型并给出预测的一种预测方法

是处理小样本预测问题的有效工具

适用条件

 (1)数据是以年份度量的非负数据（如果是月份或者季度数据一般要用时间序列模型）

 (2)数据能通过准指数规律的检验（确定原始数据是否可以使用灰色预测模型）

 (3)数据的期数较短（>3&&<=10）且其他数据之间的关联性不强（如果期数较长，一般用传统的时间序列模型比较合适）

GM(1,1)模型

只含有一个变量的一阶微分方程模型

模型评价

检验模型对原始数据的拟合程度

 代码

gm11脚本（函数文件）

function[result,x0_hat,relative_residuals,eta]=gm11(x0,predict_num)
n=length(x0);
x1=cumsum(x0);
z1=(x1(1:end-1)+x1(2:end))/2;
y=x0(2:end);x=z1;
k=((n-1)*sum(x.*y)-sum(x)*sum(y))/((n-1)*sum(x.*x)-sum(x)*sum(x));
b=(sum(x.*x)*sum(y)-sum(x)*sum(x.*y))/((n-1)*sum(x.*x)-sum(x)*sum(x));
a=-k;
x0_hat=zeros(n,1);x0_hat(1)=x0(1);
for m=1:n-1
    x0_hat(m+1)=(1-exp(a))*(x0(1)-b/a)*exp(-a*m);
end
result=zeros(predict_num,1);
for i=1:predict_num
    result(i)=(1-exp(a))*(x0(1)-b/a)*exp(-a*(n+i-1));
end
absolute_residuals=x0(2:end)-x0_hat(2:end);
relative_residuals=abs(absolute_residuals)./x0(2:end);
class_ratio=x0(2:end)./x0(1:end-1);
eta=abs(1-(1-0.5*a)/(1+0.5*a)*(1./class_ratio));
end

主脚本

%%
%输入原始数据并做出时间序列图
clear;clc
year=[2019:1:2023]';%加'表示转置，将横向量转为列向量
x0=[5503,5222,5429,5119,5175]';

%画出原始数据的时间序列图
figure(1);%因为我们的图形不止一个，因此要设置编号
plot(year,x0,'o-');grid on;%plot是用于绘制二维线图的函数；%grid on打开图形的网格线，网格线有助于在图表上更加清晰地看到数据点和趋势
set(gca,'xtick',year(1:1:end))%set(gca, 'xtick', ...) 用于设置当前图形的 x 轴刻度，gca 表示获取当前坐标轴；%year(1:1:end) 表示从 year 向量中提取所有的年份数据。1:1:end 是一种索引方式，表示从第一个元素到最后一个元素，间隔为 1。所以 year(1:1:end) 实际上就是 year 向量的所有元素。
xlabel('Years');ylabel('Number of Pet Cats in China (in 10,000s)');

%%
%对一次累加后的数据进行准指数规律的检验
disp('准指数规律检验')
n=length(x0);%计算原始数据的长度
x1=cumsum(x0);%cumsum是累加函数
rho=x0(2:end)./x1(1:end-1);%计算光滑度rho
%画出光滑度的图形，并画上0.5的直线，表示临界值
figure(2)
plot(year(2:end),rho,'o-',[year(2),year(end)],[0.5,0.5],'-');grid on;
text(year(end-1)+0.2,0.55,'临界线')%在坐标(year(end-1)+0.2,0.55)上添加文本
set(gca,'xtick',year(2:1:end))
xlabel('年份');ylabel('原始数据的光滑程度');

disp(strcat('指标1：光滑比小于0.5的数据占比为',num2str(100*sum(rho<0.5)/(n-1)),'%'))
disp(strcat('指标2：除去前两个时期外，光滑比小于0.5的数据占比为',num2str(100*sum(rho(3:end)<0.5)/(n-3)),'%'))
disp('参考指标1一般要大于60%，指标2要大于90%，你认为本例数据可以通过检验吗？')
Judge=input('你认为可以通过准指数规律的检验吗？可以通过请输入1，进行进一步的计算：');
%%
%将数据分为训练组和实验组（根据原数据量大小n来取，n<7则取最后2年作为试验组，n>7则取最后3年作为实验组）
if n>7
    test_num=3;
else
    test_num=2;
end
train_x0=x0(1:end-test_num);%训练数据
disp('训练的数据是：')
disp(mat2str(train_x0'))%mat2str可以将矩阵或向量转换为字符串显示
test_x0=x0(end-test_num+1:end);%试验数据
disp('试验数据是：')
disp(mat2str(test_x0'))
%使用GM(1,1)模型对训练数据进行训练
disp(' ')
result1=gm11(train_x0,test_num);
disp(' ')
%绘制对试验数据进行预测的图形
test_year=year(end-test_num+1:end);
figure(3)
plot(test_year,test_x0,'o-',test_year,result1,'*-');grid on;
set(gca,'xtick',year(end-test_num+1):1:year(end))
legend('试验组的真实数据','GM(1,1)预测结果')
xlabel('Years');ylabel('Number of Pet Cats in China (in 10,000s)')

predict_num=input('请输入你要往后面预测的期数：');
[result,x0_hat,relative_residuals,eta]=gm11(x0,predict_num);
%%
%输出使用模型预测出来的结果
disp('对原始数据拟合效果：')
for i=1:n
    disp(strcat(num2str(year(i)),': ',num2str(x0_hat(i))))
end
disp(strcat('往后预测',num2str(predict_num),'期的得到的结果：'))
for i=1:predict_num
    disp(strcat(num2str(year(end)+i),': ',num2str(result(i))))
end
%%
%绘制相对残差和级比偏差的图形
figure(4)
subplot(2,1,1)
plot(year(2:end),relative_residuals,'*-');grid on;%相对残差relative_residuals
legend('相对残差');xlabel('年份');
set(gca,'xtick',year(2:1:end))
subplot(2,1,2)
plot(year(2:end),eta,'o-');grid on;%级比偏差eta
legend('级比偏差');xlabel('年份');
set(gca,'xtick',year(2:1:end))
disp(' ')
disp('下面将输出对原数据拟合的评价结果')
%%
%残差检验
average_relative_residuals=mean(relative_residuals);
disp(strcat('平均相对残差为',num2str(average_relative_residuals)))
if average_relative_residuals<0.1
    disp('残差检验的结果表明：该模型对原数据的拟合程度非常不错')
elseif average_relative_residuals<0.2
    disp('残差检验的结果表明：该模型对原数据的拟合程度达到一般要求')
else
    disp('残差检验的结果表明：该模型对原数据的拟合程度不太好，建议使用其他模型进行预测')
end
%%
%级比偏差检验
average_eta=mean(eta);
disp(strcat('平均级比偏差为',num2str(average_eta)))
if average_eta<0.1
    disp('级比偏差检验的结果表明：该模型对原数据的拟合程度非常不错')
elseif average_eta<0.2
    disp('级比偏差检验的结果表明：该模型对原数据的拟合程度达到一般要求')
else
    disp('级比偏差检验的结果表明：该模型对原数据的拟合程度不太好，建议使用其他模型进行预测')
end
%%
%绘制最终的预测效果图
figure(5)
plot(year,x0,'-o',year,x0_hat,'-*m',year(end)+1:year(end)+predict_num,result,'-*b');grid on;%m洋红色 b蓝色
hold on;
plot([year(end),year(end)+1],[x0(end),result(1)],'-*b')
legend('Raw data','Fitting data','Predicted data')
set(gca,'xtick',[year(1):1:year(end)+predict_num])
xlabel('Years');ylabel('Number of Pet Dogs in China (in 10,000s)');