1、 灰色预测模型 GM(1,1)的 matlab 运行代码例 由 19902001 年中国蔬菜产量，建立模型预测 2002 年中国蔬菜产量，并对预测结果作检验。分析建模：给定原始时间 19902001 年资料序列 X (k),对 X (k)0()0(生成 1-AGO(累加)序列 X (k)及 Y 。见下表)1(nK 1 2 3 4 5 6 7 8 9 10 11 12X 19519 ,19578 ,19637,19695,16602, 25723, 30379, 34473, 38485, 40514, 42400, 48337)0(X 19519, 39097, 58734, 264605,3
2、07005,355342Y - 19578 19637 40514 42400 48337n其中X (k)= ; Y =)1( )(1i0ikn TXX)12(,)3(),2(000对上述 X (k)的 GM(1,1),得到)0( 15236.-748.0125-1639.4-10782.563-.14895-230.1215.009.8507.65.0143.25011.1)2()0(19)8(71)6(5)4(132)1()()1()1()1()1()1()1()1()()1()(11111)(1 ）（）（ ）（）（ ）（）（ ）（）（ ）（）（ ）（）（ ）（）（ ）（）（ ）（）（ ）
3、（）（ ）（）（ XXXzzzzB将 B 和 Y 代入辨识算式，有：n1.605()139.TTnaBYb得灰色 GM(1,1)模型为（1）灰微分方程 X (k)-0.1062105 Z (k)=13999.9)0( )1(2)白化方程 9.31625.)1()1( Xdt(3)白化方程的时间响应式 5.1385.12)()1( 10625.0) tat ebeXtX(4)还原为原始数据预测方程： ,即)()()(1)0 tXtttet 1625.)0(15248.9(5)残差检验：残差 error1=e1= ,这里残差有 12 个。)()00( iXi（相对残差 error2=e2= ,这里
4、相对残差有 12 个。)0()0()()( 1| Xei(6) 后验差检验：C= , 其中 , S1 为绝对误差序列的标准差。21S1)(122)0(0niSi）（, )()()(000 eiXii)(210)0(iiS2 为原始数据系列标准差， , )(12)0(02nXiSi）（ )(120)0(iXiC0.6 不合格。利用 matlab 做求解 a,b,B,并作残差分析 x0=19519,19578,19637,19695,16602,25723,30379,34473,38485,40514,42400,48337; format long; (表示设计精度) n=length(x0)
5、; （输入数据长度） x1= ; （表示 x1 是一矩阵） x1(1)=x0(1); for i=2:n;x1(i)=x1(i-1)+x0(i);end for i=1:n-1;B(i,1)=-0.5*(x1(i)+x1(i+1); （矩阵 B 的第一列）B(i,2)=1; （矩阵 B 的第二列）Y(i)=x0(i+1); （表示 Yn 数据）end alpha=(B*B)(-1)*B*Y; a=alpha(1,1); b=alpha(2,1); d=b/a; （计算时间响应函数参数） c=x1(1)-d; x2(1)=x0(1); x(1)=x0(1); for i=1:n-1;x2(i+1
6、)=c*exp(-a*i)+d; （这里 x2(i+1)相当上面所讲的 ）)1()tXx(i+1)=x2(i+1)-x2(i); （这里 x(i+1) 相当原来输入数据的预测数据 ）)()0tend for i=2: 12;x2(i)=c*exp(-a*(i-1)+d; （对上面刚引出的 x2(i)进行说明及计算）x(i)=x2(i)-x2(i-1);end for i=1:n;error(i)=x(i)-x0(i); （残差）error1(i)=abs(error(i);（计算残差，abs 表示绝对值）error2(i)=error1(i)/x0(i); （计算相对误差）end C=std(
7、error1)/std(x0); （计算后验差检验数,std 表示标准差） k=1; （k 表示预测长度，这里每次预测下一年）a =-0.106210475032772 bb =1.399996741173038e+04BB =1.0e+05 *-0.293080000000000 0.000010000000000-0.489155000000000 0.000010000000000-0.685815000000000 0.000010000000000-0.867300000000000 0.000010000000000-1.078925000000000 0.000010000000
8、000-1.359435000000000 0.000010000000000-1.683695000000000 0.000010000000000-2.048485000000000 0.000010000000000-2.443480000000000 0.000010000000000-2.858050000000000 0.000010000000000-3.311735000000000 0.000010000000000 C （求后检验数）C =0.163969348419772 x （原始数据 的对应的预测数据 ，这里也是 12 个）)(0iX)(0iXx =1.0e+04 *
9、Columns 1 through 31.951900000000000 1.695769385830782 1.885790370699694Columns 4 through 62.097104330304606 2.332097268346191 2.593422554345592Columns 7 through 92.884030883565190 3.207203594115656 3.566589717441855Columns 10 through 123.966247180534730 4.410688625094428 4.904932361001964 eroor1 （求
10、残差）eroor1 =1.0e+03 *Columns 1 through 40 2.620306141692185 0.779096293003065 1.276043303046055Columns 5 through 86.718972683461907 0.211225543455919 1.538691164348100 2.400964058843441Columns 9 through 122.819102825581445 0.851528194652696 1.706886250944284 0.712323610019637 eroor2 （求相对误差）eroor2 =Co
11、lumns 1 through 40 0.133839316666267 0.039674914345525 0.064790215945471Columns 5 through 80.404708630494031 0.008211543888968 0.050649829301429 0.069647667996503Columns 9 through 120.073251989751369 0.021018121998635 0.040256751201516 0.014736611912606a =-0.106210475032772b =1.399996741173038e+04 B
12、B =1.0e+05 *-0.293080000000000 0.000010000000000-0.489155000000000 0.000010000000000-0.685815000000000 0.000010000000000-0.867300000000000 0.000010000000000-1.078925000000000 0.000010000000000-1.359435000000000 0.000010000000000-1.683695000000000 0.000010000000000-2.048485000000000 0.000010000000000
13、-2.443480000000000 0.000010000000000-2.858050000000000 0.000010000000000-3.311735000000000 0.000010000000000方法二 程序（1）一次累加生成序列的 matlab 命令 x0=19519,19578,19637,19695,16602,25723,30379,34473,38485,40514,42400,48337; x1(1)=x0(1); x1(1)x1(1) =19519 for t=2:12;x1(t)=x1(t-1)+x0(t);endx1 回车x1 =Columns 1 thr
14、ough 819519 39097 58734 78429 95031 120754 151133 185606Columns 9 through 12224091 264605 307005 355342（2）由一次累加生成序列紧邻均值生成 Z 的 matlab 命令：)1(x0=19519,19578,19637,19695,16602,25723,30379,34473,38485,40514,42400,48337; x1(1)=x0(1); for t=2:12;x1(t)=x1(t-1)+x0(t);z1 (t)=(1/2)*(x1(t)+x1(t-1);endz1 =1.0e+05 *Columns 1 through 70 0.2931 0.4892 0.6858 0.8673 1.0789 1.3594Columns 8 through 121.6837 2.0485 2.4435 2.8581 3.3117(3) 由于 GM(1,1)的灰微分方程为 X (t)-+aZ (k)=b)0()1(设 为待估参数， ,利用最小二乘法得到ba, matlab 程序为YB1)（B=-z1(2:12),ones(11,1); Y=(x0(2:12); alpha=inv(B*B)*B*Y; alphaalpha =1.0e+04 *-0.00001.4000
