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)+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=?\begin{array}{c}?z^{(1)}(2)1\\ ?z^{(1)}(3)1\\ ??\\ ?z^{(1)}(n)1\end{array}?,Y=?\begin{array}{c}{x}^{(1)}(2)\\ x^{(1)}(3)\\ ?\\ x^{(1)}(n)\end{array}?$$
離散響應式（通過求解白化方程得到，取初值為${\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=?\begin{array}{cc}?z^{(1)}(2)& 1\\ ?z^{(1)}(3)& 1\\ ?& ?\\ ?z^{(1)}(n)& 1\end{array}?$$
第四步：算 矩陣的轉置、 乘積、求逆再算矩陣乘向量（后面馬上就來具體實操，并不難）：
$$\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

總結

以上是生活随笔為你收集整理的4行关键代码实现灰色模型GM(1, 1)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。