「 数学模型 」“灰色预测GM(1,1)模型原理及Matlab实例”讲解

一、前言

        小白将使用的灰色预测模型的Matlab代码贴在下面，代码使用了灰色模型进行两次回归预测，缺点是精度不高，预测后的残差有点大。

        灰色预测模型主要内容：

1. GM(1,1)模型（1阶1个影响因素）；
2. 离散灰色模型（解决GM(1,1)中的理论缺陷）；
3. GM(1,N)模型（多个变量影响）；
4. GM(n,h)模型（多个变量高阶影响）；
5. 多变量离散灰色模型（解决GM(n,h)中的理论缺陷）；

二、模型原理

        1. GM(1,1)模型内容

        差分方程；

        取均值作中间状态，减少误差；——背景值

        微分方程——影子方程

根据已知数据，通过软件求出a,b参数后，需要通过关联度来判断预测方程的准确程度，如下例：

         2. 灰色预测模型扩展及应用

        （1）灰色预测模型扩展

        1）残差GM(1,1)模型

        数据选取太长时，会出现估计不准情况，所以引申出残差模型。

        之所以选择时间分段，很可能是因为数据选择过长，而这一长段数据不能用同一模型描述，出现模型分段现象。这时使用CM(1,1)残差模型可增加数据预测的准确性。

这也是灰度预测模型的难点之一：到底选取多少个数据是准确的。（通常选取4个或5个，但很难准确确定某个数）

        2）GM(1,1)模型群

        （2）GM(1,1)模型应用：数列预测、区间预测(比例带/包络带/发展带)、灾变预测、波形预测。

        1）区间预测——比例带：

        2）区间预测——包络带：

        3）区间预测——发展带：

        4）灾变预测

        灾变预测实质上是异常值预测，主要是给出下一个或几个异常值出现的时刻，以便人们提前准备，采取对策。

        上灾变序列和下灾变序列统称灾变序列。灾变预测就是通过对灾变序列的研究，寻找其规律性，预测以后若干次灾变发生的日期。实例：

        然后进行GM(1,1)建模，预测下一个灾害年份。

        5）波形预测

        当原始数据频频波动且摆动幅度较大时，往往难以找到合适的模拟模型，因此可以考虑根据原始数据的波形预测未来行为数据发展变化的波形。思想是：找出等高线。

        3. 灰色预测中的累加生成建模思想

        灰色预测模型通过数据处理来分析和对待随机量，即通过数据到数据的“映射”、时间序列到时间序列的“映射”来处理和发现规律。通过映射把原来存在扰动项的微分相互抵消掉。右边图示是处理后的结果，可以看出数据的明显变化趋势。

三、Matlab代码实例

%二次拟合预测GM(1,1)模型
%预测数据
 
x=[x1,x2,x3,x4];    %x[]为输入数组
sizexd2 = size(x,2);
%求数组长度
k=0;
for y1=x
    k=k+1;
        if k>1
        x1(k)=x1(k-1)+x(k);
         %累加生成
         z1(k-1)=-0.5*(x1(k)+x1(k-1));
         %z1维数减1，用于计算B
         yn1(k-1)=x(k);
        else
         x1(k)=x(k);
        end
end
%x1,z1,k,yn1
sizez1=size(z1,2);
%size(yn1);
z2 = z1';
z3 = ones(1,sizez1)';
YN = yn1';     %转置
%YN
B=[z2 z3];
au0=inv(B'*B)*B'*YN;
au = au0';
%B,au0,au
afor = au(1);
ufor = au(2);
ua = au(2)./au(1);
%afor,ufor,ua?
%输出预测的? a u 和 u/a的值
constant1 = x(1)-ua;
afor1 = -afor;
x1t1 = 'x1(t+1)';
estr = 'exp';
tstr = 't';
leftbra = '(';
rightbra = ')';
%constant1,afor1,x1t1,estr,tstr,leftbra,rightbra
strcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,rightbra,'+',leftbra,num2str(ua),rightbra)
%输出时间响应方程
%******************************************************
%二次拟合
k2 = 0;
for y2 = x1
     k2 = k2 + 1;
     if k2 > k
     else
         ze1(k2) = exp(-(k2-1)*afor);
     end
end
%ze1
sizeze1 = size(ze1,2);
z4 = ones(1,sizeze1)';
G=[ze1' z4];
X1 = x1';
au20=inv(G'*G)*G'*X1;
au2 = au20';
%z4,X1,G,au20
Aval = au2(1);
Bval = au2(2);
%Aval,Bval
%输出预测的? A,B的值
strcat(x1t1,'=',num2str(Aval),estr,leftbra,num2str(afor1),tstr,rightbra,'+',leftbra,num2str(Bval),rightbra)
%输出时间响应方程
nfinal = sizexd2-1 + 1;
%决定预测的步骤数5? 这个步骤可以通过函数传入
%nfinal = sizexd2 - 1 + 1;
%预测的步骤数 1
for k3=1:nfinal
         x3fcast(k3) = constant1*exp(afor1*k3)+ua;
end
%x3fcast
%一次拟合累加值
for  k31=nfinal:-1:0
     if k31>1
         x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1);
     else
             if k31>0
                 x31fcast(k31+1) = x3fcast(k31)-x(1);
             else
                 x31fcast(k31+1) = x(1);
             end
     end
    
end
x31fcast
%一次拟合预测值
 
for  k4=1:nfinal
     x4fcast(k4) = Aval*exp(afor1*k4)+Bval;
end
%x4fcast
for   k41=nfinal:-1:0
     if k41>1
         x41fcast(k41+1) = x4fcast(k41)-x4fcast(k41-1);
     else
         if k41>0
             x41fcast(k41+1) = x4fcast(k41)-x(1);
         else
             x41fcast(k41+1) = x(1);
         end
      end
 
end
x41fcast,x
%二次拟合预测值
%***精度检验p C************//
k5 = 0;
for y5 = x
     k5 = k5 + 1;
         if k5 > sizexd2
         else
         err1(k5) = x(k5) - x41fcast(k5);
         end
end
%err1
%绝对误差
 
xavg = mean(x);
%xavg
%x平均值
err1avg = mean(err1);
%err1avg
%err1平均值
k5 = 0;
s1total = 0 ;
for y5 = x
     k5 = k5 + 1;
     if k5 > sizexd2
     else
         s1total = s1total + (x(k5) - xavg)^2;
     end
end
s1suqare = s1total ./ sizexd2;
s1sqrt = sqrt(s1suqare);
%s1suqare,s1sqrt
%s1suqare? 残差数列x的方差? s1sqrt 为x方差的平方根S1
k5 = 0;
s2total = 0 ;
for y5 = x
     k5 = k5 + 1;
         if k5 > sizexd2
         else
         s2total = s2total + (err1(k5) - err1avg)^2;
         end
end
s2suqare = s2total ./ sizexd2;
%s2suqare?? 残差数列err1的方差S2
Cval = sqrt(s2suqare ./ s1suqare);
Cval
%nnn = 0.6745 * s1sqrt
%Cval? C检验值
k5 = 0;
pnum = 0 ;
for y5 = x
     k5 = k5 + 1;
     if abs( err1(k5) - err1avg ) < 0.6745 * s1sqrt
     pnum = pnum + 1;
         %ppp = abs( err1(k5) - err1avg )?????
     else
     end
end
pval = pnum ./ sizexd2;
pval
%p检验值
%arr1 = x41fcast(1:6)
%预测结果为区间范围? 预测步长和数据长度可调整程序参数进行改进

%二次拟合预测GM(1,1)模型
%预测数据
 
x=[4,5,6,7]; %写入的数据
sizexd2 = size(x,2);
%求数组长度
k=0;
for y1=x
    k=k+1;
        if k>1
        x1(k)=x1(k-1)+x(k);
         %累加生成
         z1(k-1)=-0.5*(x1(k)+x1(k-1));
         %z1维数减1，用于计算B
         yn1(k-1)=x(k);
        else
         x1(k)=x(k);
        end
end
%x1,z1,k,yn1
sizez1=size(z1,2);
%size(yn1);
z2 = z1';
z3 = ones(1,sizez1)';
YN = yn1';     %转置
%YN
B=[z2 z3];
au0=inv(B'*B)*B'*YN;
au = au0';
%B,au0,au
afor = au(1);
ufor = au(2);
ua = au(2)./au(1);
%afor,ufor,ua?
%输出预测的? a u 和 u/a的值
constant1 = x(1)-ua;
afor1 = -afor;
x1t1 = 'x1(t+1)';
estr = 'exp';
tstr = 't';
leftbra = '(';
rightbra = ')';
%constant1,afor1,x1t1,estr,tstr,leftbra,rightbra
strcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,rightbra,'+',leftbra,num2str(ua),rightbra)
%输出时间响应方程
%******************************************************
%二次拟合
k2 = 0;
for y2 = x1
     k2 = k2 + 1;
     if k2 > k
     else
         ze1(k2) = exp(-(k2-1)*afor);
     end
end
%ze1
sizeze1 = size(ze1,2);
z4 = ones(1,sizeze1)';
G=[ze1' z4];
X1 = x1';
au20=inv(G'*G)*G'*X1;
au2 = au20';
%z4,X1,G,au20
Aval = au2(1);
Bval = au2(2);
%Aval,Bval
%输出预测的? A,B的值
strcat(x1t1,'=',num2str(Aval),estr,leftbra,num2str(afor1),tstr,rightbra,'+',leftbra,num2str(Bval),rightbra)
%输出时间响应方程
nfinal = sizexd2-1 + 1;
%决定预测的步骤数5? 这个步骤可以通过函数传入
%nfinal = sizexd2 - 1 + 1;
%预测的步骤数 1
for k3=1:nfinal
         x3fcast(k3) = constant1*exp(afor1*k3)+ua;
end
%x3fcast
%一次拟合累加值
for  k31=nfinal:-1:0
     if k31>1
         x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1);
     else
             if k31>0
                 x31fcast(k31+1) = x3fcast(k31)-x(1);
             else
                 x31fcast(k31+1) = x(1);
             end
     end
    
end
x31fcast
%一次拟合预测值
 
for  k4=1:nfinal
     x4fcast(k4) = Aval*exp(afor1*k4)+Bval;
end
%x4fcast
for   k41=nfinal:-1:0
     if k41>1
         x41fcast(k41+1) = x4fcast(k41)-x4fcast(k41-1);
     else
         if k41>0
             x41fcast(k41+1) = x4fcast(k41)-x(1);
         else
             x41fcast(k41+1) = x(1);
         end
      end
 
end
x41fcast,x
%二次拟合预测值
%***精度检验p C************//
k5 = 0;
for y5 = x
     k5 = k5 + 1;
         if k5 > sizexd2
         else
         err1(k5) = x(k5) - x41fcast(k5);
         end
end
%err1
%绝对误差
 
xavg = mean(x);
%xavg
%x平均值
err1avg = mean(err1);
%err1avg
%err1平均值
k5 = 0;
s1total = 0 ;
for y5 = x
     k5 = k5 + 1;
     if k5 > sizexd2
     else
         s1total = s1total + (x(k5) - xavg)^2;
     end
end
s1suqare = s1total ./ sizexd2;
s1sqrt = sqrt(s1suqare);
%s1suqare,s1sqrt
%s1suqare? 残差数列x的方差? s1sqrt 为x方差的平方根S1
k5 = 0;
s2total = 0 ;
for y5 = x
     k5 = k5 + 1;
         if k5 > sizexd2
         else
         s2total = s2total + (err1(k5) - err1avg)^2;
         end
end
s2suqare = s2total ./ sizexd2;
%s2suqare?? 残差数列err1的方差S2
Cval = sqrt(s2suqare ./ s1suqare);
Cval
%nnn = 0.6745 * s1sqrt
%Cval? C检验值
k5 = 0;
pnum = 0 ;
for y5 = x
     k5 = k5 + 1;
     if abs( err1(k5) - err1avg ) < 0.6745 * s1sqrt
     pnum = pnum + 1;
         %ppp = abs( err1(k5) - err1avg )?????
     else
     end
end
pval = pnum ./ sizexd2;
pval
%p检验值
%arr1 = x41fcast(1:6)
%预测步长和数据长度可调整程序参数进行改进

最后

以上就是拉长板凳为你收集整理的「 数学模型 」“灰色预测GM(1,1)模型原理及Matlab实例”讲解的全部内容，希望文章能够帮你解决「 数学模型 」“灰色预测GM(1,1)模型原理及Matlab实例”讲解所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。