灰度模型GM(1,1)-学习笔记

忘了_

已于 2023-11-21 08:54:44 修改

阅读量566

点赞数

文章标签：学习笔记

于 2023-11-16 15:56:55 首次发布

本文链接：https://blog.csdn.net/weixin_42282506/article/details/134437932

版权

灰色模型 GM(1,1) Gray Model

Gray Model的本质：一种利用一阶微分方程拟合一维变量的数据实现预测的算法

正文：

0-写在前面

本文以GM(1,1)为例进行介绍，在GM(1,1)中第一个1代表方程的阶数，第二个1代表预测对象的数量。

GM(1,1)模型的本质是通过对原始数据做生成处理来寻找系统变动的规律，利用一阶微分方程拟合一维变量的数据实现预测。为什么用微分方程来拟合？是因为创始人受到了其他教授的影响，认为微分方程很好。详见知乎问答https://www.zhihu.com/question/345098724/answer/1252671607
在这里插入图片描述

白色系统：系统的信息是完全明确的。不仅知道系统的输入-输出关系，还知道实现输入-输出的具体方式。

黑色系统：系统的内部信息是完全未知的。只知道一些输入和相应的输出。

灰色系统：介于黑色和白色之间，系统的部分信息已知，部分信息未知。

后续算法过程及代码均为笔者本人参考网页、论文后梳理消化，自行写作，如有错误，恳请批评指正。

1-GM(1,1)的特点

1.1 能利用少部分数据进行预测，并且预测性能较好，简单常用；

1.2适用于原始离散的非负序列；

1.3由于数据量较少，预测结果易受到数据的影响。

2-GM(1,1)预测计算步骤

2.1原始序列级比检验

2.1.1为什么要对原始序列进行级比检验？

因为需要判断是原始数据是否可以直接用于GM（1,1）模型，如果不符合使用条件，需要对原始数据进行平移操作，以达到使用GM(1,1)模型的条件。

2.1.2如何对原始数列进行级比检验？

设原始数列为x(0)= 在这里插入图片描述
定义级比：

如果所有级别都满足区间

，则原始数据可以使用灰度模型进行预测，否则需要对原始数据进行平移变换在这里插入图片描述

，使其满足上述技术比要求。

2.2 建立GM模型

设原始数列为
在这里插入图片描述

2.2.1对原始数据进行累加，得到一次累加序列X⁽¹⁾_k

累加的目的是弱化随机序列的波动性和随机性。一次累加上标为1，二次累加上标为2，以此类推。原始数据上标为0说明累加的次数为0，即没做过累加。累加公式如公式（1）。
在这里插入图片描述
因此得到累加序列

2.2.2对累加序列进行平滑，得到数列Z⁽¹⁾

平滑的目的是原始数据中间的点是观测不到的，并且累加序列之间的点是不知道的，需要通过平滑在已知的各点之间增加点，平滑公式如（2）所示。

在这里插入图片描述
权值α可以在调试中多多尝试，哪个效果好就用那个。当α=0.5时，Z⁽¹⁾_k

为紧邻均值生成数列。

2.2.3 建立离散微分方程，求a和b

定义离散微分方程为

在这里插入图片描述
在方程（3）中，X⁽⁰⁾_k已知，是原始数据，有观测器观测得到；Z⁽¹⁾_k已知，是由一次累加序列经过平滑得到。a是发展系数，b是灰作用量，均未知。

利用矩阵的思想以及最小二乘法进行求解：

对公式（3）移项可得
在这里插入图片描述
设
则有X=BU，根据最小二乘法

则有在这里插入图片描述
X，即可得到估计的 a_hat和b_hat 的值。

2.2.4 建立连续微分方程，求一次累加序列

在这里插入图片描述
带入初值k=1，即
又因为k=1时，有

现在的得到的是一次累加数列的预测公式。

2.3逆用累加序列定义公式，求原始数据序列

上一个阶段的到的是一次累加数列的预测公式，若想得到是原始数据的预测公式，需要逆用累加公式来将一次累加数列还原成原始数据。

已知累加数列公式
在这里插入图片描述
当1≦k≦n时，所求值为拟合值，当k>n时，所求值为预测值。至此，用模型预测的阶段已经完成。

2.4 对模型进行检验

在模型构建好后，需要检验一下预测是否合理，只有通过检验的模型才可以用。常用的检验方法有很多，如相对残差检验法、方差比检验法、小误差概率检验法等。本文将对相对残差检验法展开介绍。

相对残差检验法的检验原理是看真实值和拟合值之间的差距大不大，如果大的话说明拟合效果不好，模型不行，所以k的取值是1~n，利用已知数据判断模型性能。

绝对残差的定义：
在这里插入图片描述当e<0.2时，模型预测效果一般；当e<0.1时，模型预测效果较好。

平均相对残差：
在这里插入图片描述
当e_bar<0.2时，模型预测效果一般；当e_ba<0.1时，模型预测效果较好。

2.5 残差模型修正

如果用原始序列经过变换所构建的GM（1,1）模型的预测效果不好，需要对所构建的GM（1,1）模型进行残差修正以提高模型的预测精度。

修正的方法是建立GM(1,1)的残差模型，也就是以GM(1,1)模型的残差序列为基础，再次建立GM(1,1)模型。

2.5.1构建残差模型的原始数据序列

设
在这里插入图片描述
其中，
为X⁽¹⁾的残差序列，即一次累加序列数据与它的估计值的差。如果存在k0，满足

在这里插入图片描述
则包含 e⁽⁰⁾_(k)及其之后序列的序列段可以用于GM(1,1)残差模型的建立。当k≦ n时，是对GM(1,1)拟合序列的修正；当k>n时，是对GM(1,1)预测序列的修正。

也就是说，只要存在一个满足条件的k0 ，那么在这里插入图片描述
就可以用来作为GM(1,1)残差模型的原始数据进行建模。

同时，称
在这里插入图片描述
为可建模残差尾端，仍记为

需要注意的是，如果残差序列不满足“符号一致”的条件，那么可以采用文献【基于ln(x+c)变换的GM(1,1)模型及其变形预测，邱利军、张波、张京奎】中的为加常数的方法进行序列变换后，建模预测，最后将预测结果反向还原即可。

2.5.2 对残差原始数据进行累加，得到残差数据的一次累加序列

生成可建模残差尾端的累加序列
在这里插入图片描述

2.5.3对累加序列进行平滑，得到数列Z⁽¹⁾。

在这里插入图片描述
一般情况下，α=0.5。

2.5.4构建离散微分方程，得到系数a_e和b_e

定义离散微分方程为
在这里插入图片描述
在上述方程中，e⁽⁰⁾_k已知，是原始数据，由GM（1,1）模型得到；Z⁽⁰⁾_k已知，是由残差的一次累加序列经过平滑得到。 a_e和b_e均为未知系数。

利用矩阵的思想以及最小二乘法进行求解：

对公式（3）移项可得
在这里插入图片描述

设
在这里插入图片描述
则有B=XU，根据最小二乘法

则有在这里插入图片描述
即可得到估计值 a_e_hat和b_e_hat。

2.5.5构建连续微分方程，得到残差数据一次累加序列的预测公式

在这里插入图片描述
现在的得到的是残差的一次累加数列的预测公式。

2.5.6通过残差系数、一次累加数列预测公式和残差一次累加数列预测公式的导数，共同构成修正后的一次累加数列预测公式

设残差系数由m表示，
在这里插入图片描述

2.5.7逆用累加序列定公式，求修正后的原始数据预测数据

上一个阶段的到的是残差的一次累加数列的预测公式，若想得到是残差原始数据的预测公式，需要逆用累加公式来将一次累加数列还原成原始数据。

已知累加数列公式
在这里插入图片描述
当-1≦k≦n时，所求值为残差的拟合值，当k>n时，所求值为残差预测值。至此，用残差模型预测的阶段已经完成。

2.5.8 最后用修正的点去替换误差过大的点即可

3-MATLAB程序

3.1 主函数main.m

clc
clear
close all

x0=[331.2692	331.1492	331.0294	330.9100	330.7908	330.6719	330.5532	330.4349	330.3168	330.1990];
% 
% %GM(1,1)
%%%1 [y1,y0,U]=Gary_model(x0);%使用Gary_model没发生误差过大的情况，因此又列了一组新的数据y0、y1来调试残差检验error_correction的代码
%运行时，%%%1和%%%2代码块二选一即可

%%%2
y0=[331.1048 330.1906 330.8765 330.7624 330.6483]
y1=[331.2692 660.3741 993.3647 1324.2413 1655.0037]
U=[0.000344890668966504 331.276205806469]
%%%

%残差检验
num=length(y0)

%判断是否有误差过大的情况出现
r=zeros(1,num);
e0=zeros(1,num)
%for i=2:num-1
for i=1:num-1
    e0(i+1)=y0(i+1)-y0(i);
    e0(i+1)=e0(i+1)*100;
    if e0(i+1)>20
       r(i)=1; %有误差过大的情况，标记r=1，对该点进行修正
    else
        continue
    end 
end

%求残差序列
for i=1:num-1
    e(i+1)=y1(i+1)-y1(i);
end

y2_e=zeros(1,num);
if sum(r)>0
    y2_e=error_correction(e,x0,U);%调用修正函数
end

y=y0;
for i=1:num
    if r(i)==1
        y(i)=y2_e(i);
    end
end

3.2 Gary_model.m

%GM(1,1)
function [y1,y0,a]=Gary_model(x0)
%得到一次累加序列
    num=length(x0);
    x1=zeros(1,num);
    x1(1)=x0(1);
    for i=2:num
        x1(i)=x1(i-1)+x0(i);
    end

    %得到均值平滑序列
    z1=zeros(1,num);
    for i=2:num
        z1(i)=0.5*x1(i)+0.5*x1(i-1);%z1下标从2开始
    end

    %调用函数-最小二乘法解a和b
    a=f(z1,x0);

    %调用函数-求一次累加序列预测序列
    y1=f1(a,x0);

    %还原原始数据的预测序列
    for i=1:num-1
        y0(i+1)=y1(i+1)-y1(i);
    end

end 


%最小二乘法解a和b
function [U]=f(z,x0)
    num=length(z);
    X=zeros(num,1);
    for i=1:num-1
        X(i+1)=[x0(i+1)]; %X下标从2开始
    end
    B=zeros(num,2);
    %X(1,:)=[z(1),1];
%     X(1,1)=z(1);
%     X(1,2)=1;
    for i=1:num-1
        B(i+1,:)=[-z(i+1) 1]; %B下标从2开始
    end
    U=inv((B')*B)*(B')*X;
end

%求一次累加序列的预测序列
function  [y1]=f1(U,x0)
    num=length(x0);
    for i=0:num-1
       temp=U(2)/U(1);
       y1(i+1)=temp+(x0(1)-temp)*exp(-U(1)*(i));
    end
end

3.3error_correction.m

function [y1_e]=error_correction(e,x0,U)
%输入e：残差序列
%输入x0：原始序列
%输出y2:修正后的原始数据预测序列

%首先选择合适的K0

%判断残差序列符号是否一致
k0=symbol_same_or_not(e,1);

%判断是否满足n-k0>=4这一条件
num=length(e);
if ~(num-k0>=4)
    k0=symbol_same_or_not(e,k0);
end


%找到了符合条件的k0
%求残差系数flag
for i=1:k0-1
    flag(i)=0;
end
for i=k0:num
    flag(i)=1;
end

%求残差序列的累加序列
    x2=zeros(1,num);
    x2(k0)=e(k0);
    for i=k0+1:num
        x2(i)=x2(i-1)+x2(i);
    end

    %得到均值平滑序列
    z2=zeros(1,num);
    for i=k0+1:num
        z2(i)=0.5*x2(i)+0.5*x2(i-1);%z1下标从k0+1开始
    end

    %调用函数-最小二乘法解a和b
    a1=f(k0,z2,e);

 %调用函数-求残差序列的修正预测序列
    y1_e=f1(k0,U,a1,e,x0,flag);
    

    %还原原始数据的预测序列
    for i=1:num-1
        y0_e(i+1)=y1_e(i+1)-y1_e(i);
    end

end



%判断残差序列符号是否一致
function [k0]=symbol_same_or_not(e,t) 
num=length(e);
x=zeros(1,num);
for i=t:num
    if e(i)>0
        x(i)=1;
    else
        x(i)=0;
    end  
end

for i=1:num
    %m=ones(1,num-i+1);
    %n=zeros(1,num-i+1);
    if isequal(x(1,i:end),ones(1,num-i+1))==1 || isequal(x(1,i:end),zeros(1,num-i+1))==1
       k0=i;
       break
       %如果满足if条件，说明从i开始往后的序列符号一样，那么就记录当前的i，并跳出循环
    else
         %如果不满足if条件，继续循环找到i
       continue
    end  
end
end


%最小二乘法解a和b
function [Ue]=f(k0,z,e)
    num=length(z);
    X=zeros(num,1);
    for i=k0:num-1
        X(i+1)=[e(i+1)]; %X下标从k0+1开始
    end
    B=zeros(num,2);
    %X(1,:)=[z(1),1];
%     X(1,1)=z(1);
%     X(1,2)=1;
    for i=k0:num-1
        B(i+1,:)=[-z(i+1) 1]; %X下标从2开始
    end
    Ue=inv((B')*B)*(B')*X;
end

%求残差序列的修正预测序列
function  [y1_e]=f1(k0,U,Ue,e,x0,flag)
    num=length(e);
    for i=k0-1:num
       temp=U(2)/U(1);
       temp_e=Ue(2)/Ue(1);
       y1_e(i+1)=(-Ue(1))*(e(k0)-temp_e)*exp(-Ue(1)*(i-k0))*(1-exp(Ue(1)))*flag(i)+ (x0(1)-temp)*exp(-U(1)*(i-1))*(1-exp(U(1)));
    end
end

忘了_

关注

0
点赞
踩
5

收藏

觉得还不错? 一键收藏
1
评论
灰度模型GM(1,1)-学习笔记

上一个阶段的到的是一次累加数列的预测公式，若想得到是原始数据的预测公式，需要逆用累加公式来将一次累加数列还原成原始数据。已知累加数列公式当1≦k≦n时，所求值为拟合值，当k>n时，所求值为预测值。至此，用模型预测的阶段已经完成。clcclearclose all% %GM(1,1)%使用Gary_model没发生误差过大的情况，因此又列了一组新的数据y0、y1来调试残差检验error_correction的代码%运行时，%%%1和%%%2代码块二选一即可%%%2%%%%残差检验。
复制链接

扫一扫