GM（1n） matlab代码

%灰色预测模型GM（1，n）模型的matlab源代码，包括预测模型的建立，以及模型的精度检验指标cp的计算
%假设预测3步，N=3
%如在命令窗口键入：
%gm=huise_gm1n（[201.9	219	221.1	324.5	346.97	332.9	359.53	406.98  440.69	475.96	479.25	482][30.68	42.86	47.44	61.87	79.5	115.7  156.6	193.5	218.58	237.84	272.96	298.65][13.5	13.8	14.1	15.1	16.01	16.1	16.1	16.1	16.2	16.3	16.4	16.5]5）;
function GM=huise_gm1n（data1data2data3N） %data1:垃圾，data2:环保，data3:收入
data1=[201.9	219	221.1	324.5	346.97	332.9	359.53	406.98	440.69	475.96	479.25	482];
data2=[30.68	42.86	47.44	61.87	79.5	115.7	156.6	193.5	218.58	237.84	272.96	298.65];
%data3=[20905.68	22759.92	24940.68	25935.84	27596.4	21494.4	22567.08	24301.38	26729.31	29244.52	32380	36505];
data3=[13.5	13.8	14.1	15.1	16.01	16.1	16.1	16.1	16.2	16.3	16.4	16.5];
N=5;
T=length（data1）;  
PYX1=data1;
PYX2=data2;
PYX3=data3;
%进行数据预处理，这里用初值化
X0_1=PYX1./PYX1（1）;
X0_2=PYX2./PYX2（1）;
X0_3=PYX3./PYX3（1）;
%用AGO生成一阶累加生成模块
X1_1（1）=X0_1（1）;
X1_2（1）=X0_2（1）;
X1_3（1）=X0_3（1）;
for i=2:T  
   X1_1（i）=X1_1（i-1）+X0_1（i）; 
   X1_2（i）=X1_2（i-1）+X0_2（i）;
   X1_3（i）=X1_3（i-1）+X0_3（i）;
end 
%构造累加矩阵B
for i=1:T-1 
   M1（i）=（0.5*（X1_1（i）+X1_1（i+1）））;
   M2（i）=（0.5*（X1_2（i）+X1_2（i+1）））;
   M3（i）=（0.5*（X1_3（i）+X1_3（i+1）））;
end 
 
B1=zeros（T-13）; 
for i=1:（T-1） 
    B1（i1）=-M1（i）;   %-（X1_1（i）+X1_1（i+1）））/2; 
    B1（i2）=X1_2（i+1）; 
    B1（i3）=X1_3（i+1）;
end
B2=zeros（T-12）; 
for i=1:（T-1） 
    B2（i1）=-M2（i）;   %-（X1_2（i）+X1_2（i+1）））/2; 
    B2（i2）=X1_3（i+1）; 
end
B3=zeros（T-12）; 
for i=1:（T-1） 
    B3（i1）=-M3（i）;   %-（X1_3（i）+X1_3（i+1）））/2; 
    B3（i2）=1; 
end
save B1 B1;
save B2 B2;
save B3 B3;
%构造常数项向量Y
for i=2:T                          
    Y1（i-1）=X0_1（i）; 
    Y2（i-1）=X0_2（i）;
    Y3（i-1）=X0_3（i）;
end 
HCS1=inv（B1‘*B1）*B1‘*Y1‘;               %用最小二乘法求灰参数HCS1 
H1=HCS1‘;                            %H1=[ab2b3]——即最终方程中系数
HCS2=inv（B2‘*B2）*B2‘*Y2‘;               %用最小二乘法求灰参数HCS2 
H2=HCS2‘;                            %H2=[ab3]  ?b2
HCS3=inv（B3‘*B3）*B3‘*Y3‘;               %用最小二乘法求灰参数HCS3 
H3=HCS3‘;                            %H3=[ba]   ?[ab]——用来预测Y3的，类比于GM（11）

%计算出X3的累加序列
for i=1:T+N                         
YCX13（i）=（X0_3（1）-H3（2）/H3（1））*exp（-1*H3（1）*（i-1））+H3（2）/H3（1）; 
en