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

灰色预测模型GM(1,1)的matlab 运行代码

例 由1990—2001年中国蔬菜产量，建立模型预测2002年中国蔬菜产量，并对预测结果作检验。

分析建模：给定原始时间1990—2001年资料序列X )0((k),对X )0((k)生成1-AGO(累加)序列X )1((k)及Y n 。见下表

X )1((k)=)(1

i )

0(i X

k

∑=; Y n =T X X X )]12(,),3(),2([)0()0()0(

对上述X )0((k)的GM(1,1),得到

[

][][][][][][][][][][]

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

⎤⎢

⎢⎢

⎢⎢⎢

⎢⎢⎢

⎢⎢⎢

⎢⎢⎢

⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢

⎢⎢⎢

⎢⎢⎢⎢⎢

⎢⎢⎢

⎢⎣⎡+-+-+-+-+-+-+-+-+-+-+-==⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-----------=15236.2-1331173.5-1244348.01204848.5-1168369.5-1135943.5-1107892.5-186730.0-168581.5-148915.5-129308.0-112115.0111105.011095.01985.01875.01765.01655.01545.01435.01325.01215.01)12(1)11(1)10(1)9(1)8(1)

7(1)6(1)5(1)4(1)3(1)2()1()1()

1()1()1()1()

1()1()1()1()1()1()1()1()1()1()1()1()1()1()

1()1()

1()1()

1()1()1()

1()1()

1()1()

1()1(）（）（

）（）（）（）（）（）（）

（）（）

（）（）

（）（）（）（）

（）（）（）（）（）（X X X X X X X X X X X X X X X X X X X X X X z z z z z z z z z z z B 将 B 和Y n 代入辨识算式，有：

10.1062105()13999.9T T

n a B B B Y b α--⎡⎤⎡⎤==•=⎢⎥⎢⎥⎣⎦⎣⎦

得灰色GM(1,1)模型为

（1）灰微分方程X )0((k)-0.1062105 Z )1((k)=13999.9

(2)白化方程

9.139991062105.0)

1()1(=-X dt

dX (3)白化方程的时间响应式

5.1318135.151332)1()1(ˆ1062105.0)0()1(-=+⎥⎦⎤⎢⎣

⎡-=+--t at e a b e a b X t X (4)还原为原始数据预测方程：)(ˆ)1(ˆ)1(ˆ)1()1()0(t X t X t X

-+=+,即 t e t X

1062105.0)0(15248.968)1(ˆ=+ (5)残差检验： 残差error1=e1=

)()ˆ)0()0(i X i X -（,这里残差有12个。

相对残差error2=e2=)0()

0()0()0(1

)(|)()(ˆ|X

e i X i X i X =-,这里相对残差有12个。

(6) 后验差检验：

C=

2

1

S S , 其中1

])([12

1

2

)0(01-∆-∆

=∑=n i S i ）

（, S1为绝对误差序列的标准差。

1)()(ˆ)()0()0()

0(e i X i X i =-=∆, )(121121

)0()0(i i ∑=∆=∆ S2为原始数据系列标准差，1

])([12

1

2

)0(02--=

∑=n X i X

S i ）

（ , )(121121

)

0()

0(i X X

i ∑== C<0.35 好；C<0.5 合格；C>0.6 不合格。

利用matlab 做求解a,b,B,并作残差分析

>> x0=[19519,19578,19637,19695,16602,25723,30379,34473,38485,40514,42400,48337]; >> format long; (表示设计精度)

>> n=length(x0); （输入数据长度） >> 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)=c*exp(-a*i)+d; （这里x2(i+1)相当上面所讲的)1(ˆ)

1(+t X

）

x(i+1)=x2(i+1)-x2(i); （这里x(i+1) 相当原来输入数据的预测数据)1(ˆ)

0(+t X

）

end

>> 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(error1)/std(x0); （计算后验差检验数,std 表示标准差） >> k=1; （k 表示预测长度，这里每次预测下一年） a =

-0.106210475032772 >> b b =

1.399996741173038e+04 B B =