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

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

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

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

灰色GM(1, 1)模型：
$$x^{(0)}(k)+a{z}^{(1)}(k)=b$$

其白化方程通常写为:
$$\frac{d{x}^{(1)}(t)}{dt}+a{x}^{(1)}(t)=b$$
这里$z^{(1)}(k)$ 在邓和刘的书中通常称为紧邻均值序列，有时也称为背景值，它的公式为:
$$z^{(1)}(k)=\frac{1}{2}\left(x^{(1)}(k)+x^{(1)}(k-1)\right)$$
这里 $x^{(1)}(k)$ 即是原始数据的累加值，表示为
$$x^{(1)}(k)=\sum _{j=1}^k{x}^{(0)}(j)$$
参数估计式：
$$\left(\begin{array}{c}\hat{a}\\ \hat{b}\end{array}\right)={\left(B^{T}B\right)}^{-1}{B}^{T}Y$$
其中：
$$B=\left(\begin{array}{c}-z^{(1)}(2)1\\ -z^{(1)}(3)1\\ ⋮⋮\\ -z^{(1)}(n)1\end{array}\right),Y=\left(\begin{array}{c}{x}^{(1)}(2)\\ x^{(1)}(3)\\ ⋮\\ x^{(1)}(n)\end{array}\right)$$
离散响应式（通过求解白化方程得到，取初值为${\hat{x}}^{(1)}(1)=x^{(1)}(1)=x^{(0)}(1)$):
$${\hat{x}}^{(1)}(k)=\left(x^{(0)}(1)-\frac{\hat{a}}{\hat{b}}\right)e^{-\hat{a}(k-1)}+\frac{\hat{a}}{\hat{b}}$$
最后计算还原式：
$${\hat{x}}^{(0)}(k)={\hat{x}}^{(1)}(k)-{\hat{x}}^{(1)}(k-1)$$

2、重新梳理计算步骤（重点）

许多人通常在看书的时候就会觉得，一会又是矩阵，一会又是方程，感觉有点麻烦。在实操过程中稍遇到些问题，一想直背后太多知识点就放弃，十分可惜。

下面，我们单纯以实现这个模型的计算为目的，重新梳理其整个建模步骤。

数据准备：假设我们有原始数据
$$\left[x^{(0)}(1),x^{(0)}(2),\dots ,x^{(0)}(n),\dots ,x^{(0)}(n+p)\right]$$
第一步：计算其累加值：
$$x^{(1)}(k)=\sum _{j=1}^k{x}^{(0)}(j)$$
第二步: 再算背景值：
$$z^{(1)}(k)=\frac{1}{2}\left(x^{(1)}(k)+x^{(1)}(k-1)\right)$$
第三步：构造矩阵$B$ (在具体实操中，我们并不需要专门构造矩阵$Y$)：
$$B=\left(\begin{array}{cc}-z^{(1)}(2)& 1\\ -z^{(1)}(3)& 1\\ ⋮& ⋮\\ -z^{(1)}(n)& 1\end{array}\right)$$
第四步：算 矩阵的转置、 乘积、求逆再算矩阵乘向量（后面马上就来具体实操，并不难）：
$$\left(\begin{array}{c}\hat{a}\\ \hat{b}\end{array}\right)={\left(B^{T}B\right)}^{-1}{B}^{T}Y$$
第五步：算时间响应式和还原值（取 $k=1,2,\dots ,n,\dots ,n+p$）：
$${\hat{x}}^{(1)}(k)=\left(x^{(0)}(1)-\frac{\hat{a}}{\hat{b}}\right)e^{-\hat{a}(k-1)}+\frac{\hat{a}}{\hat{b}}$$

$${\hat{x}}^{(0)}(k)={\hat{x}}^{(1)}(k)-{\hat{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;

第一步：计算其累加值：

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()即是矩阵广义逆的函数，当矩阵行大于列时，这个函数计算出的值即是：
$${\left(B^{T}B\right)}^{-1}{B}^T$$
所以我们所看到的这一大坨矩阵计算，一个函数就可以搞定。

第五步: 计算响应式和还原式

这里你以为我要用循环了？（图样图桑破！）看代码：

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$ 的值，通常的思路是写一个FOR循环让$k$取遍1到10. 不过MATLAB提供了十

分灵活的矩阵计算方式。用第二行代码算出的值直接就是响应式所有值构成的列向量，因此，不需要循环。

​ (2) 在计算还原式的时候，采用了2个技巧。第1个仍然是分块矩阵的运用，这一点不用多说。第2个技巧是MATLAB对下标的操作。这里 hat_x1(2:end) 指的是这一向量中第2个到最后一个元素，同理hat_x1(1:end-1)即是第1到倒数第二个元素。这两组相减，刚好每个元素的值就是${\hat{x}}^{(0)}(k)={\hat{x}}^{(1)}(k)-{\hat{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