%灰色預測模型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