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

一、前言

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

灰色预测模型主要内容：

GM(1,1)模型（1阶1个影响因素）；
离散灰色模型（解决GM(1,1)中的理论缺陷）；
GM(1,N)模型（多个变量影响）；
GM(n,h)模型（多个变量高阶影响）；
多变量离散灰色模型（解决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=xk=k+1;if k>1x1(k)=x1(k-1)+x(k);%累加生成z1(k-1)=-0.5*(x1(k)+x1(k-1));%z1维数减1，用于计算Byn1(k-1)=x(k);elsex1(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 = x1k2 = k2 + 1;if k2 > kelseze1(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:nfinalx3fcast(k3) = constant1*exp(afor1*k3)+ua;
end
%x3fcast
%一次拟合累加值
for  k31=nfinal:-1:0if k31>1x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1);elseif k31>0x31fcast(k31+1) = x3fcast(k31)-x(1);elsex31fcast(k31+1) = x(1);endendend
x31fcast
%一次拟合预测值for  k4=1:nfinalx4fcast(k4) = Aval*exp(afor1*k4)+Bval;
end
%x4fcast
for   k41=nfinal:-1:0if k41>1x41fcast(k41+1) = x4fcast(k41)-x4fcast(k41-1);elseif k41>0x41fcast(k41+1) = x4fcast(k41)-x(1);elsex41fcast(k41+1) = x(1);endendend
x41fcast,x
%二次拟合预测值
%***精度检验p C************//
k5 = 0;
for y5 = xk5 = k5 + 1;if k5 > sizexd2elseerr1(k5) = x(k5) - x41fcast(k5);end
end
%err1
%绝对误差xavg = mean(x);
%xavg
%x平均值
err1avg = mean(err1);
%err1avg
%err1平均值
k5 = 0;
s1total = 0 ;
for y5 = xk5 = k5 + 1;if k5 > sizexd2elses1total = 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 = xk5 = k5 + 1;if k5 > sizexd2elses2total = 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 = xk5 = k5 + 1;if abs( err1(k5) - err1avg ) < 0.6745 * s1sqrtpnum = pnum + 1;%ppp = abs( err1(k5) - err1avg )?????elseend
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=xk=k+1;if k>1x1(k)=x1(k-1)+x(k);%累加生成z1(k-1)=-0.5*(x1(k)+x1(k-1));%z1维数减1，用于计算Byn1(k-1)=x(k);elsex1(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 = x1k2 = k2 + 1;if k2 > kelseze1(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:nfinalx3fcast(k3) = constant1*exp(afor1*k3)+ua;
end
%x3fcast
%一次拟合累加值
for  k31=nfinal:-1:0if k31>1x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1);elseif k31>0x31fcast(k31+1) = x3fcast(k31)-x(1);elsex31fcast(k31+1) = x(1);endendend
x31fcast
%一次拟合预测值for  k4=1:nfinalx4fcast(k4) = Aval*exp(afor1*k4)+Bval;
end
%x4fcast
for   k41=nfinal:-1:0if k41>1x41fcast(k41+1) = x4fcast(k41)-x4fcast(k41-1);elseif k41>0x41fcast(k41+1) = x4fcast(k41)-x(1);elsex41fcast(k41+1) = x(1);endendend
x41fcast,x
%二次拟合预测值
%***精度检验p C************//
k5 = 0;
for y5 = xk5 = k5 + 1;if k5 > sizexd2elseerr1(k5) = x(k5) - x41fcast(k5);end
end
%err1
%绝对误差xavg = mean(x);
%xavg
%x平均值
err1avg = mean(err1);
%err1avg
%err1平均值
k5 = 0;
s1total = 0 ;
for y5 = xk5 = k5 + 1;if k5 > sizexd2elses1total = 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 = xk5 = k5 + 1;if k5 > sizexd2elses2total = 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 = xk5 = k5 + 1;if abs( err1(k5) - err1avg ) < 0.6745 * s1sqrtpnum = pnum + 1;%ppp = abs( err1(k5) - err1avg )?????elseend
end
pval = pnum ./ sizexd2;
pval
%p检验值
%arr1 = x41fcast(1:6)
%预测步长和数据长度可调整程序参数进行改进