目录
- 灰色预测
- 一阶灰色方程GM(1,1)
- 建模步骤
- 应用及其求解步骤
- 求级比
- 一次累加序列
- 求参数矩阵 u u u
- 时间响应式
- 求预测序列
- 模型检验
- 实际值与预测值比较及可视化
- 二阶灰色方程GM(2,1)
灰色预测
灰色预测模型是通过少量的、不完全的信息，建立数学模型做出预测的预测方法，基于客观事物的过去和现在的发展规律，借助于科学的方法对未来的发展趋势和状况进行描述和分析，并形成科学的假设和判断。
一阶灰色方程GM(1,1)
对于原始数据序列 x ( 0 ) = [ x ( 0 ) ( 1 ) , x ( 0 ) ( 2 ) , . . . , x ( 0 ) ( n ) ] x^{(0)}=[x^{(0)}(1),x^{(0)}(2),...,x^{(0)}(n)] x(0)=[x(0)(1),x(0)(2),...,x(0)(n)]
建模步骤
- 计算级比
已知序列的级比为
σ ( k ) = x ( k ) x ( k − 1 ) , k = ( 2 , 3 , . . . , n ) \sigma(k)=\frac{x(k)}{x(k-1)},k=(2,3,...,n) σ(k)=x(k−1)x(k),k=(2,3,...,n)
若所有级比 σ ( k ) \sigma(k) σ(k)都落在区间 ( e − 2 n + 1 , e 2 n + 1 ) (e^{-\frac{2}{n+1}},e^{\frac{2}{n+1}}) (e−n+12,en+12)内，则序列 x ( 0 ) x^{(0)} x(0)可使用GM(1,1)进行灰色预测，否则应对序列数据 x ( 0 ) x^{(0)} x(0)进行平移变换
x ( 0 ) ( k ) = x ( 0 ) ( k ) + c , k = 1 , 2 , . . . , n x^{(0)}(k)=x^{(0)}(k)+c,k=1,2,...,n x(0)(k)=x(0)(k)+c,k=1,2,...,n
遗留的问题：原始序列进行累加后为什么不会影响到最终结果？
- 累加序列及其均值生成序列
一次累加序列
x ( 1 ) = ( x ( 1 ) ( 1 ) , x ( 1 ) ( 2 ) , . . . , x ( 1 ) ( 7 ) ) = ( x ( 0 ) ( 1 ) , x ( 0 ) ( 1 ) + x ( 0 ) ( 2 ) , . . . , x ( 0 ) ( 1 ) + . . . + x ( 0 ) ( 7 ) ) \begin{aligned} x^{(1)}&=(x^{(1)}(1),x^{(1)}(2),...,x^{(1)}(7))\\ &=(x^{(0)}(1),x^{(0)}(1)+x^{(0)}(2),...,x^{(0)}(1)+...+x^{(0)}(7)) \end{aligned} x(1)=(x(1)(1),x(1)(2),...,x(1)(7))=(x(0)(1),x(0)(1)+x(0)(2),...,x(0)(1)+...+x(0)(7))
一次累加序列的均值生成序列为
z ( 1 ) = ( z ( 1 ) ( 2 ) , z ( 1 ) ( 3 ) , . . . , z ( 1 ) ( 7 ) ) 其中 z ( 1 ) ( k ) = α x ( 1 ) ( k ) + α x ( 1 ) ( k − 1 ) , k = ( 2 , 3 , . . . , 7 ) z^{(1)}=(z^{(1)}(2),z^{(1)}(3),...,z^{(1)}(7))\\ 其中z^{(1)}(k)=\alpha x^{(1)}(k)+\alpha x^{(1)}(k-1),k=(2,3,...,7) z(1)=(z(1)(2),z(1)(3),...,z(1)(7))其中z(1)(k)=αx(1)(k)+αx(1)(k−1),k=(2,3,...,7)
注：均值生成序列中的权重参数取值范围为 0 ≤ α ≤ 1 0\leq \alpha \leq 1 0≤α≤1，一般取0.5
- 建立模型
白化方程
d x ( 1 ) ( t ) d t − a z ( 1 ) ( t ) = b \frac{dx^{(1)}(t)}{dt}-az^{(1)}(t)=b dtdx(1)(t)−az(1)(t)=b
灰微分方程
x ( 0 ) ( k ) − a z k ( 1 ) = b 即 − a z k ( 1 ) + b = x ( 0 ) ( k ) (1) x^{(0)}(k)-az_k^{(1)}=b\\ 即-az_k^{(1)}+b=x^{(0)}(k)\tag{1} x(0)(k)−azk(1)=b即−azk(1)+b=x(0)(k)(1)
其实在这个地方还是有一些遗留问题的，比如白化方程、灰微分方程各是怎么得到的？建立这两个方程背后的数学原理？
将(1)式化为矩阵形式为
[
−
z
(
1
)
(
2
)
1
−
z
(
1
)
(
3
)
1
⋮
1
−
z
(
1
)
(
k
)
1
]
[
a
b
]
=
[
x
(
0
)
(
2
)
x
(
0
)
(
3
)
⋮
x
(
0
)
(
k
)
]
\begin{bmatrix} -z^{(1)}(2)&1\\ -z^{(1)}(3)&1\\ \vdots&1\\ -z^{(1)}(k)&1\\ \end{bmatrix} \begin{bmatrix} a\\b \end{bmatrix}=\begin{bmatrix} x^{(0)}(2)\\ x^{(0)}(3)\\ \vdots\\ x^{(0)}(k) \end{bmatrix}
−z(1)(2)−z(1)(3)⋮−z(1)(k)1111
[ab]=
x(0)(2)x(0)(3)⋮x(0)(k)
即
X
u
=
Y
Xu=Y
Xu=Y,根据最小二乘法得到参数矩阵
u
u
u的估计值为
u
=
(
X
T
X
)
−
1
X
T
Y
u=(X^TX)^{-1}X^TY
u=(XTX)−1XTY
代入数据求得参数矩阵
u
u
u，由此得到白化方程与灰微分方程的参数
a
,
b
a,b
a,b
那么白化方程的时间响应式：
x
(
1
)
(
k
+
1
)
=
(
x
(
0
)
(
1
)
−
b
a
)
e
−
a
k
+
b
a
，其中
k
=
0
,
1
,
.
.
.
x^{(1)}(k+1)=(x^{(0)}(1)-\frac{b}{a})e^{-ak}+\frac{b}{a}，其中k=0,1,...
x(1)(k+1)=(x(0)(1)−ab)e−ak+ab，其中k=0,1,...
将k代入上式，得到一次累减序列
x
m
i
n
u
s
(
1
)
=
[
x
m
i
n
u
s
(
1
)
(
1
)
,
x
m
i
n
u
s
(
1
)
(
2
)
,
.
.
.
,
x
m
i
n
u
s
(
1
)
(
n
)
]
x^{(1)}_{minus}=[x^{(1)}_{minus}(1),x_{minus}^{(1)}(2),...,x_{minus}^{(1)}(n)]
xminus(1)=[xminus(1)(1),xminus(1)(2),...,xminus(1)(n)]
那么预测序列为
x
p
r
e
d
(
0
)
=
[
x
p
r
e
d
(
0
)
(
1
)
,
x
p
r
e
d
(
0
)
(
2
)
,
.
.
.
,
x
p
r
e
d
(
0
)
(
k
)
]
,
k
=
1
,
2
,
.
.
.
其中
x
p
r
e
d
(
0
)
(
k
)
=
x
m
i
n
u
s
(
0
)
(
k
)
−
x
m
i
n
u
s
(
0
)
(
k
−
1
)
x
p
r
e
d
(
0
)
(
1
)
=
x
m
i
n
u
s
(
1
)
(
1
)
x^{(0)}_{pred}=[x^{(0)}_{pred}(1),x^{(0)}_{pred}(2),...,x^{(0)}_{pred}(k)],k=1,2,...\\其中x^{(0)}_{pred}(k)=x^{(0)}_{minus}(k)-x^{(0)}_{minus}(k-1)\\ x^{(0)}_{pred}(1)=x^{(1)}_{minus}(1)
xpred(0)=[xpred(0)(1),xpred(0)(2),...,xpred(0)(k)],k=1,2,...其中xpred(0)(k)=xminus(0)(k)−xminus(0)(k−1)xpred(0)(1)=xminus(1)(1)
5. 误差检验
（1）相对误差检验
λ
(
k
)
=
∣
x
(
0
)
(
k
)
−
x
p
r
e
d
(
0
)
(
k
)
∣
x
(
0
)
(
k
)
,
k
=
1
,
2
,
.
.
.
,
n
\lambda(k)=\frac{|x^{(0)}(k)-x^{(0)}_{pred}(k)|}{x^{(0)}(k)},k=1,2,...,n
λ(k)=x(0)(k)∣x(0)(k)−xpred(0)(k)∣,k=1,2,...,n
当
λ
(
k
)
<
0.2
时
\lambda (k)<0.2时
λ(k)<0.2时，可认为模型达到了一般要求
当
λ
(
k
)
<
0.1
时
\lambda (k)<0.1时
λ(k)<0.1时，可认为模型达到了较高要求
（2）级比偏差值检验
ϵ
(
k
)
=
∣
1
−
(
1
−
0.5
a
1
+
0.5
a
)
σ
(
k
)
∣
,
k
=
2
,
3
,
.
.
.
,
n
\epsilon(k)=|1-(\frac{1-0.5a}{1+0.5a})\sigma (k)|,k=2,3,...,n
ϵ(k)=∣1−(1+0.5a1−0.5a)σ(k)∣,k=2,3,...,n
当
ϵ
(
k
)
<
0.2
时
\epsilon(k)<0.2时
ϵ(k)<0.2时，可认为模型达到了一般要求
当
ϵ
(
k
)
<
0.1
时
\epsilon(k)<0.1时
ϵ(k)<0.1时，可认为模型达到了较高要求
应用及其求解步骤
已知
x
(
0
)
=
(
25723
,
30379
,
34473
,
38485
,
40514
,
42400
,
48337
)
,
试建立
G
M
(
1
,
1
)
模型
已知x^{(0)}=(25723,30379,34473,38485,40514,42400,48337),试建立GM(1,1)模型
已知x(0)=(25723,30379,34473,38485,40514,42400,48337),试建立GM(1,1)模型
建模开始
求级比
import numpy as np
x0=np.array([25723,30379,34473,38485,40514,42400,48337])
sigma=x0[:-1]/x0[1:]
>>> array([0.84673623, 0.88124039, 0.89575159, 0.94991855, 0.95551887,0.87717484])
bound1=[sigma.min(),sigma.max()]
>>> [0.8467362322657098, 0.9555188679245283]
bound2=[np.exp(-2/(x0.size+1)),np.exp(2/(x0.size+1))]
>>> [0.7788007830714049, 1.2840254166877414]
由于bound1 ∈ \in ∈bound2，故级比符合要求
一次累加序列
x ( 1 ) = ( x ( 1 ) ( 1 ) , x ( 1 ) ( 2 ) , . . . , x ( 1 ) ( 7 ) ) = ( x ( 0 ) ( 1 ) , x ( 0 ) ( 1 ) + x ( 0 ) ( 2 ) , . . . , x ( 0 ) ( 1 ) + . . . + x ( 0 ) ( 7 ) ) \begin{aligned} x^{(1)}&=(x^{(1)}(1),x^{(1)}(2),...,x^{(1)}(7))\\ &=(x^{(0)}(1),x^{(0)}(1)+x^{(0)}(2),...,x^{(0)}(1)+...+x^{(0)}(7))\\ \end{aligned} x(1)=(x(1)(1),x(1)(2),...,x(1)(7))=(x(0)(1),x(0)(1)+x(0)(2),...,x(0)(1)+...+x(0)(7))
x1=np.cumsum(x0)
>>> array([ 25723, 56102, 90575, 129060, 169574, 211974, 260311],dtype=int32)
一次累加序列的均值生成序列为
z
(
1
)
=
(
z
(
1
)
(
2
)
,
z
(
1
)
(
3
)
,
.
.
.
,
z
(
1
)
(
7
)
)
其中
z
(
1
)
(
k
)
=
0.5
x
(
1
)
(
k
)
+
0.5
x
(
1
)
(
k
−
1
)
,
k
=
(
2
,
3
,
.
.
.
,
7
)
z^{(1)}=(z^{(1)}(2),z^{(1)}(3),...,z^{(1)}(7))\\ 其中z^{(1)}(k)=0.5x^{(1)}(k)+0.5x^{(1)}(k-1),k=(2,3,...,7)
z(1)=(z(1)(2),z(1)(3),...,z(1)(7))其中z(1)(k)=0.5x(1)(k)+0.5x(1)(k−1),k=(2,3,...,7)
z1=(x1[:-1]+x1[1:])/2.0
>>> array([ 40912.5, 73338.5, 109817.5, 149317. , 190774. , 236142.5])
求参数矩阵 u u u
X = [ − z ( 1 ) ( 2 ) 1 − z ( 1 ) ( 3 ) 1 ⋮ 1 − z ( 1 ) ( k ) 1 ] Y = [ x ( 0 ) ( 2 ) x ( 0 ) ( 3 ) ⋮ x ( 0 ) ( k ) ] X=\begin{bmatrix} -z^{(1)}(2)&1\\ -z^{(1)}(3)&1\\ \vdots&1\\ -z^{(1)}(k)&1\\ \end{bmatrix}Y=\begin{bmatrix} x^{(0)}(2)\\ x^{(0)}(3)\\ \vdots\\ x^{(0)}(k) \end{bmatrix} X= −z(1)(2)−z(1)(3)⋮−z(1)(k)1111 Y= x(0)(2)x(0)(3)⋮x(0)(k)
X=np.vstack([-z1,np.ones(x0.size-1)]).T
Y=x0[1:].T
u=np.matmul(np.matmul((np.linalg.inv(np.matmul(X.T,X))),X.T),Y)
>>> array([-8.42648092e-02, 2.78584508e+04])
时间响应式
于是，白化方程的时间响应式为
x
(
1
)
(
k
+
1
)
=
(
x
(
0
)
(
1
)
−
b
a
)
e
−
a
k
+
b
a
=
(
25723
−
27858.45
−
0.0843
)
e
0.0843
k
−
27858.45
−
0.0843
=
356328.99
e
0.0843
k
−
330605.99
\begin{aligned} x^{(1)}(k+1)&=(x^{(0)}(1)-\frac{b}{a})e^{-ak}+\frac{b}{a}\\ &=(25723-\frac{27858.45}{-0.0843})e^{0.0843k}-\frac{27858.45}{-0.0843}\\ &=356328.99e^{0.0843k}-330605.99 \end{aligned}
x(1)(k+1)=(x(0)(1)−ab)e−ak+ab=(25723−−0.084327858.45)e0.0843k−−0.084327858.45=356328.99e0.0843k−330605.99
求预测序列
#生成匿名函数
x_m=lambda k:(x0[0]-(u[1]/u[0]))*np.e**(-u[0]*k)+u[1]/u[0]
# 一次累减序列
x_m1=x_m(np.arange(x0.size+1))
# 预测序列
x_pred=np.diff(x_m1)
>>> array([31327.3567122 , 34081.56224489, 37077.90911705, 40337.68565577,
43884.05179286, 47742.20361062, 51939.55235392])
模型检验
相对误差
λ
(
k
)
=
∣
x
(
0
)
(
k
)
−
x
p
r
e
d
(
0
)
(
k
)
∣
x
(
0
)
(
k
)
,
k
=
1
,
2
,
.
.
.
,
n
\lambda(k)=\frac{|x^{(0)}(k)-x^{(0)}_{pred}(k)|}{x^{(0)}(k)},k=1,2,...,n
λ(k)=x(0)(k)∣x(0)(k)−xpred(0)(k)∣,k=1,2,...,n
lambda_k=abs(x0[1:]-x_pred[:-1])/x0[1:]
>>> array([0.03121751, 0.01135491, 0.03656206, 0.00435194, 0.03500122,
0.0123052 ])
相对误差值 λ ( k ) < 0.1 \lambda(k)<0.1 λ(k)<0.1，可认为模型达到了较高要求
级比偏差值
ϵ
(
k
)
=
∣
1
−
(
1
−
0.5
a
1
+
0.5
a
)
σ
(
k
)
∣
,
k
=
2
,
3
,
.
.
.
,
n
\epsilon(k)=|1-(\frac{1-0.5a}{1+0.5a})\sigma (k)|,k=2,3,...,n
ϵ(k)=∣1−(1+0.5a1−0.5a)σ(k)∣,k=2,3,...,n
epsilon_k=abs(1-((1-0.5*(-0.0834))/(1+0.5*(-0.0834)))*sigma)
>>> array([0.07957306, 0.04206604, 0.02629194, 0.03258912, 0.03867683,
0.04648542])
级比偏差值 ϵ ( k ) < 0.1 \epsilon(k)<0.1 ϵ(k)<0.1，可认为模型达到了较高要求
实际值与预测值比较及可视化
import matplotlib
import matplotlib.pyplot as plt
x_pred1=x_pred.tolist()
x_pred1.insert(0,x0[0])
matplotlib.rcParams['font.family']='SimHei'
plt.plot(range(x0.size),x0,'r-',label='原始序列')
plt.plot(range(x0.size),x_pred1[:-1],'g--',label='预测序列')
plt.legend()
plt.show()
运行结果
二阶灰色方程GM(2,1)
To be continue…
