灰色预测模型的matlab运行代码(讲解)

1、灰色预测模型 GM(1,1)的matlab运行代码例 由19902001年中国蔬菜产量，建立模型预测 2002年中国蔬菜产量，并对预测结果作检验。分析建模：给定原始时间19902001年资料序列X(0)(k),对X(0)(k)生成1-AGO(累加)序列X(1)(k)及Yn。见下表K 123 45 67891011 12X 19519 ,19578 ,19637,19695,16602, 25723, 30379, 34473, 38485, 40514, 42400, 48337X(1) 19519, 39097, 58734,264605,307005,355342Yn19578 19637

2、40514 42400 48337其中X (k):x(0)(i)；i 1Yn = X(0)(2),X(0)(3), ,X(0)(12)TX (0) (k)的 GM(1,1),得到z(1 2)10.5 X (1)X(。1-29308.0z(3)10.5 X (2)X "(3)1-48915.5z(1 )(4)10.5 X (3)X(14)1-68581.5z(5)10.5 X(4)X(15)1-86730.0z(1 )(6)10.5 X (5)X1-107892.5z(7)10.5 X (6)X( "(7)1-135943.5z(8)10.5 X (7)X "(8)

3、1-168369.5z(9)10.5 X (8)X(9)1-204848.5z(10)10.5 X (9)X (10)1244348.0z(11)10.5 X (10)X( "(11)1-331173.5z(12)10.5 X (11)X( "(12)1-5236.21111111111将B和Yn代入辨识算式,有:对上述ka(BtB) 1BT?Ynb得灰色GM(1,1)模型为0.106210513999.91(1)灰微分方程 X(0) (k)-0.1062105 Z (k)=13999.9白化方程 罟 0.1062105X白化方程的时间响应式乂亿 1) X(0)(1) ba

4、13999.9(4)还原为原始数据预测方程:e at)?(0)(t1)151332.5e 0.1062105t131813.5X(1)(t 1)刃(t),即B =刃亿 1)15248.968e0.1062105t(5)残差检验：残差 error1=e1 =)?(i)X(0)(i)，这里残差有12个。相对残差error2=e2=|)?(0) (i)X(0)(i)|X(0)(i)e1XW,这里相对残差有12个。(6)后验差检验:C=S其中S1?(0)(0)12(0)(i) F2,S1为绝对误差序列的标准差。刃 X(0)(i) e1, -(0)丄 12 (0)12 i 1(i)12X( 0)(i)X

5、(0)21 12S2为原始数据系列标准差，S2L , X(0)丄 X(0)(i)n 112 i 1C<0.35好；C<0.5合格；00.6不合格。利用matlab做求解a,b,B,并作残差分析>> x0=19519,19578,19637,19695,16602,25723,30379,34473,38485,40514,42400,48337;>> format long;(表示设计精度)>> n=le ngth(x0);(输入数据长度)>> x1=; (表示x1是一矩阵)>> x1(1)=x0(1);>>

6、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>> alp ha=(B'*B)(-1)*B'*Y'>> a=alpha(1,1);>> b=alpha(2,1);>> d=b/a; (计算时间响应函数参数)>> c=x1(1)-d;>> x2(1)=x0(1);&g

7、t;> x(1)=x0(1);>> for i=1:n-1;x2(i+1)=c*exp(-a*i)+d;(这里x2(i+1 )相当上面所讲的 乂(t 1)x(i+1)=x2(i+1)-x2(i);这里 x(i+1) 相当原来输入数据的预测数据X?(0)(t1)end>> for i=2: 12; x2(i)=c*exp(-a*(i-1)+d; x(i)=x2(i)-x2(i-1);end对上面刚引出的 x2(i) 进行说明及计算)>> for i=1:n;error(i)=x(i)-x0(i); (残差)error1(i)=abs(error(i);

8、(计算残差， abs 表示绝对值) error2(i)=error1(i)/x0(i); (计算相对误差) end>> C=std(error1)/std(x0); (计算后验差检验数 ,std 表示标准差)>> k=1; ( k 表示预测长度，这里每次预测下一年)-0.1772>> b1.399996741173038e+041.0e+05 *-0.29300.0000-0.48900.0000-0.68580.0000-0.86730.0000-1.00000.0000-1.3594350000000000.0000-1.6836950000000000

9、.0000-2.00000.0000-2.4434800000000000.0000-2.85800.0000-3.3117350000000000.0000>> C （求后检验数）0.163969348419772 >> x （原始数据X（0）（i）的对应的预测数据 刃（i），这里也是12个）1.0e+04 *Columns 1 through 31.9519 1.695769385830782 1.885790370699694Columns 4 through 62.4606 2.3321 2.593422554345592Columns 7 through 92

10、.8840 3.2656 3.566589717441855Columns 10 through 123.966247180534730 4.4128 4.9964 >> eroor1 （求残差） eroor1 =1.0e+03 *Columns 1 through 40 2.6285 0.7795 1.2765Columns 5 through 86.7907 0.211225543455919 1.538691164348100 2.4441Columns 9 through 122.8191 0.8596 1.7284 0.712323610019637>> er

11、ror2 （求相对误差）error2 =Columns 1 through 40 0.6267 0.5525 0.5471Columns 5 through 80.4031 0.8968 0.14290.6503Columns 9 through 120.1369 0.8635 0.15160.606-0.17721.399996741173038e+04>> B1.0e+05 *-0.29300.0000-0.48900.0000-0.68580.0000-0.86730.0000-1.00000.0000-1.3594350000000000.0000-1.6836950000

12、000000.0000-2.00000.0000-2.4434800000000000.0000-2.85800.0000-3.3117350000000000.0000方法二 程序matlab 命令(1)一次累加生成序列的>> 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 =Co

13、lumns 1 through 81951939097587347842995031120754151133185606Columns 9 through 122240912646053070053553422)由一次累加生成序列紧邻均值生成Z (1) 的 matlab 命令： 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); end z1 =1.0e+05 *Columns 1 through 70 0.29310.48920.68580.8673 1.0789 1.3594Columns 8 through 121.6837 2.04852.44352.85813.3117(1) (k)=b(3) 由于 G