4行关键代码实现灰色模型GM(1, 1)

文章目录

4行关键代码实现灰色模型GM(1, 1)
- 1、灰色模型GM(1, 1)
- 2、重新梳理计算步骤（重点）
- 3、MATLAB代码手把手实现（以下高能）
- 4、福利完整代码：
- 5、进一步讨论（核能继续）
- 6、小结及后记

1、灰色模型GM(1, 1)

先抄书（嫌烦可略过，建议略过）：

灰色GM(1, 1)模型：
x(0)(k)+az(1)(k)=bx^{(0)}(k)+az^{(1)}(k)=b x(0)(k)+az(1)(k)=b

其白化方程通常写为:
dx(1)(t)dt+ax(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)=12(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=1kx(0)(j)x^{(1)}(k) = \sum^k_{j=1}x^{(0)}(j) x(1)(k)=j=1∑kx(0)(j)
参数估计式：
(a^b^)=(BTB)−1BTY\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=1kx(0)(j)x^{(1)}(k) = \sum_{j=1}^{k}x^{(0)}(j) x(1)(k)=j=1∑kx(0)(j)
第二步: 再算背景值：
z(1)(k)=12(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))
第三步：构造矩阵BBB (在具体实操中，我们并不需要专门构造矩阵YYY)：
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^)=(BTB)−1BTY\begin{pmatrix} \hat{a}\\\hat{b} \end{pmatrix} =\left(B^TB\right)^{-1}B^TY (a^b^)=(BTB)−1BTY
第五步：算时间响应式和还原值（取 k=1,2,…,n,…,n+pk=1,2,\dots,n,\dots,n+pk=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=2n=8,p=2n=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()即是矩阵广义逆的函数，当矩阵行大于列时，这个函数计算出的值即是：
(BTB)−1BT\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,...,10k = 1, 2,...,10k=1,2,...,10 的值，通常的思路是写一个FOR循环让kkk取遍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