4行关键代码实现灰色模型GM(1, 1)
文章目录
1、灰色模型GM(1, 1)
先抄书(嫌烦可略过,建议略过):
灰色GM(1, 1)模型:
x
(
0
)
(
k
)
+
a
z
(
1
)
(
k
)
=
b
x^{(0)}(k)+az^{(1)}(k)=b
x(0)(k)+az(1)(k)=b
其白化方程通常写为:
d
x
(
1
)
(
t
)
d
t
+
a
x
(
1
)
(
t
)
=
b
\frac{dx^{(1)}(t)}{dt}+ax^{(1)}(t)=b
dtdx(1)(t)+ax(1)(t)=b
这里
z
(
1
)
(
k
)
z^{(1)}(k)
z(1)(k) 在邓和刘的书中通常称为紧邻均值序列,有时也称为背景值,它的公式为:
z
(
1
)
(
k
)
=
1
2
(
x
(
1
)
(
k
)
+
x
(
1
)
(
k
−
1
)
)
z^{(1)}(k) = \frac{1}{2}\left( x^{(1)}(k) + x^{(1)}(k-1)\right)
z(1)(k)=21(x(1)(k)+x(1)(k−1))
这里
x
(
1
)
(
k
)
x^{(1)}(k)
x(1)(k) 即是原始数据的累加值,表示为
x
(
1
)
(
k
)
=
∑
j
=
1
k
x
(
0
)
(
j
)
x^{(1)}(k) = \sum^k_{j=1}x^{(0)}(j)
x(1)(k)=j=1∑kx(0)(j)
参数估计式:
(
a
^
b
^
)
=
(
B
T
B
)
−
1
B
T
Y
\begin{pmatrix} \hat{a}\\\hat{b} \end{pmatrix} =\left(B^TB\right)^{-1}B^TY
(a^b^)=(BTB)−1BTY
其中:
B
=
(
−
z
(
1
)
(
2
)
1
−
z
(
1
)
(
3
)
1
⋮
⋮
−
z
(
1
)
(
n
)
1
)
,
Y
=
(
x
(
1
)
(
2
)
x
(
1
)
(
3
)
⋮
x
(
1
)
(
n
)
)
B = \begin{pmatrix} -z^{(1)}(2) \quad 1 \\ -z^{(1)}(3) \quad 1 \\ \vdots \quad \vdots\\-z^{(1)}(n) \quad 1 \end{pmatrix} , Y = \begin{pmatrix} x^{(1)}(2) \\ x^{(1)}(3) \\ \vdots \\x^{(1)}(n) \end{pmatrix}
B=⎝⎜⎜⎜⎛−z(1)(2)1−z(1)(3)1⋮⋮−z(1)(n)1⎠⎟⎟⎟⎞,Y=⎝⎜⎜⎜⎛x(1)(2)x(1)(3)⋮x(1)(n)⎠⎟⎟⎟⎞
离散响应式(通过求解白化方程得到,取初值为
x
^
(
1
)
(
1
)
=
x
(
1
)
(
1
)
=
x
(
0
)
(
1
)
\hat{x}^{(1)}(1)=x^{(1)}(1)=x^{(0)}(1)
x^(1)(1)=x(1)(1)=x(0)(1)):
x
^
(
1
)
(
k
)
=
(
x
(
0
)
(
1
)
−
a
^
b
^
)
e
−
a
^
(
k
−
1
)
+
a
^
b
^
\hat{x}^{(1)}(k)= \left(x^{(0)}(1) -\frac{\hat{a}}{\hat{b}} \right)e^{-\hat{a}(k-1)}+\frac{\hat{a}}{\hat{b}}
x^(1)(k)=(x(0)(1)−b^a^)e−a^(k−1)+b^a^
最后计算还原式:
x
^
(
0
)
(
k
)
=
x
^
(
1
)
(
k
)
−
x
^
(
1
)
(
k
−
1
)
\hat{x}^{(0)}(k)=\hat{x}^{(1)}(k)-\hat{x}^{(1)}(k-1)
x^(0)(k)=x^(1)(k)−x^(1)(k−1)
2、重新梳理计算步骤(重点)
许多人通常在看书的时候就会觉得,一会又是矩阵,一会又是方程,感觉有点麻烦。在实操过程中稍遇到些问题,一想直背后太多知识点就放弃,十分可惜。
下面,我们单纯以实现这个模型的计算为目的,重新梳理其整个建模步骤。
数据准备:假设我们有原始数据
[
x
(
0
)
(
1
)
,
x
(
0
)
(
2
)
,
…
,
x
(
0
)
(
n
)
,
…
,
x
(
0
)
(
n
+
p
)
]
\left[ x^{(0)}(1),x^{(0)}(2),\dots,x^{(0)}(n),\dots,x^{(0)}(n+p)\right]
[x(0)(1),x(0)(2),…,x(0)(n),…,x(0)(n+p)]
第一步:计算其累加值:
x
(
1
)
(
k
)
=
∑
j
=
1
k
x
(
0
)
(
j
)
x^{(1)}(k) = \sum_{j=1}^{k}x^{(0)}(j)
x(1)(k)=j=1∑kx(0)(j)
第二步: 再算背景值:
z
(
1
)
(
k
)
=
1
2
(
x
(
1
)
(
k
)
+
x
(
1
)
(
k
−
1
)
)
z^{(1)}(k) = \frac{1}{2}\left( x^{(1)}(k) + x^{(1)}(k-1)\right)
z(1)(k)=21(x(1)(k)+x(1)(k−1))
第三步:构造矩阵
B
B
B (在具体实操中,我们并不需要专门构造矩阵
Y
Y
Y):
B
=
(
−
z
(
1
)
(
2
)
1
−
z
(
1
)
(
3
)
1
⋮
⋮
−
z
(
1
)
(
n
)
1
)
B = \begin{pmatrix} -z^{(1)}(2) & 1 \\ -z^{(1)}(3) &1 \\ \vdots & \vdots \\-z^{(1)}(n) &1 \end{pmatrix}
B=⎝⎜⎜⎜⎛−z(1)(2)−z(1)(3)⋮−z(1)(n)11⋮1⎠⎟⎟⎟⎞
第四步:算 矩阵的转置、 乘积、求逆再算矩阵乘向量(后面马上就来具体实操,并不难):
(
a
^
b
^
)
=
(
B
T
B
)
−
1
B
T
Y
\begin{pmatrix} \hat{a}\\\hat{b} \end{pmatrix} =\left(B^TB\right)^{-1}B^TY
(a^b^)=(BTB)−1BTY
第五步:算时间响应式和还原值(取
k
=
1
,
2
,
…
,
n
,
…
,
n
+
p
k=1,2,\dots,n,\dots,n+p
k=1,2,…,n,…,n+p):
x
^
(
1
)
(
k
)
=
(
x
(
0
)
(
1
)
−
a
^
b
^
)
e
−
a
^
(
k
−
1
)
+
a
^
b
^
\hat{x}^{(1)}(k)= \left(x^{(0)}(1) -\frac{\hat{a}}{\hat{b}} \right)e^{-\hat{a}(k-1)}+\frac{\hat{a}}{\hat{b}}
x^(1)(k)=(x(0)(1)−b^a^)e−a^(k−1)+b^a^
x ^ ( 0 ) ( k ) = x ^ ( 1 ) ( k ) − x ^ ( 1 ) ( k − 1 ) \hat{x}^{(0)}(k)=\hat{x}^{(1)}(k)-\hat{x}^{(1)}(k-1) x^(0)(k)=x^(1)(k)−x^(1)(k−1)
3、MATLAB代码手把手实现(以下高能)
在上一节里我们将模型的主要步骤分成了5个步骤,实际上真正的计算过程就是按这个顺序逐步执行的。接下来,我们就来一步步写出各步的代码。
在以下内容中,你将看到熟悉MATLAB编程特性的情况下,实现这一模型是何其简单。
数据准备:
x0 = exp(0.1*[1:10])' + rand(10,1);
注: 这里是生成了一组10*1的矩阵,以指数函数为真实值,加入0到1之间的噪声。这里我们采用列向量方便后面计算。
另外,我们设 n = 8 , p = 2 n=8,p=2 n=8,p=2
n = 8;
p = 2;
第一步:计算其累加值:
x1 = cumsum(x0);
注: cumsum()
是MATLAB内置函数,直接求出某一组列向量的累加值,返回一个和原数据同样的矩阵。那么也就是说,这段代码直接实现了累加的所有计算,不需要任何循环。
第二步: 计算背景值:
z1 = 0.5 * ( x1(1:n-1) + x1(2:n) );
注: 如果上述操作是小高能,那么这一步就是高能操作。即是利用MATLAB中对矩阵元素的取出方式,简单实现了背景值的计算。同样,不需要任何循环。
第三步: 构造矩阵:
B = [ -z1 ones(n-1,1)];
Y = x0(2:n);
注:如果上面是高能,那么这里就算是小超高能。利用MATLAB分块矩阵的记法,轻松实现两个矩阵的构造,一句循环都不需要,同时代码简洁易懂。
第四步:计算矩阵乘积、求逆、、、、(不想打了,直接看代码):
val_a_b = pinv(B)*Y;
a = val_a_b(1);
b = val_a_b(2);
注: 不明真相的群众需要补充一些知识,即广义逆矩阵(也称为伪逆矩阵)。见百科:
https://zh.wikipedia.org/wiki/%E5%B9%BF%E4%B9%89%E9%80%86%E9%98%B5zh.wikipedia.org
这里: 函数pinv()
即是矩阵广义逆的函数,当矩阵行大于列时,这个函数计算出的值即是:
(
B
T
B
)
−
1
B
T
\left(B^TB\right)^{-1}B^T
(BTB)−1BT
所以我们所看到的这一大坨矩阵计算,一个函数就可以搞定。
第五步: 计算响应式和还原式
这里你以为我要用循环了?(图样图桑破!)看代码:
k = [1:10]';
hat_x1 = ( x0(1) - b/a ) * exp(-a*k) + b/a;
hat_x0 = [ x0(1);hat_x1(2:end)-hat_x1(1:end-1) ];
注:这里属于核能。
这段代码有几个技巧需要说明:
(1) 在计算响应式的时候,我们预先知道需要算出 k = 1 , 2 , . . . , 10 k = 1, 2,...,10 k=1,2,...,10 的值,通常的思路是写一个FOR循环让 k k k取遍1到10. 不过MATLAB提供了十
分灵活的矩阵计算方式。用第二行代码算出的值直接就是响应式所有值构成的列向量,因此,不需要循环。
(2) 在计算还原式的时候,采用了2个技巧。第1个仍然是分块矩阵的运用,这一点不用多说。第2个技巧是MATLAB对下标的操作。这里 hat_x1(2:end)
指的是这一向量中第2个到最后一个元素,同理hat_x1(1:end-1)
即是第1到倒数第二个元素。这两组相减,刚好每个元素的值就是
x
^
(
0
)
(
k
)
=
x
^
(
1
)
(
k
)
−
x
^
(
1
)
(
k
−
1
)
\hat{x}^{(0)}(k)=\hat{x}^{(1)}(k)-\hat{x}^{(1)}(k-1)
x^(0)(k)=x^(1)(k)−x^(1)(k−1) ,所以这里,我们仍然不需要任何循环。
4、福利完整代码:
x0 = exp(0.1*[1:10])' + rand(10,1);
n = 8;
p = 2;
x1 = cumsum(x0);
z1 = 0.5 * ( x1(1:n-1) + x1(2:n) );
B = [ -z1 ones(n-1,1)];
Y = x0(2:n);
val_a_b = pinv(B)*Y;
a = val_a_b(1);
b = val_a_b(2);
k = [1:10]';
hat_x1 = ( x0(1) - b/a ) * exp(-a*k) + b/a;
hat_x0 = [ x0(1);hat_x1(2:end)-hat_x1(1:end-1) ];
至于模型的评估,计算误差这些问题,自己简单写一写即可。
提示:仍然全部使用MATLAB自带函数,如mean()
,abs()
以及向量的表达式。
5、进一步讨论(核能继续)
上述代码其实仍有简化的空间,有几个地方其实是完全没必要的。
(1)矩阵B和Y并没有必要显式构造出来,同时背景值也只用来做了一下矩阵构造,所以,这三行可以不要。
(2)上述三行代码只是用来算参数,因此把这三行合并到这一行来即可。
即:
z1 = 0.5 * ( x1(1:n-1) + x1(2:n) );
B = [ -z1 ones(n-1,1)];
Y = x0(2:n);
val_a_b = pinv(B)*Y;
缩减到一行:
val_a_b = pinv([ -0.5 * ( x1(1:n-1) + x1(2:n) ) ones(n-1,1)] )*x0(2:n);
(3)单独用变量a,b,k
来存储一些值,只是为了让代码更好懂,其实可以全删了。直接用val_a_b(1)
, val_a_b(2)
分别代替a,b
; 后面的k直接写成[1:n+p]'
即可,这里就不拆了。
直接上最终简洁版:
x0 = exp(0.1*[1:10])' + rand(10,1);
n = 8;
p = 2;
x1 = cumsum(x0);
val_a_b = pinv([ -0.5 * ( x1(1:n-1) + x1(2:n) ) ones(n-1,1)] )*x0(2:n);
hat_x1 = ( x0(1) - val_a_b(2)/val_a_b(1) ) * exp(-val_a_b(1)*[1:n+p]') + val_a_b(2)/val_a_b(1);
hat_x0 = [ x0(1);hat_x1(2:end)-hat_x1(1:end-1) ];
仔细看看,前3行是必要的数据准备和参数设置。真正的计算过程4行代码即可。
6、小结及后记
读完本文,如果认真读的话,可能你花得不止10分钟。
读完本文,相信你已经能学会灰色模型的实现,换上你自己的数据试一试,或许会有更多发现。
通过本文,相信你已经能体会到一些MATLAB编程的乐趣,尤其是亲自动手实现过这个模型或者类似模型的朋友。
更多关于灰色预测模型的最新成果,可以参考以下链接:
https://www.researchgate.net/profile/Sifeng_Liuwww.researchgate.netBo Zeng | Chongqing Technology and Business University, Chongqingwww.researchgate.net
https://www.researchgate.net/profile/Xin_Ma58
更多代码也可在MATHWORKS社区中找到:
ng | Chongqing Technology and Business University, Chongqingwww.researchgate.net[外链图片转存中…(img-8ZbaeKRB-1583045245538)]](https://link.zhihu.com/?target=https%3A//www.researchgate.net/profile/Bo_Zeng15)
更多代码也可在MATHWORKS社区中找到:
The kernel-based grey system model - File Exchange - MATLAB Centralww2.mathworks.cn