【学习笔记】Matlab实现GM(1,1)模型

参考资料：
清风：数学建模算法、编程和写作培训的视频课程以及Matlab等软件教学

概念

GM（1,1）模型是使用原始的离散非负数据列，通过一次累加生成削弱随机性、较有规律的离散数据列，然后通过建立微分方程模型，得到在离散点处的解经过累减生成的原始数据的近似估计值，从而预测原始数据的后续发展

准指数规律的检验

原理

求解部分

解析解

GM(1,1)模型的评价

残差检验

级比偏差检验

结合具体预测场景来判断临界值，一般认为<0.1则说明拟合效果较好，<0.2则说明拟合效果达到一般要求，

代码部分

主函数

%%  输入数据
clear;clc
Time=input('请输入数据时期：');
Data=input('请输入数据：');
if size(Time,1)==1
    Time=Time';
end
if size(Data,1)==1
    Data=Data';
end
n=size(Time,1);
Model={'全数据GM(1,1)','新信息GM(1,1)','新陈代谢(GM1,1)'};
GraphCount=1;   %   统计需要绘制的图的个数
%判断数据是否存在负值或小于3期
if sum(Data<0)>0
    error('数据存在负值，程序终止');
end
if n<4
    error('数据过少，无法进行灰色预测');
end
if n>10
    error('数据过多，无法进行灰色预测');
end

%画出原始数据的时间序列图
fprintf('\n');
disp('------------------------------------------------------------');
disp('下面开始绘制原始数据的时间序列图');
LabelX=input('请输入X轴标签：','s');
LabelY=input('请输入Y轴标签：','s');
figure(GraphCount);
GraphCount=GraphCount+1;
plot(Time,Data,'o-');
grid on;
set(gca,'xtick',Time(1:1:end));
xlabel(LabelX);ylabel(LabelY);
fprintf('\n');
disp('------------------------------------------------------------');

%%  准指数规律检验
AccData=cumsum(Data);
Smooth=Data(2:end,1)./AccData(1:end-1,1);

%   统计ρ(k)∈(0,0.5)的占比（指标1）
count=0;
for i=1:size(Smooth,1)
    if Smooth(i,1)<0.5
        count=count+1;
    end
end
SmoPer1=count/size(Smooth,1)*100;
disp(['指标1：光滑比小于0.5的数据占比为',num2str(SmoPer1),'%'])

%   统计ρ(k)∈(0,0.5)的占比（指标2）
count=0;
for i=3:size(Smooth,1)
    if Smooth(i,1)<0.5
        count=count+1;
    end
end
SmoPer2=count/(size(Smooth,1)-2)*100;
disp(['指标2：去除前两个时期后，光滑比小于0.5的数据占比为',num2str(SmoPer2),'%'])
Judge=input('该数据是否通过准指数规律检验？请判断（输入1为通过，输入0为不通过）：');
fprintf('\n');
disp('------------------------------------------------------------');
if Judge==0
    error('数据未通过准指数规律检验，无法进行灰色预测，程序终止');
end

fprintf('\n');
disp('------------------------------------------------------------');
if n==4
    N=input('请输入要预测的期数：');
    [Result_Original,RelaResi_Original,EtaOriginal]=Original(Data,N);
    [Result_New]=New(Data,N);
    [Result_Metabolism,RelaResi_Metabolism,EtaMetabolism]=Metabolism(Data,N);
    FinalResult=(Result_Original+Result_New+Result_Metabolism)/3;
else
    disp('由于数据过多，故取部分数据作为试验组，其余数据为训练组，以提高模型精度。');
    if n>7
        TrainNum=3;
    else
        TrainNum=2;
    end
    TrainData=Data(1:end-TrainNum,1);
    [Result_Original,RelaResi_Original,EtaOriginal,DataHat]=Original(TrainData,TrainNum);
    [Result_New]=New(TrainData,TrainNum);
    [Result_Metabolism]=Metabolism(TrainData,TrainNum);
    %   计算SSE
    SSE1=sum(Data(end-TrainNum+1:end)-Result_Original);
    SSE2=sum(Data(end-TrainNum+1:end)-Result_New);
    SSE3=sum(Data(end-TrainNum+1:end)-Result_Metabolism);
    %   选取误差最小的模型
    if SSE1<SSE2
        if SSE1<SSE3
            ModelType=1;
        else
            ModelType=3;
        end
    else
        if SSE2>SSE3
            ModelType=3;
        else
            ModelType=2;
        end
    end
    string=strcat('经计算，',Model(1,ModelType),'精度最高，故以此模型来预测未来数据');
    fprintf('\n');
    disp('------------------------------------------------------------');
    disp(char(string))
    
    %   预测结果
    N=input('请输入要预测的期数：');
    fprintf('\n\n');
    [FinalResult,RelaResi,Eta,DataHat]=Original(Data,N);
    if ModelType==2
        FinalResult=New(Data,N);
    else
        FinalResult=Metabolism(Data,N);
    end
    disp(['往后预测',num2str(N),'的结果为：']);
    disp(mat2str(FinalResult));
    
    %   残差检验
    AverageRelaResi=mean(RelaResi)*100;
    disp(['该模型的平均相对残差为',num2str(AverageRelaResi),'%']);
    if AverageRelaResi<10
        disp('经残差检验，该模型的预测精度较高');
    elseif AverageRelaResi<20
        disp('经残差检验，该模型的预测精度一般');
    else
        disp('经残差检验，该模型的预测精度较低，建议更换其他模型');
    end
    
    %   级比偏差检验
    AverageEta=mean(Eta)*100;
    disp(['该模型的平均级比偏差为',num2str(AverageEta),'%']);
    if AverageEta<10
        disp('经级比偏差检验，该模型的预测精度较高');
    elseif AverageEta<20
        disp('经级比偏差检验，该模型的预测精度一般');
    else
        disp('经级比偏差检验，该模型的预测精度较低，建议更换其他模型');
    end
    
    %   绘制最终数据图
    figure(GraphCount)
    plot(Time,Data,'-o',  Time,DataHat,'-*m',  Time(end)+1:Time(end)+N,FinalResult,'-*b' );
    grid on;
    hold on;
    plot([Time(end),Time(end)+1],[Data(end),FinalResult(1)],'-*b')
    legend('原始数据','拟合数据','预测数据') 
    set(gca,'xtick',[Time(1,1):1:Time(end,1)+N])
    xlabel(LabelX);  ylabel(LabelY);
end
            
            

全数据模型

function [Result,RelativeResidual,Eta,DataHat]=Original(Data,N)
n=size(Data,1);
AccData=cumsum(Data);
Z=(AccData(1:end-1,1)+AccData(2:end,1))/2;
Data1=Data(2:end,1);
%%   计算灰作用量及发展系数
k=((n-1)*sum(Data1.*Z)-sum(Z)*sum(Data1))/((n-1)*sum(Z.*Z)-(sum(Z)*sum(Z)));
b=(sum(Z.*Z)*sum(Data1)-sum(Z)*sum(Z.*Data1))/((n-1)*sum(Z.*Z)-sum(Z)*sum(Z));
a=-k;

%%   预测
Result=zeros(N,1);
DataHat=zeros(n,1);
DataHat(1,1)=Data(1,1);
for i=2:n-1
    DataHat(i,1)=(1-exp(a))*(Data(1,1)-b/a)*exp(-a*(i-1));
end
for i=1:N
    Result(i,1)=(1-exp(a))*(Data(1,1)-b/a)*exp(-a*(n+i-1));
end

%%   计算残差及级比偏差
%   相对残差
RelativeResidual=(Data(2:end,1)-DataHat(2:end,1))./Data(2:end,1);
%   级比偏差
Temp=Data(2:end,1)./Data(1:end-1,1);
Eta=abs(1-(1-0.5*a)/(1+0.5*a)*(1./Temp));
end

新信息模型

function [Result]=New(Data,N)
Result=zeros(N,1);
for i=1:N
    Result(i,1)=Original(Data,1);
    Data=[Data;Result(i,1)];
end
end

此函数中用到了全数据模型中的求解步骤

新陈代谢模型

function [Result]=Metabolism(Data,N)
Length=length(Data);
Result=zeros(N,1);
for i=1:N
    Result(i,1)=Original(Data,1);
    for j=1:Length-1
        Data(j,1)=Data(j+1,1);
    end
    Data(Length,1)=Original(Data,1);
end
end

此函数中用到了全数据模型中的求解步骤