【MATLAB】GM(1,1) 灰色预测模型及算法

一、灰色预测模型概念

灰色预测是一种对含有不确定因素的系统进行预测的方法。 灰色预测通过鉴别系统因素之间发展趋势的相异程度，即进行关联分析，并对原始数据进行生成处理来寻找系统变动的规律，生成有较强规律性的数据序列，然后建立相应的微分方程模型，从而预测事物未来发展趋势的状况。 其用等时距观测到的反映预测对象特征的一系列数量值构造灰色预测模型，预测未来某一时刻的特征量，或达到某一特征量的时间。

二、使用步骤

1.使用条件

这个应该是很清晰的，数据量少，四五个，然后类指数形式增长的，什么是类指数增长：就是一组数据累加后排列是否呈类指数增长，与且预测的时间不是很长，就可以用。 分为GM (1,1)，GM (1,m)，GM (n,m)分别用于一个自变量一个因变量，多个自变量一个因变量，多个自变量多个因变量。 灰色预测就是尽可能使用数据中含有的信息。 假设你有十组数据，需要预测接下来的三组数据，程序跑完，一般是有求残差的过程，看一看是不是小于0.1，如果每个数据点都是小于0.1，那这次灰色预测就是很好的。

注：

数据量较少的情况下使用；
数据呈类指数增长；
预测时间较短
要先进行级比值检验
后验差比检验
模型拟合检验
模型残差检验

2.以GM（1,1）为例
（1）级比值检验

 

（注：级比值介于区间[0.982,1.0098]时说明数据适合模型构建。）

从上表可知，针对某数据进行GM(1,1)模型构建，结果显示：级比值的最大值为1.010，在适用范围区间[0.982,1.0098]之外，意味着本数据进行GM(1,1)可能得不到满意的模型。但从数据来看，1.01非常接近于1.0098，因此有理由接着进行建模。

（2）后验差比检验

后验差比C值用于模型精度等级检验，该值越小越好，一般C值小于0.35则模型精度等级好，C值小于0.5说明模型精度合格，C值小于0.65说明模型精度基本合格，如果C值大于0.65，则说明模型精度等级不合格。

三、总结

 算法（1）

clc,clear
a = [2.874,3.278,3.337,3.39,3.679];%初始数据
n = length(a); %数据的长度
k=1:n;
k2=1:(n+2);  %【】要预测未来几年的或其他的，改这里，再只改下面b4和b5的n,改成一样的，b5存的就是预测的值%其他都不需要改%(n+2)意为预测未来2年的  
%进行累加
b=zeros(1,n);  %对存放累加的矩阵进行初始化
b(1)=a(1);   %累加的第一项
for i=2:nb(i)=b(i-1)+a(i);
end%计算光滑比
b1=zeros(1,n-1); %初始化
for i=1:(n-1)b1(i)=a(i+1)/b(i);  %即b1(1) = a(2)/a(1)  b1(2) = a(3)/(a(1)+a(2)) ....
endc = 0.7;%光滑阈值
count = 0; %计数
for i=1:(n-1)if b1(i)<0.6   %这里可以变count=count+1;end
end
if(count/(n-1))>0.7 %cdisp('数据是光滑的')
end%计算模型背景值
b2=zeros(1,n-1);
for i=2:nb2(i-1) = 0.5*b(i-1)+0.5*b(i);
end%计算最小二乘估计
A=[(-b2)',ones(n-1,1)];
Y=a(2:n)'  %a(2:5)指的是取原始数据第2到5个，再进行转置
b3=zeros(2,1);
b3=inv(A'*A)*(A'*Y)  %最终得到两个数a和b%计算一次累加的预测值
b4=zeros(1,n+2);
b4(1)=b(1);
for i=2:n+2b4(i)=(b(1)-b3(2)/b3(1))*exp(-b3(1)*(i-1))+b3(2)/b3(1)
end%最终计算预测值
b5=zeros(1,n+2);
b5(1)=b4(1);
for i=2:n+2b5(i)=b4(i)-b4(i-1)
endscatter(k,a,'r') %画散点图    初始数据
hold on %可使两个图在一个窗口上显示
scatter(k2,b5,'b')    %预测值

算法（2）

clc;clear;      %建立符号变量a(发展系数)和b(灰作用量)
syms a b;
c = [a b]';    A = [1 4 6 9 10 12 21 34];    %输入需要预测的数据T1=length(A);
T2=100;                             %输入需要预测数据个数 
t1=1:T1;
t2=1:T1+T2;
n = T1;m = length(A);
JiBi = ones(1,m-1);
for i =2:mJiBi(i-1) = A(i-1)/A(i);
end
max1 = max(JiBi);
min1 = min(JiBi);
FanWei = exp(2/(n+2))-exp(-2/(n+1))
if max1 - min1<FanWeidisp(['数据通过级别检验']);
else disp(['数据不通过级比检验']);
end%对原始数列 A 做累加得到数列 B
B = cumsum(A);
%对数列 B 做紧邻均值生成
for i = 2:nC(i) = (B(i) + B(i - 1))/2; 
end
C(1) = [];%构造数据矩阵 
B = [-C;ones(1,n-1)];
Y = A; Y(1) = []; Y = Y';      %Y进行了转置，C的公式求法与百度文库 发生了一些变化
% 使用最小二乘法计算参数 a(发展系数)和b(灰作用量)
c = inv(B*B')*B*Y;          %核心公式  
c = c';
a = c(1); b = c(2);
disp(['发展系数：',num2str(a)]);
disp(['灰色作用量：',num2str(b)]);
%预测后续数据
F = []; F(1) = A(1);for i = 2:T1+T2F(i) = (A(1)-b/a)*exp(-a*(i-1))+ b/a;
end
%对数列 F 累减还原,得到预测出的数据
G = []; G(1) = A(1);
for i = 2:T1+T2G(i) = F(i) - F(i-1); %得到预测出来的数据
end
disp(['预测数据为：',num2str(G)]);%模型检验
H = G(1:T1);
epsilon = A - H;               %计算残差序列disp(['残差检验：',num2str(epsilon)]);
%法一：计算相对误差Q
delta = abs(epsilon./A);
Q = mean(delta);
disp(['相对残差Q检验：',num2str(Q)]);%法二：方差比C检验
C = std(epsilon, 1)/std(A, 1);  %方差函数std   按照列分
disp(['方差比C检验：',num2str(C)]);%法三：小误差概率P检验
S1 = std(A, 1);
tmp = find(abs(epsilon - mean(epsilon))< 0.6745 * S1);
P = length(tmp)/n;
disp(['小误差概率P检验：',num2str(P)])
%级比偏差和相对残差
check(A)%绘制曲线图
plot(t1, A,'ro'); hold on;
plot(t2, G, 'g-');
grid on;