【数学建模常用算法】之灰色预测模型GM

  • Post author:
  • Post category:其他


作者:張張張張

github地址:https://github.com/zhanghekai

【转载请注明出处,谢谢!】



一、灰色预测模型GM(1,1)


灰色预测是用灰色模型GM(1,1)来进行定量分析的,通常分为以下几类:


(1) 灰色时间序列预测。用等时距观测到的反映预测对象特征的一系列数量(如产量、销量、人口数量、存款数量、利率等)构造灰色预测模型,预测未来某一时刻的特征量,或者达到某特征量的时间。

(2) 畸变预测(灾变预测)。通过模型预测异常值出现的时刻,预测异常值什么时候出现在特定时区内。

(3) 波形预测,或称为拓扑预测,它是通过灰色模型预测事物未来变动的轨迹。

(4) 系统预测,对系统行为特征指标建立一族相互关联的灰色预测理论模型,在预测系统整体变化的同时,预测系统各个环节的变化。


灰色模型预测流程

  1. 数据检验与数据预处理
  2. 模型构建
  3. 模型检验



1、数据检验与数据预处理



1.1 构建模型前检验




\qquad











给定观测数据列



x

(

0

)

=

{

x

1

(

0

)

,

x

2

(

0

)

,


 

,

x

n

(

0

)

}

x^{(0)}=\{x_1^{(0)},x_2^{(0)},\cdots,x_n^{(0)}\}







x











(


0


)












=








{




x










1









(


0


)



















,





x










2









(


0


)



















,











,





x










n









(


0


)



















}





一般采用



x

(

0

)

x^{(0)}







x











(


0


)













的级比



λ

k

(

0

)

\lambda^{(0)}_k







λ










k









(


0


)






















的大小与所属区间来判断。




\qquad











级比定义:



λ

k

(

0

)

=

x

k

1

(

0

)

x

k

(

0

)

,

k

=

2

,

3

,


 

,

n

\lambda^{(0)}_k=\frac{x^{(0)}_{k-1}}{x^{(0)}_k},\quad k=2,3,\cdots,n







λ










k









(


0


)





















=





















x










k









(


0


)


































x











k





1










(


0


)






































,






k




=








2


,




3


,











,




n









\qquad











若满足



λ

k

(

0

)

(

e

2

n

+

1

,

e

2

n

+

1

)

\lambda^{(0)}_k\in (e^{-\frac{2}{n+1}},e^{\frac{2}{n+1}})







λ










k









(


0


)






























(



e


























n


+


1
















2





























,





e























n


+


1
















2





























)





,则认为



x

(

0

)

x^{(0)}







x











(


0


)













是可作为



G

M

(

1

,

1

)

GM(1,1)






G


M


(


1


,




1


)





建模的。如果原始序列不满足,需要对原始序列做必要的变换处理,使其落入可容覆盖内。即取适当的常数



C

C






C





,做平移变换:





y

(

0

)

=

x

i

(

0

)

+

C

,

i

=

1

,

2

,


 

,

n

y^{(0)}=x^{(0)}_i+C,\quad i=1,2,\cdots,n







y











(


0


)












=









x










i









(


0


)





















+








C


,






i




=








1


,




2


,











,




n










\qquad











使序列



y

(

0

)

=

(

y

1

(

0

)

,

y

2

(

0

)

,


 

,

y

n

(

0

)

)

y^{(0)}=(y^{(0)}_1,y^{(0)}_2,\cdots,y^{(0)}_n)







y











(


0


)












=








(



y










1









(


0


)



















,





y










2









(


0


)



















,











,





y










n









(


0


)



















)





的级比



λ

y

(

k

)

=

y

k

1

(

0

)

y

(

k

)

\lambda_y^{(k)}=\frac{y^{(0)}_{k-1}}{y^{(k)}}







λ










y









(


k


)





















=





















y











(


k


)

























y











k





1










(


0


)









































落入上述区间范围。



1.2 数据生成



(1) 累加生成(AGO)




\qquad











设原始序列为



x

(

0

)

=

(

x

1

(

0

)

,

x

2

(

0

)

,


 

,

x

n

(

0

)

)

x^{(0)}=(x^{(0)}_1,x^{(0)}_2,\cdots,x^{(0)}_n)







x











(


0


)












=








(



x










1









(


0


)



















,





x










2









(


0


)



















,











,





x










n









(


0


)



















)





,令:




x

k

(

1

)

=

i

=

1

k

x

i

(

0

)

,

k

=

1

,

2

,


 

,

n

x^{(1)}_k=\sum_{i=1}^{k}x^{(0)}_i,\quad k=1,2,\cdots ,n







x










k









(


1


)





















=

















i


=


1



















k





















x










i









(


0


)



















,






k




=








1


,




2


,











,




n






有:



x

(

1

)

=

(

x

1

(

1

)

,

x

2

(

1

)

,


 

,

x

n

(

1

)

)

x^{(1)}=(x^{(1)}_1,x^{(1)}_2,\cdots,x^{(1)}_n)







x











(


1


)












=








(



x










1









(


1


)



















,





x










2









(


1


)



















,











,





x










n









(


1


)



















)





,称所得到的新数列



x

(

1

)

x^{(1)}







x











(


1


)













为数列



x

(

0

)

x^{(0)}







x











(


0


)













的1次累加生成的数列。



(2)累减生成(IAGO)




\qquad











通过累加数列得到的新数列,可以通过累减生成还原出原始数列,设原始数列



x

(

1

)

=

(

x

1

(

1

)

,

x

2

(

1

)

,


 

,

x

n

(

1

)

)

x^{(1)}=(x^{(1)}_1,x^{(1)}_2,\cdots,x^{(1)}_n)







x











(


1


)












=








(



x










1









(


1


)



















,





x










2









(


1


)



















,











,





x










n









(


1


)



















)





,令:





x

k

(

0

)

=

(

x

k

(

1

)

x

k

1

(

1

)

)

,

k

=

2

,

3

,

n

x^{(0)}_k=(x^{(1)}_k-x^{(1)}_{k-1}),\quad k=2,3,\cdots n







x










k









(


0


)





















=








(



x










k









(


1


)































x











k





1










(


1


)



















)


,






k




=








2


,




3


,









n






称所得到的数列



x

(

0

)

x^{(0)}







x











(


0


)

















x

(

1

)

x^{(1)}







x











(


1


)













的1次累减生成数列。



(3)加权邻值生成




\qquad















z

(

1

)

z^{(1)}







z











(


1


)

















x

(

1

)

x^{(1)}







x











(


1


)













的紧邻均值生成序列:





z

(

1

)

=

(

z

2

(

1

)

,

z

3

(

1

)

,


 

,

z

n

(

1

)

)

z^{(1)}=(z^{(1)}_2,z^{(1)}_3,\cdots,z^{(1)}_n)







z











(


1


)












=








(



z










2









(


1


)



















,





z










3









(


1


)



















,











,





z










n









(


1


)



















)






其中,





z

k

(

1

)

=

α

x

k

(

1

)

+

(

1

α

)

x

k

1

(

1

)

,

k

=

2

,

3

,


 

,

n

z^{(1)}_k=\alpha x^{(1)}_k+(1-\alpha)x^{(1)}_{k-1},\quad k=2,3,\cdots,n







z










k









(


1


)





















=








α



x










k









(


1


)





















+








(


1













α


)



x











k





1










(


1


)



















,






k




=








2


,




3


,











,




n









α

\alpha






α





称为生成系数,当



α

=

0.5

\alpha = 0.5






α




=








0


.


5





时,则称该数列为均值生成数,也称为等权邻值生成数。



2 、灰色模型GM(1,1)构建



2.1 定义



x

(

1

)

x^{(1)}







x











(


1


)













的灰导数为





d

k

=

x

k

(

0

)

=

x

k

(

1

)

x

k

1

(

1

)

d_k=x^{(0)}_k=x^{(1)}_k-x^{(1)}_{k-1}







d










k




















=









x










k









(


0


)





















=









x










k









(


1


)































x











k





1










(


1


)
























于是GM(1,1)的灰微分方程模型为:






d

k

+

a

z

k

(

1

)

=

b

x

k

(

0

)

+

a

z

k

(

1

)

=

b

(

1

)

d_k+a z^{(1)}_k = b \quad或 \quad x^{(0)}_k+a z^{(1)}_k = b\qquad\qquad\qquad(1)







d










k




















+








a



z










k









(


1


)





















=








b










x










k









(


0


)





















+








a



z










k









(


1


)





















=








b








(


1


)







其中,



x

k

(

0

)

x^{(0)}_k







x










k









(


0


)






















称为灰导数,



α

\alpha






α





称为发展系数,



z

k

(

1

)

z^{(1)}_k







z










k









(


1


)






















称为白化背景值,



b

b






b





称为灰作用量。将时刻



k

=

2

,

3

,


 

,

n

k=2,3,\cdots,n






k




=








2


,




3


,











,




n





代入上式,有:





{

x

2

(

0

)

+

a

z

2

(

1

)

=

b

x

3

(

0

)

+

a

z

3

(

1

)

=

b

x

n

(

0

)

+

a

z

n

(

1

)

=

b

\begin{cases}x^{(0)}_2+a z^{(1)}_2=b\\ x^{(0)}_3+a z^{(1)}_3=b\\ \qquad\quad\cdots\\ x^{(0)}_n+a z^{(1)}_n=b \end{cases}






















































































































x










2









(


0


)





















+




a



z










2









(


1


)





















=




b









x










3









(


0


)





















+




a



z










3









(


1


)





















=




b






















x










n









(


0


)





















+




a



z










n









(


1


)





















=




b



























引入矩阵向量记号:





μ

=

[

a

b

]

,

Y

=

[

x

2

(

0

)

x

3

(

0

)

x

n

(

0

)

]

,

B

=

[

z

2

(

1

)

1

z

3

(

1

)

1

z

n

(

1

)

1

]

\mu=\begin{bmatrix}a\\b\end{bmatrix} ,Y=\begin{bmatrix}x^{(0)}_2\\x^{(0)}_3\\ \vdots\\x^{(0)}_n\end{bmatrix},B=\begin{bmatrix}-z^{(1)}_2\quad 1\\-z^{(1)}_3\quad 1\\\vdots\\-z^{(1)}_n\quad 1 \end{bmatrix}






μ




=










[













a








b




















]






,




Y




=





























































































x










2









(


0


)


























x










3









(


0


)







































x










n









(


0


)
















































































































,




B




=
































































































z










2









(


1


)





















1












z










3









(


1


)





















1

























z










n









(


1


)





















1



































































































于是GM(1,1)模型可以表示为



Y

=

B

μ

Y=B\mu






Y




=








B


μ







那么现在的问题就是求



a

a






a









b

b






b





的值,可以用最小二乘法求它们的估计值为:





μ

=

[

a

b

]

=

(

B

T

B

)

1

B

T

Y

\mu=\begin{bmatrix}a\\b\end{bmatrix}=(B^TB)^{-1}B^TY






μ




=










[













a








b




















]






=








(



B










T









B



)














1











B










T









Y







2. 2 GM(1,1)的白化模型(也叫影子模型)




\qquad











对于



G

M

(

1

,

1

)

GM(1,1)






G


M


(


1


,




1


)





的灰微分方程。如果将时刻



k

=

2

,

3

,


 

,

n

k=2,3,\cdots,n






k




=








2


,




3


,











,




n





视为连续变量



t

t






t





,则之前的



x

(

1

)

x^(1)







x










(









1


)





视为时间



t

t






t





函数,于是灰导数



x

(

0

)

x^{(0)}







x











(


0


)













变为连续函数的导数



d

x

t

(

1

)

d

t

\frac{dx^{(1)}_t}{dt}


















d


t
















d



x










t









(


1


)









































,白化背景值



z

k

(

1

)

z^{(1)}_k







z










k









(


1


)






















对应于导数



x

t

(

1

)

x^{(1)}_t







x










t









(


1


)






















。于是GM(1,1)的灰微分方程对应于的白微分方程为:





d

x

t

(

1

)

d

t

+

a

x

t

(

1

)

=

b

(

2

)

\frac{dx^{(1)}_t}{dt}+ax^{(1)}_t=b\qquad\qquad\qquad(2)

















d


t














d



x










t









(


1


)







































+








a



x










t









(


1


)





















=








b








(


2


)







解为:





x

t

(

1

)

=

(

x

1

(

0

)

b

a

)

e

a

(

t

1

)

+

b

a

x^{(1)}_t=(x^{(0)}_1-\frac{b}{a})e^{-a(t-1)}+\frac{b}{a}







x










t









(


1


)





















=








(



x










1









(


0


)









































a














b




















)



e














a


(


t





1


)












+



















a














b
























注1:以往的文献中并没有对这种从离散形式到连续形式的直接跳跃做出解释。

于是得到预测值:





x

^

k

+

1

(

1

)

=

(

x

1

(

0

)

b

a

)

e

a

k

+

b

a

,

k

=

1

,

2

,


 

,

n

1

\hat{x}^{(1)}_{k+1}=(x^{(0)}_1-\frac{b}{a})e^{-ak}+\frac{b}{a},\qquad k=1,2,\cdots,n-1















x







^
















k


+


1










(


1


)





















=








(



x










1









(


0


)









































a














b




















)



e














a


k












+



















a














b




















,






k




=








1


,




2


,











,




n













1






从而得到相应的预测值:




x

^

k

+

1

(

0

)

=

x

^

k

+

1

(

1

)

x

^

k

(

1

)

,

k

=

1

,

2

,


 

,

n

1

\hat{x}^{(0)}_{k+1}=\hat{x}^{(1)}_{k+1}-\hat{x}^{(1)}_k,\qquad k=1,2,\cdots,n-1















x







^
















k


+


1










(


0


)





















=

















x







^
















k


+


1










(


1


)







































x







^















k









(


1


)



















,






k




=








1


,




2


,











,




n













1






注2:原始序列数据不一定要全部使用,相应建立的模型也会不同,即a和b不同。



注3:原始序列数据必须要等时间间隔、不间断。



3、检验预测值



3.1 残差检验:计算相对残差





ξ

k

=

x

k

(

0

)

x

^

k

(

0

)

x

k

(

0

)

,

k

=

1

,

2

,


 

,

n

\xi_k=\frac{x^{(0)}_k-\hat{x}^{(0)}_k}{x^{(0)}_k},\qquad k=1,2,\cdots,n







ξ










k




















=




















x










k









(


0


)
































x










k









(


0


)



































x







^















k









(


0


)





































,






k




=








1


,




2


,











,




n






如果对所有的



ξ

k

<

0.1

|\xi_k|<0.1










ξ










k























<








0


.


1





,则认为达到较高的要求;

如果对所有的



ξ

k

&lt;

0.2

|\xi_k|&lt;0.2










ξ










k























<








0


.


2





,则认为达到一般的要求。



3.2 级比偏差检验

根据前面计算出来的级比



λ

k

\lambda_k







λ










k





















和发展系数



a

a






a





,计算相应的级比偏差:





ρ

k

=

1

1

0.5

a

1

+

0.5

a

λ

k

\rho_k=1-\frac{1-0.5a}{1+0.5a}\lambda_k







ρ










k




















=








1
























1




+




0


.


5


a














1









0


.


5


a





















λ










k


























ρ

k

&lt;

0.1

|\rho_k|&lt;0.1










ρ










k























<








0


.


1





,则认为达到较高要求;





ρ

k

&lt;

0.2

|\rho_k|&lt;0.2










ρ










k























<








0


.


2





,则认为达到一般要求;



二、使用灰色模型进行预测案例

北方某城市1986~1992年道路交通噪声平均声级数据见表1,为近年来交通噪声数据【dB(A)】

序号年份 eq L
1986 71.1
1987 72.4
1988 72.4
1989 72.1
1990 71.4
1991 72.0
1992 71.6


第一步:级比检验





\quad











建立交通噪声平均声级数据时间序列:




x

(

0

)

=

(

x

1

(

0

)

,

x

2

(

0

)

,


&ThinSpace;

,

x

7

(

0

)

)

=

(

71.1

,

72.4

,

72.4

,

72.1

,

71.4

,

72.0

,

71.6

)

x^{(0)}=(x^{(0)}_1,x^{(0)}_2,\cdots,x^{(0)}_7)=(71.1,72.4,72.4,72.1,71.4,72.0,71.6)







x











(


0


)












=








(



x










1









(


0


)



















,





x










2









(


0


)



















,











,





x










7









(


0


)



















)




=








(


7


1


.


1


,




7


2


.


4


,




7


2


.


4


,




7


2


.


1


,




7


1


.


4


,




7


2


.


0


,




7


1


.


6


)








\quad











(1)求级比



λ

k

\lambda_k







λ










k

























λ

k

=

x

k

1

(

0

)

x

k

(

0

)

\lambda_k=\frac{x^{(0)}_{k-1}}{x^{(0)}_k}







λ










k




















=





















x










k









(


0


)


































x











k





1










(


0


)














































λ

=

(

λ

2

,

λ

3

,


&ThinSpace;

,

λ

7

)

=

(

0.982

,

1

,

1.0042

,

1.0098

,

0.9917

,

1.0056

)

\lambda=(\lambda_2,\lambda_3,\cdots,\lambda_7)=(0.982,1,1.0042,1.0098,0.9917,1.0056)






λ




=








(



λ










2


















,





λ










3


















,











,





λ










7


















)




=








(


0


.


9


8


2


,




1


,




1


.


0


0


4


2


,




1


.


0


0


9


8


,




0


.


9


9


1


7


,




1


.


0


0


5


6


)








\quad











(2)级比判断:



λ

k

(

0

)

(

e

2

n

+

1

,

e

2

n

+

1

)

\lambda^{(0)}_k\in (e^{-\frac{2}{n+1}},e^{\frac{2}{n+1}})







λ










k









(


0


)






























(



e


























n


+


1
















2





























,





e























n


+


1
















2





























)






由于所有的



λ

k

[

0.982

,

1.0098

]

,

k

=

2

,

3

,


&ThinSpace;

,

7

\lambda_k\in[0.982,1.0098],k=2,3,\cdots,7







λ










k





























[


0


.


9


8


2


,




1


.


0


0


9


8


]


,




k




=








2


,




3


,











,




7





则可以用



x

(

0

)

x^{(0)}







x











(


0


)













进行GM(1,1)建模。


第二步:GM(1,1)建模





\quad











(1)对原始数据



x

(

0

)

x^{(0)}







x











(


0


)













作一次累加





x

(

1

)

=

(

71.1

,

143.5

,

215.9

,

288

,

359.4

,

431.4

,

503

)

x^{(1)}=(71.1,143.5,215.9,288,359.4,431.4,503)







x











(


1


)












=








(


7


1


.


1


,




1


4


3


.


5


,




2


1


5


.


9


,




2


8


8


,




3


5


9


.


4


,




4


3


1


.


4


,




5


0


3


)








\quad











(2)构造数据矩阵B及数据向量Y





B

=

[

1

2

(

x

1

(

1

)

+

x

2

(

1

)

)

1

1

2

(

x

2

(

1

)

+

x

3

(

1

)

)

1

1

2

(

x

6

(

1

)

+

x

7

(

1

)

)

1

]

,

Y

=

[

x

2

(

0

)

x

3

(

0

)

x

7

(

0

)

]

B=\begin{bmatrix}-\frac{1}{2}(x^{(1)}_1+x^{(1)}_2)\qquad 1\\[1ex]-\frac{1}{2}(x^{(1)}_2+x^{(1)}_3)\qquad 1\\[1ex]\vdots\\[1ex]-\frac{1}{2}(x^{(1)}_6+x^{(1)}_7)\qquad 1\end{bmatrix},\quad Y=\begin{bmatrix}x^{(0)}_2\\[1ex]x^{(0)}_3\\[1ex]\vdots\\[1ex] x^{(0)}_7\end{bmatrix}






B




=





























































































































2
















1





















(



x










1









(


1


)





















+





x










2









(


1


)



















)




1























2
















1





















(



x










2









(


1


)





















+





x










3









(


1


)



















)




1




































2
















1





















(



x










6









(


1


)





















+





x










7









(


1


)



















)




1

















































































































,






Y




=















































































































x










2









(


0


)


























x










3









(


0


)







































x










7









(


0


)






































































































































\quad











(3)计算



μ

\mu






μ










μ

=

[

a

b

]

=

(

B

T

B

)

1

B

T

Y

=

[

0.0023

72.6573

]

\mu=\begin{bmatrix}a\\b\end{bmatrix}=(B^TB)^{-1}B^TY=\begin{bmatrix}0.0023\\72.6573\end{bmatrix}






μ




=










[













a








b




















]






=








(



B










T









B



)














1











B










T









Y




=










[













0


.


0


0


2


3








7


2


.


6


5


7


3




















]










\quad











(4)建立模型





d

x

(

1

)

d

t

+

0.0023

x

(

1

)

=

72.6573

\frac{dx^{(1)}}{dt}+0.0023x^{(1)}=72.6573

















d


t














d



x











(


1


)






























+








0


.


0


0


2


3



x











(


1


)












=








7


2


.


6


5


7


3










\quad\qquad













求解得





x

^

k

+

1

(

1

)

=

(

x

1

(

0

)

b

a

)

e

a

k

+

b

a

=

30929

e

0.0023

k

+

31000

\hat{x}^{(1)}_{k+1}=(x^{(0)}_1-\frac{b}{a})e^{-ak}+\frac{b}{a}=-30929e^{-0.0023k}+31000















x







^
















k


+


1










(


1


)





















=








(



x










1









(


0


)









































a














b




















)



e














a


k












+



















a














b






















=











3


0


9


2


9



e














0


.


0


0


2


3


k












+








3


1


0


0


0








\quad











(5)求生成数列值



x

^

k

+

1

(

1

)

\hat{x}^{(1)}_{k+1}















x







^
















k


+


1










(


1


)






















及模型还原值



x

^

k

+

1

(

0

)

\hat{x}^{(0)}_{k+1}















x







^
















k


+


1










(


0


)






















,令:



k

=

1

2

6

k=1,2,\cdots ,6






k




=








1





2















6





,由上面的时间响应函数可算得



x

^

(

1

)

\hat{x}^{(1)}















x







^
















(


1


)













,其中取




x

^

1

(

1

)

=

x

^

1

(

0

)

=

x

1

(

0

)

=

71.1

\hat{x}^{(1)}_1=\hat{x}^{(0)}_1=x^{(0)}_1=71.1















x







^















1









(


1


)





















=

















x







^















1









(


0


)





















=









x










1









(


0


)





















=








7


1


.


1










x

^

k

(

0

)

=

x

^

k

(

1

)

x

^

k

1

(

1

)

\hat{x}^{(0)}_k=\hat{x}^{(1)}_k-\hat{x}^{(1)}_{k-1}















x







^















k









(


0


)





















=

















x







^















k









(


1


)







































x







^
















k





1










(


1


)






















,取



k

=

2

,

3

,


&ThinSpace;

,

7

k=2,3,\cdots,7






k




=








2


,




3


,











,




7





,得





x

^

(

0

)

=

(

x

^

1

(

0

)

x

^

2

(

0

)

x

^

7

(

0

)

)

=

71.1

72.4

72.2

72.1

71.9

71.7

71.6

\hat{x}^{(0)}=(\hat{x}^{(0)}_1,\hat{x}^{(0)}_2,\cdots,\hat{x}^{(0)}_7)=(71.1,72.4,72.2,72.1,71.9,71.7,71.6)















x







^
















(


0


)












=








(











x







^















1









(


0


)































x







^















2









(


0


)









































x







^















7









(


0


)



















)




=











7


1


.


1





7


2


.


4





7


2


.


2





7


2


.


1





7


1


.


9





7


1


.


7





7


1


.


6









第三步:模型检验





\quad











模型的各种检验指标的计算结果如下表:

年份 原始值 模型值 残差 相对误差 级比偏差
1986 71.1 71.1
1987 72.4 72.4 -0.0057 0.01% 0.0023
1988 72.4 72.2 0.1638 0.23% 0.0203
1989 72.1 72.1 0.0329 0.05% -0018
1990 71.4 71.9 -0.4984 0.7% -0.0074
1991 72.0 71.7 0.2699 0.37% 0.0107
1992 71.6 71.6 0.0378 0.05% -0.0032




\quad











经验证,该模型的精度较高,可进行预测和预报。


【参考文献】

  • 徐华锋,方志耕《优化白化方程的GM(1, 1)模型》
  • 灰色预测模型GM(1,1) 与例题分:https://blog.csdn.net/qq547276542/article/details/77865341
  • 【数学建模】灰色预测及Python实现:https://www.jianshu.com/p/a35ba96d852b



版权声明:本文为qq_41995574原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。