灰度模型与灰色预测

概念

灰色系统是相对于黑色系统和白色系统而言的。

白色系统：系统的内部特征是完全已知的，即系统的信息是完全充分的。

黑色系统：一个系统的内部信息对外界来说是一无所知的，只能通过他与外界的联系来加以观测研究。

灰色系统：一部分信息是已知的，另一部分信息是未知的，系统内各因素之间具有不确定关系。其特点是‘少数据建模’，着重研究‘外延明确，内涵不明确’的对象。

灰色系统具有相对性与广泛性。指系统对于不同对象的灰度不一样。作为实际问题，灰色系统在大千世界中是大量存在的，绝对的白色或黑色系统是很少的。

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

基本思想

基本思想是用原始数据组成原始序列(0)，经累加生成法生成序列(1)，它可以弱化原始数据的随机性，使其呈现出较为明显的特征规律。对生成变换后的序列(1)
建立微分方程型的模型即GM模型。

GM(1,1) 模型表示1阶的、1个变量的微分方程模型。GM(1,1) 模型群中，新陈代谢模型是最理想的模型。这是因为任何一个灰色系统在发展过程中，随着时间的推移，将会不断地有一些随即扰动和驱动因素进入系统，使系统的发展相继地受其影响。用GM(1,1) 模型进行预测，精度较高的仅仅是原点数据(0)(n) 以后的1到2个数据，即预测时刻越远预测的意义越弱。而新陈代谢GM(1,1)模型的基本思想为越接近的数据，对未来的影响越大。也就是说，在不断补充新信息的同时，去掉意义不大的老信息，这样的建模序列更能动态地反映系统最新的特征，这实际上是一种动态预测模型。

优点

1、不需要大量样本。
2、样本不需要有规律性分布。
3、计算工作量小。
4、定量分析结果与定性分析结果不会不一致。
5、可用于短期、中长期预测。
6、灰色预测准确度高。

Codes

Aiq灰度预测模型：
clc
clear
clf
x0=[79,72.9550561800000,92.7802197800000,110.482352900000,82.2777777800000,66.9021739100000,81.2333333300000,103.168539300000,91.8750000000000,84.9120879100000,106.511111100000,88.6666666700000,83.4022988500000,64.7826087000000,76.8426966300000,92.6976744200000,73.5384615400000,71.2247191000000,84.6956521700000,77.4772727300000,89.0444444400000,63.4666666700000,88.9890109900000,90.7159090900000,76.8351648400000,68.6153846200000,89.5326087000000,81.6444444400000,80.6263736300000,68.9560439600000,84.6043956000000,77.1098901100000,71.6966292100000,63.1978022000000,85.1304347800000,82.8651685400000,70.8571428600000,62.5000000000000,84.4065934100000,85.5280898900000,80.3222222200000,71.9890109900000,87.9673913000000,79.8222222200000,78.8222222200000,68.1739130400000,81.3296703300000,85.3111111100000,75.4888888900000,66.3225806500000,83.5333333300000,226.384615400000,206.637893375000,186.891171350000,167.144449325000,147.397727300000,106.879120900000,88.8409090900000,146.113636400000,191.745454500000]
n=length(x0);%求x0的长度
x1=zeros(1,n);%生成与x0等长度的零向量
x1(1)=x0(1);%将x0首位赋值给x1
for i=2:n        %计算累加序列 x1    
    x1(i)=x1(i-1)+x0(i); %x1为x0的累加序列
end
i=2:n;     %对原始数列平行移位并赋给y 
y(i-1)=x0(i);
y=y' %将 y 变成列向量 
i=1:n-1;   %计算数据矩阵 B 的第一列数据
c(i)=-0.5*(x1(i)+x1(i+1));  %两个数平均值取负数，总感觉饶了一大圈
B=[c' ones(n-1,1)];%构造矩阵 B ，转置后找个1伴随
au=inv(B'*B)*B'*y;%计算参数 au 矩阵 
i=1:n+1;       %计算预测累加数列的值 
ago(i)=(x0(1)-au(2)/au(1))*exp(-au(1)*(i-1))+ au(2)/au(1); 
yc(1)=ago(1);
i=1:n-1;         %还原数列的值 
yc(i+1)=ago(i+1)-ago(i); 
i=2:n;
error(i)=yc(i)-x0(i);  %计算残差值 yc(1)=ago(1);
i=1:n-1;           %修正的还原数列的值 ，我怎么感觉没修正的样子
yc(i+1)=ago(i+1)-ago(i);
c=std(error)/std(x0);  %计算后验差比，也就是残差标准差和原数值标准差的比值
p=0; 
for i=2:n  
    if(abs(error(i)-mean(error))<0.6745*std(x0))      
        p=p+1; 
    end
end%看残差中有几个数值和残差均值的相差不大
p=p/(n-1);%p分配到每个参与的残差上，即为小误差概率的值
w1=min(abs(error)); 
w2=max(abs(error)); %计算残差的最大最小位
i=1:n;                %计算关联度 w
w(i)=(w1+0.5*w2)./(abs(error(i))+0.5*w2);
w=sum(w)/(n-1);
au             %输出参数 a,u 的值
x0            %输出原始序列值
ago   %输出累加数列 ago 的值

plot(1:60,x0,'+',1:60,yc,'*'); 
hold on;
plot(1:60,x0,1:60,yc); 
hold on;
yc              %输出预测的值 
grid on; 
error            %输出残差的值 
xlabel('季度');
c               %输出后验差比的值 
ylabel(' AQI');
p               %输出小误差概率的值
title(' 郑州空气质量指数灰色模型预测拟和曲线');
w               %输出关联度 w
legend(' 实测值','预测值',4); 


i=1:104;       %预测2015年至2025年
p_ago(i)=(x0(1)-au(2)/au(1))*exp(-au(1)*(i-1))+ au(2)/au(1); 
p_yc(1)=p_ago(1);
i=1:102;         %预测
p_yc(i+1)=p_ago(i+1)-p_ago(i); 
plot(p_yc)

Aqi随季度变化曲线图：
clc
clear
clf
y_season=[79 72.95505618 92.78021978 110.4823529 82.27777778 66.90217391 81.23333333 103.1685393...
91.875 84.91208791 106.5111111 88.66666667 83.40229885 64.7826087 76.84269663 92.69767442 73.53846154...
71.2247191 84.69565217 77.47727273 89.04444444 63.46666667 88.98901099 90.71590909 76.83516484 68.61538462 89.5326087...
81.64444444 80.62637363 68.95604396 84.6043956 77.10989011 71.69662921 63.1978022 85.13043478...
82.86516854 70.85714286 62.5 84.40659341 85.52808989 80.32222222 71.98901099 87.9673913...
79.82222222 78.82222222 68.17391304 81.32967033 85.31111111 75.48888889 66.32258065 83.53333333...
226.3846154...
147.3977273 106.8791209 88.84090909 146.1136364 191.7454545]%2010年第二季度到2015年第一季度的aqi的数值
x=[1:52 56:60]
xi=[1:60]
yi_season=interp1(x,y_season,xi,'linear')%采用线性插值方式，插值补充缺失的点
plot(xi,yi_season,'-mo','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor',[.49 1 .63], 'MarkerSize',6,'Color',[1 0 1])
box off;
set(gca,'XTick',[])%修改横坐标标志
text(1,48,'\fontname{30}2000年'),text(9,48,'\fontname{30}2002年'),text(17,48,'\fontname{30}2004年'),text(25,48,'\fontname{30}2006年'),text(33,48,'\fontname{30}2008年'),text(41,48,'\fontname{30}2010年')...
    ,text(49,48,'\fontname{30}2012年'),text(57,48,'\fontname{30}2014年')
title('2010年第二季度至2015年第一季度郑州市平均AQI变化曲线','FontWeight','bold','FontSize',20)
axis tight
ylabel('AQI指数')

灰色预测：
clc
clear
clf
x0=[20.0 54.7 73.3 84.0 89.3 100 106.7 121.1 133.3 143.3]
gm(x0)
子代码：
function [e,c,w,yc]=gm(x0)   %定义函数 gm(x0)
n=length(x0);%求x0的长度
x1=zeros(1,n);%生成与x0等长度的零向量
x1(1)=x0(1);%将x0首位赋值给x1
for i=2:n        %计算累加序列 x1    
    x1(i)=x1(i-1)+x0(i); %x1为x0的累加序列
end
i=2:n;     %对原始数列平行移位并赋给y 
y(i-1)=x0(i);
y=y' %将 y 变成列向量 
i=1:n-1;   %计算数据矩阵 B 的第一列数据
c(i)=-0.5*(x1(i)+x1(i+1));  %两个数平均值取负数，总感觉饶了一大圈
B=[c' ones(n-1,1)];%构造矩阵 B ，转置后找个1伴随
au=inv(B'*B)*B'*y;%计算参数 au 矩阵 
i=1:n+1;       %计算预测累加数列的值 
ago(i)=(x0(1)-au(2)/au(1))*exp(-au(1)*(i-1))+ au(2)/au(1); 
yc(1)=ago(1);
i=1:n-1;         %还原数列的值 
yc(i+1)=ago(i+1)-ago(i); 
i=2:n;
error(i)=yc(i)-x0(i);  %计算残差值 yc(1)=ago(1);
i=1:n-1;           %修正的还原数列的值 ，我怎么感觉没修正的样子
yc(i+1)=ago(i+1)-ago(i);
c=std(error)/std(x0);  %计算后验差比，也就是残差标准差和原数值标准差的比值
p=0; 
for i=2:n  
    if(abs(error(i)-mean(error))<0.6745*std(x0))      
        p=p+1; 
    end
end%看残差中有几个数值和残差均值的相差不大
p=p/(n-1);%p分配到每个参与的残差上，即为小误差概率的值
w1=min(abs(error)); 
w2=max(abs(error)); %计算残差的最大最小位
i=1:n;                %计算关联度 w
w(i)=(w1+0.5*w2)./(abs(error(i))+0.5*w2);
w=sum(w)/(n-1);
v=[1979,1988,20,160];  
au             %输出参数 a,u 的值
axis(v);
x0            %输出原始序列值
ago   %输出累加数列 ago 的值

plot([1979:n+1978],x0,'+',[1979:n+1978],yc,'*'); 
hold on;
plot([1979:n+1978],x0,[1979:n+1978],yc); 
yc              %输出预测的值 
grid on; 
error            %输出残差的值 
xlabel(' 时序');
c               %输出后验差比的值 
ylabel(' 沉降量(mm)');
p               %输出小误差概率的值
title(' 地面沉降灰色模型预测拟和曲线');
w               %输出关联度 w
legend(' 实测值','预测值',4); 

灰度预测：
x = [5999,5903,5848,5700,7884];
gm1(x);
子代码：
%二次拟合预测GM(1,1)模型
function  gmcal=gm1(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)


%预测结果为区间范围  预测步长和数据长度可调整程序参数进行改进

灰度预测完美模板：
clc
clear
clf
x0=[79,72.9550561800000,92.7802197800000,110.482352900000,82.2777777800000,66.9021739100000,81.2333333300000,103.168539300000,91.8750000000000,84.9120879100000,106.511111100000,88.6666666700000,83.4022988500000,64.7826087000000,76.8426966300000,92.6976744200000,73.5384615400000,71.2247191000000,84.6956521700000,77.4772727300000,89.0444444400000,63.4666666700000,88.9890109900000,90.7159090900000,76.8351648400000,68.6153846200000,89.5326087000000,81.6444444400000,80.6263736300000,68.9560439600000,84.6043956000000,77.1098901100000,71.6966292100000,63.1978022000000,85.1304347800000,82.8651685400000,70.8571428600000,62.5000000000000,84.4065934100000,85.5280898900000,80.3222222200000,71.9890109900000,87.9673913000000,79.8222222200000,78.8222222200000,68.1739130400000,81.3296703300000,85.3111111100000,75.4888888900000,66.3225806500000,83.5333333300000,226.384615400000,206.637893375000,186.891171350000,167.144449325000,147.397727300000,106.879120900000,88.8409090900000,146.113636400000,191.745454500000]
m=10;
n=length(x0);%求x0的长度
x1=zeros(1,n);%生成与x0等长度的零向量
x1(1)=x0(1);%将x0首位赋值给x1
for i=2:n        %计算累加序列 x1    
    x1(i)=x1(i-1)+x0(i); %x1为x0的累加序列
end
i=2:n;     %对原始数列平行移位并赋给y 
y(i-1)=x0(i);
y=y' %将 y 变成列向量 
i=1:n-1;   %计算数据矩阵 B 的第一列数据
c(i)=-0.5*(x1(i)+x1(i+1));  %两个数平均值取负数，总感觉饶了一大圈
B=[c' ones(n-1,1)];%构造矩阵 B ，转置后找个1伴随
au=inv(B'*B)*B'*y;%计算参数 au 矩阵 
i=1:n+1+m;       %计算预测累加数列的值 
ago(i)=(x0(1)-au(2)/au(1))*exp(-au(1)*(i-1))+ au(2)/au(1); 
yc(1)=ago(1);
i=1:n-1;         %还原数列的值 
yc(i+1)=ago(i+1)-ago(i); 
i=2:n;
error(i)=yc(i)-x0(i);  %计算残差值 yc(1)=ago(1);
i=1:n-1+m;           %修正的还原数列的值 ，我怎么感觉没修正的样子
yc(i+1)=ago(i+1)-ago(i);
c=std(error)/std(x0);  %计算后验差比，也就是残差标准差和原数值标准差的比值
p=0; 
for i=2:n  
    if(abs(error(i)-mean(error))<0.6745*std(x0))      
        p=p+1; 
    end
end%看残差中有几个数值和残差均值的相差不大
p=p/(n-1);%p分配到每个参与的残差上，即为小误差概率的值
w1=min(abs(error)); 
w2=max(abs(error)); %计算残差的最大最小位
i=1:n;                %计算关联度 w
w(i)=(w1+0.5*w2)./(abs(error(i))+0.5*w2);
w=sum(w)/(n-1);

au             %输出参数 a,u 的值
x0            %输出原始序列值
ago   %输出累加数列 ago 的值          
yc              %输出预测的值 
error            %输出残差的值
c     %输出后验差比的值
p    %输出小误差概率的值
w    %输出关联度 w



plot(1:n,x0,'-.r+',1:m+n,yc,'--k*'); 
xlabel('季度');
ylabel(' AQI');
title(' 郑州空气质量指数灰色模型预测拟和曲线');              
legend(' 实测值','预测值',4); 

一个应用：看病综合指数m代码
clc
clear
x0=[3.374145257 4.803207169 6.337084736 7.812989934  9.241120107 10.31063977 11.15319239 11.94645246 13.26184711]
n=length(x0);
x1=zeros(1,n);
x1(1)=x0(1);
for i=2:n      
    x1(i)=x1(i-1)+x0(i); 
end
i=2:n;   
y(i-1)=x0(i);
y=y'  
i=1:n-1;   
c(i)=-0.5*(x1(i)+x1(i+1));  
B=[c' ones(n-1,1)];
au=inv(B'*B)*B'*y;
i=1:n+1;       
ago(i)=(x0(1)-au(2)/au(1))*exp(-au(1)*(i-1))+ au(2)/au(1); 
yc(1)=ago(1);
i=1:n-1;         
yc(i+1)=ago(i+1)-ago(i); 
i=2:n;
error(i)=yc(i)-x0(i);  
i=1:n-1;          
yc(i+1)=ago(i+1)-ago(i);
c=std(error)/std(x0);  
p=0; 
for i=2:n  
    if(abs(error(i)-mean(error))<0.6745*std(x0))      
        p=p+1; 
    end
end
p=p/(n-1);
w1=min(abs(error)); 
w2=max(abs(error)); 
i=1:n;                
w(i)=(w1+0.5*w2)./(abs(error(i))+0.5*w2);
w=sum(w)/(n-1);
au            
x0           
ago   %

plot(1:n,x0,'+',1:n,yc,'*'); 
hold on;
plot(1:n,x0,1:n,yc); 
hold on;
yc              
grid on; 
error            
xlabel('年份');
c                
ylabel(' 看病难易综合指标');
p               
title(' 中国看病难易程度综合指标灰色模型预测拟和曲线');
w               
legend(' 实测值','预测值',4); 


i=1:20;       
p_ago(i)=(x0(1)-au(2)/au(1))*exp(-au(1)*(i-1))+ au(2)/au(1); 
p_yc(1)=p_ago(1);
i=1:18;        
p_yc(i+1)=p_ago(i+1)-p_ago(i); 
plot(p_yc)