灰色预测- 灰色模型GM(1,N)
上篇博客说到了 G M ( 1 , 1 ) GM(1,1) GM(1,1)，即只有单个变量，且模型是1阶的，此次我们来说说 G M ( 1 , N ) GM(1,N) GM(1,N),即 1 1 1 阶的，包含有 N N N 个变量的灰色模型。如果 G M ( 1 , 1 ) GM(1,1) GM(1,1)弄明白了，那么 G M ( 1 , N ) GM(1,N) GM(1,N)就简单了。
GM(1,N)模型
类似地，定义
N
N
N个变量的原始序列为：
x
i
(
0
)
=
(
x
i
(
0
)
(
1
)
,
x
i
(
0
)
(
2
)
,
.
.
.
,
,
x
i
(
0
)
(
n
)
)
,
i
=
1
,
2
,
.
.
.
,
n
(1)
x^{(0)}_i=(x^{(0)}_i(1), x^{(0)}_i(2),...,, x^{(0)}_i(n) ),i=1,2,...,n \tag{1}
xi(0)=(xi(0)(1),xi(0)(2),...,,xi(0)(n)),i=1,2,...,n(1)
对应的一次累加序列
x
i
(
1
)
x^{(1)}_i
xi(1)为：
x
i
(
1
)
=
(
x
i
(
1
)
(
1
)
,
x
i
(
1
)
(
2
)
,
.
.
.
,
,
x
i
(
1
)
(
n
)
)
,
i
=
1
,
2
,
.
.
.
,
n
(2)
x^{(1)}_i=(x^{(1)}_i(1), x^{(1)}_i(2),...,, x^{(1)}_i(n) ),i=1,2,...,n \tag{2}
xi(1)=(xi(1)(1),xi(1)(2),...,,xi(1)(n)),i=1,2,...,n(2)
其中，
x
i
(
1
)
=
∑
j
=
α
k
x
i
(
0
)
(
j
)
,
k
=
α
,
α
+
1
,
.
.
.
n
(3)
x^{(1)}_i=\sum_{j=\alpha}^k x^{(0)}_i(j),k=\alpha,\alpha+1,...n \tag{3}
xi(1)=j=α∑kxi(0)(j),k=α,α+1,...n(3)
其中
α
≤
n
\alpha \le n
α≤n，且为正整数,
(
3
)
(3)
(3)式中，取
α
=
1
\alpha = 1
α=1,称为一般累加过程，记作
1
−
A
G
O
1-AGO
1−AGO。定义
x
i
(
1
)
(
k
)
x^{(1)}_i(k)
xi(1)(k)的灰导数（实际上就是累减），为：
d
i
(
k
)
=
x
i
(
0
)
(
k
)
=
x
i
(
1
)
(
k
)
−
x
i
(
1
)
(
k
−
1
)
(4)
d_i (k)=x^{(0)}_i(k)=x^{(1)}_i(k)- x^{(1)}_i(k-1) \tag{4}
di(k)=xi(0)(k)=xi(1)(k)−xi(1)(k−1)(4)
取
x
i
(
1
)
x^{(1)}_i
xi(1)的等权重紧邻值，
z
i
(
1
)
(
k
)
=
0.5
x
i
(
1
)
(
k
)
+
0.5
x
i
(
1
)
(
k
−
1
)
，
k
=
2
,
3
,
.
.
.
,
n
(5)
z^{(1)}_i(k)=0.5 x^{(1)}_i(k)+0.5 x^{(1)}_i(k-1) ，k=2,3,...,n \tag{5}
zi(1)(k)=0.5xi(1)(k)+0.5xi(1)(k−1)，k=2,3,...,n(5)
定义灰微分方程
G
M
(
1
,
N
)
GM(1,N)
GM(1,N)：
d
1
(
k
)
+
a
z
1
(
1
)
(
k
)
=
∑
i
=
2
N
b
i
x
i
(
1
)
(
k
)
,
k
=
2
,
3
,
.
.
.
,
n
(6)
d_1(k)+az_1^{(1)}(k)=\sum _{i=2}^N b_i x_i^{(1)} (k), k=2,3,...,n \tag{6}
d1(k)+az1(1)(k)=i=2∑Nbixi(1)(k),k=2,3,...,n(6)
即：
x
1
(
0
)
(
k
)
+
a
z
1
(
1
)
(
k
)
=
∑
i
=
2
N
b
i
x
i
(
1
)
(
k
)
,
k
=
2
,
3
,
.
.
.
,
n
(7)
x^{(0)}_1(k)+az_1^{(1)}(k)=\sum _{i=2}^N b_i x_i^{(1)} (k), k=2,3,...,n \tag{7}
x1(0)(k)+az1(1)(k)=i=2∑Nbixi(1)(k),k=2,3,...,n(7)
注意： 虽然我们在上面定义了
d
i
(
k
)
d_i (k)
di(k),和
z
i
(
1
)
z_i^{(1)}
zi(1)
(
i
=
1
,
2
,
.
.
.
,
n
)
( i=1,2,...,n )
(i=1,2,...,n),但是在灰微分方程中仅仅用到了
i
=
1
i=1
i=1。
同样地，将
k
=
2
,
3
,
.
.
.
,
n
k=2,3,...,n
k=2,3,...,n的数据带入上面的
(
8
)
(8)
(8)式,得到：
{
x
1
(
0
)
(
2
)
+
a
z
1
(
1
)
(
2
)
=
∑
i
=
2
N
b
i
x
i
(
1
)
(
2
)
x
1
(
0
)
(
3
)
+
a
z
1
(
1
)
(
3
)
=
∑
i
=
2
N
b
i
x
i
(
1
)
(
3
)
.
.
.
.
x
1
(
0
)
(
n
)
+
a
z
1
(
1
)
(
n
)
=
∑
i
=
2
N
b
i
x
i
(
1
)
(
n
)
(8)
\left\{ \begin{matrix} x^{(0)}_1(2)+az^{(1)}_1(2)= \sum _{i=2}^Nb_ix_i^{(1)} (2) \\ x^{(0)}_1(3)+az^{(1)}_1(3)= \sum _{i=2}^Nb_ix_i^{(1)} (3) \\ .... \\ x^{(0)}_1(n)+az^{(1)}_1(n)= \sum _{i=2}^Nb_ix_i^{(1)} (n) \\ \end{matrix} \right. \tag{8}
⎩⎪⎪⎪⎨⎪⎪⎪⎧x1(0)(2)+az1(1)(2)=∑i=2Nbixi(1)(2)x1(0)(3)+az1(1)(3)=∑i=2Nbixi(1)(3)....x1(0)(n)+az1(1)(n)=∑i=2Nbixi(1)(n)(8)
将上面写成矩阵的形式，令:
Y
=
(
x
1
(
0
)
(
2
)
,
x
1
(
0
)
(
3
)
,
.
.
.
,
x
1
(
0
)
(
n
)
)
T
Y = (x^{(0)}_1(2),x^{(0)}_1(3),...,x^{(0)}_1(n))^T
Y=(x1(0)(2),x1(0)(3),...,x1(0)(n))T
B
=
(
−
z
1
(
1
)
(
2
)
x
2
(
1
)
(
2
)
…
x
N
(
1
)
(
2
)
−
z
1
(
1
)
(
3
)
x
2
(
1
)
(
3
)
…
x
N
(
1
)
(
3
)
⋮
⋮
⋮
⋮
−
z
1
(
1
)
(
n
)
x
2
(
1
)
(
n
)
…
x
N
(
1
)
(
n
)
)
B =\begin{pmatrix} -z^{(1)}_1(2) & x^{(1)}_2(2) & \dots & x^{(1)}_N(2) \\ -z^{(1)}_1(3) & x^{(1)}_2(3) & \dots& x^{(1)}_N(3) \\ \vdots &\vdots & \vdots& \vdots \\ -z^{(1)}_1(n) & x^{(1)}_2(n) & \dots & x^{(1)}_N(n) \\ \end{pmatrix}
B=⎝⎜⎜⎜⎜⎛−z1(1)(2)−z1(1)(3)⋮−z1(1)(n)x2(1)(2)x2(1)(3)⋮x2(1)(n)……⋮…xN(1)(2)xN(1)(3)⋮xN(1)(n)⎠⎟⎟⎟⎟⎞
u
=
(
a
b
2
b
3
…
b
N
)
T
u = \begin{pmatrix} a &b_2 &b_3 &\dots&b_N \end{pmatrix}^T
u=(ab2b3…bN)T
称
Y
Y
Y为数据向量，
B
B
B 为数据矩阵，
u
u
u 为参数向量，则
G
M
(
1
,
N
)
GM(1,N)
GM(1,N)模型可以表示为矩阵方程:
Y = B u (9) Y = Bu \tag{9} Y=Bu(9)
(
8
)
(8)
(8)式可以看出，参数个数为
N
N
N个，方程个数
(
n
−
1
)
(n-1)
(n−1)一般大于
N
N
N，因此利用最小二乘求解待解未知数
u
u
u。最小二乘（使得
J
(
u
^
)
=
(
Y
−
B
u
^
)
T
(
Y
−
B
u
^
)
J(\hat u)=(Y-B\hat u)^T(Y-B\hat u)
J(u^)=(Y−Bu^)T(Y−Bu^)最小）解为：
u
^
=
(
a
^
b
^
2
b
^
3
…
b
^
N
)
T
=
(
B
T
B
)
−
1
B
T
Y
(10)
\hat u = \begin{pmatrix} \hat a &\hat b_2 & \hat b_3 &\dots&\hat b_N \end{pmatrix}^T=(B^TB)^{-1}B^TY \tag{10}
u^=(a^b^2b^3…b^N)T=(BTB)−1BTY(10)
需要注意的是， ( 10 ) (10) (10)式有解，必须满足 ( B T B ) (B^TB) (BTB)是非奇异矩阵，否则， ( B T B ) − 1 (B^TB)^{-1} (BTB)−1就不存在。但注意到 u ^ \hat u u^ 的元素实际上是各子因素对主因素影响大小的反映，因此，我们可以引入加权矩阵 W = d i a g ( w 1 , w 2 , . . . , w N ) W = diag(w_1,w_2 ,...,w_N ) W=diag(w1,w2,...,wN)，使对各因素的未来发展趋势进行调整控制。对于未来发展减弱趋势的因素赋予较大的权值，而对于未来增强趋势的因素赋予较小的权值，使之更好地反映未来的实际情况。此时，向量 u ^ \hat u u^ 可以表示为：
u ^ = ( a ^ b ^ 2 b ^ 3 … b ^ N ) T = W − 1 B T ( B W − 1 B T ) − 1 Y (11) \hat u = \begin{pmatrix} \hat a &\hat b_2 &\hat b_3 &\dots&\hat b_N \end{pmatrix}^T=W^{-1}B^T(BW^{-1}B^T)^{-1}Y \tag{11} u^=(a^b^2b^3…b^N)T=W−1BT(BW−1BT)−1Y(11)
GM(1,N)的白化型
对于
G
M
(
1
,
N
)
GM(1,N)
GM(1,N)的灰微分方程
(
7
)
(7)
(7)，如果将
x
i
(
0
)
(
k
)
x^{(0)}_i (k)
xi(0)(k) 的时刻
k
=
2
,
3
,
.
.
.
n
k = 2,3,... n
k=2,3,...n视为连续连续的变量 t ，则数列
x
i
(
1
)
x^{(1)}_i
xi(1)就可以视为时间
t
t
t 的函数，记为
x
i
(
1
)
=
x
i
(
1
)
(
t
)
,
i
=
2
,
3
,
.
.
.
N
x^{(1)}_i = x^{(1)}_i (t),i=2,3,...N
xi(1)=xi(1)(t),i=2,3,...N ，并让灰导数
x
1
(
0
)
(
k
)
x^{(0)}_1 (k)
x1(0)(k) 对应于导数
d
1
x
(
1
)
d
t
\frac{d^{x(1)}_1}{dt}
dtd1x(1)，背景值
z
1
(
1
)
(
k
)
z^{(1)}_1 (k)
z1(1)(k) 对应于
x
1
(
1
)
(
t
)
x^{(1)} _1(t)
x1(1)(t) 。于是得到
G
M
(
1
,
1
)
GM(1,1)
GM(1,1)的灰微分方程对应的白微分方程为(
G
M
(
1
,
N
)
GM(1,N)
GM(1,N)的白化型):
d
1
x
(
1
)
d
t
+
a
x
1
(
1
)
(
t
)
=
∑
i
=
2
N
b
i
x
i
(
1
)
(
t
)
(12)
\frac{d^{x(1)}_1}{dt} +ax^{(1)}_1(t)=\sum _{i=2}^Nb_ix_i^{(1)} \tag{12} (t)
dtd1x(1)+ax1(1)(t)=i=2∑Nbixi(1)(t)(12)
需要注意的是，
G
M
(
1
,
N
)
GM(1,N)
GM(1,N)的白化型并不是
G
M
(
1
,
N
)
GM(1,N)
GM(1,N)的灰微分方程离散得到的，仅仅是一种类推。
G
M
(
1
,
N
)
GM(1,N)
GM(1,N)的白化型是一个真正的微分方程，如果白化型模型精度高，则表明所用数列建立的模型
G
M
(
1
,
N
)
GM(1,N)
GM(1,N)与真正的微分方程模型吻合较好，反之亦然。
G
M
(
1
,
N
)
GM(1,N)
GM(1,N)的白化型为一阶 N 个变量的微分方程。
GM(1,N)预测
通过将(
G
M
(
1
,
N
)
GM(1,N)
GM(1,N)的白化型方程离散，得到
x
1
(
1
)
x^{(1)}_1
x1(1)的预测值
x
^
1
(
1
)
\hat x^{(1)}_1
x^1(1),
x
1
(
1
)
x^{(1)}_1
x1(1)是由
x
1
(
0
)
x^{(0)}_1
x1(0)累加得到的，因此可以通过累减
x
^
1
(
1
)
\hat x^{(1)}_1
x^1(1)得到
x
1
(
0
)
x^{(0)}_1
x1(0)的预测值
x
^
1
(
0
)
\hat x^{(0)}_1
x^1(0),即：
x
^
1
(
0
)
(
k
+
1
)
=
x
^
1
(
1
)
(
k
+
1
)
−
x
^
1
(
1
)
(
k
)
(10)
\hat x^{(0)} _1(k+1)= \hat x^{(1)}_1 (k+1)-\hat x^{(1)}_1 (k) \tag{10}
x^1(0)(k+1)=x^1(1)(k+1)−x^1(1)(k)(10)