數學建模方法-灰色預測法

一、引言

  大家好，今天是第三天，前兩篇博文用到的數學建模方法應用的數學知識還比較少，在這一篇章，需要的數學知識變多了，大致需要一些高等數學中導數部分的知識，大家如果忘了的話記得自己去補哈。好言歸正傳，今天要講的到底是什么方法呢，就是那個很出名的“灰色預測模型”啦。好了開講咯！

二、灰色預測模型

2.1概念

2.1.1 灰色系統、白色系統、黑色系統

  在講灰色預測模型的概念時我想先跟大家介紹三種系統。

• 白色系統：指一個系統的內部特征的完全已知的，即系統信息是完全公開的。
• 黑色系統：與白色系統相反，指一個系統的內部信息對外界來說是一無所知的，只能通過它與外界的聯系來加以觀測研究。
• 灰色系統：介於白色系統和黑色系統之間，即系統內一部分信息是已知的，但另一部分信息是未知的，系統內各因素間有不確定的關系。

  一般我們在生活的過程中，大多遇到的系統是灰色系統。因此，我們就需要利用灰色系統內那部分公開的信息，來建模推測出隱藏的信息。這就是灰色預測模型。

2.1.2 灰色預測法

• 灰色預測法是一種預測灰色系統的預測方法。

  目前常用的一些預測方法（如回歸分析等），需要較大的樣本，對於少樣本的情況就會造成比較大的誤差，使預測目標失效。而灰色預測模型所需的建模信息少、運算方便、建模精度高，因此在各種預測領域都有着廣泛的應用，是處理小樣本預測問題的有效工具。

• 灰色預測通過鑒別系統因素之間發展趨勢的相異程度，即進行關聯分析，並對原始數據進行灰色生成來尋找系統變動的規律，生成有較強規律性的數據序列，然后建立相應的微分方程模型，從而預測事物未來發展趨勢的狀況。

2.2 關聯分析

  首先什么是關聯分析呢，實際上，關聯分析時動態過程發展態勢的量化比較分析。所謂發展態勢比較，也就是系統各時期有關統計數據的集合關系的比較。

舉例來講，某地區2010-2018年的總收入與養雞和養豬收入資料如下表格。

	2010	2011	2012	2013	2014	2015	2016	2017	2018
總收入（萬元）	16	21	22	21	21	20	20	19	25
養豬（萬元）	12	15	14	16	14	16	13	10	15
養雞（萬元）	4	6	8	5	7	4	7	9	10

  根據上述表格，可以畫出曲線圖，如下：

 

  由上圖易看出，藍色曲線 (總收入)與橙色曲線(養豬收入)發展趨勢比較接近，而與灰色曲線(養雞收入)相差較大，因此可以判斷，該地區對總收入影響較直接的是養豬業， 而不是養雞業。

  很顯然，幾何形狀越接近，關聯程度也就越大。但是，對於一些稍微復雜的問題，如果不容易畫出圖，那就很難通過直觀的手法來分析。這時就需要我們的灰色生成方法了。

2.3 灰色生成

  原始數據列中的數據，按某種要求作數據處理，這樣的一種方法叫做生成。

  灰色系統理論認為，盡管客觀表表象復制，但必定蘊含某種內在規律的，關鍵如何去發現這樣的一個規律。灰色生成就是為了從灰色序列中去找出內在的規律，從而能判斷事物的發展趨勢。

  常用的灰色生成法有：累加生成、累減生成和加權累加生成。接下來我們就來以此介紹它們。

（1）累加生成（AGO）

  假定原始數列為$x^{(0)} = \left (x^{(0)}(1), x^{(0)}(2), \cdot \cdot \cdot ,x^{(0)}(n) \right )$，令

$x^{(1)}(k) = \sum_{i=1}^{k} x^{(0)}(i) \ , \ k = 1, 2,\cdot \cdot \cdot, n$

$x^{(1)} = \left ( x^{1} (1), x^{1} (2), \cdot \cdot \cdot , x^{1} (n) \right )$

  則稱$x^{(1)}$為數列$x^{(0)}$的1次累加生成數列。若繼續推廣，有

$x^{(r)}(k) = \sum_{i =1}^{k} x^{(r-1)}(i) \ , \ k = 1, 2,\cdot \cdot \cdot, n, \ r \geq 1$

  $x^{(r)}(k)$稱為$x^{(0)}(k)$的r次累加生成數列。

（2）累減生成（IAGO）

  假定原始數列為$x^{(1)} = \left (x^{(1)}(1), x^{(1)}(2), \cdot \cdot \cdot ,x^{(1)}(n) \right )$，令

$x^{(0)}(k) = x^{(1)}(k) - x^{(1)}(k-1) \ , \ k = 2, 3,\cdot \cdot \cdot, n$

  則稱$x^{(0)}$為數列$x^{(1)}$的1次累減生成數列。可以看出，通過累加數列得到的新數列，可以通過累減生成還原出原始數列。

（3）加權鄰值生成

  假定原始數列為$x^{(1)} = \left (x^{(1)}(1), x^{(1)}(2), \cdot \cdot \cdot ,x^{(1)}(n) \right )$，稱任意一對相鄰元素$x^{(0)}(k-1), x^{(0)}(k)$互為鄰值。對於常數$\alpha \in [0, 1]$, 令

$z^{(0)}(k) = \alpha x^{(0)}(k) + (1 - \alpha)x^{(0)} (k-1), \ k = 2, 3, \cdot \cdot \cdot, n$

  由此得到的數列稱為鄰值生成數，權$\alpha$也稱為生成系數。特別地，當生成系數$\alpha = 0.5$，則稱該數列為均值生成數，也稱為等權鄰值生成數。

2.4 灰色模型GM(1,1)

  灰色模型是利用離散隨機數經過生成變為隨機性被顯著削弱而且較有規律的生成數，建立起的微分方程形式的模型，這樣便於對其變化過程進行研究和描述。

  G表示grey（灰色），M表示Model（模型）；

  首先，定義$x^{(1)}$的灰導數為

$d(k) = x^{(0)}(k)  = x^{(1)}(k) - x^{(1)}(k-1)$

  令$z^{(1)}$為數列$x^{(1)}$的鄰值生成數列，即

$z^{(1)}(k) = \alpha x^{(1)}(k) + (1 - \alpha)x^{(1)} (k-1),$

  好的。接下來我們要定義最重要的一個方程了，就是GM(1,1)的灰微分方程模型：

$d(k) + a z^{(1)}(k) = b$ 或 $x^{(0)}(k) + a z^{(1)}(k) = b$

 其中，$x^{(0)}(k)$稱為灰倒數，$a$稱為發展系數，$z^{(1)}(k)$稱為白化背景紙，b稱為灰作用量。

  將$k = 2,3,...,n$代入上式，有

 

$x^{(0)}(2) + a z^{(1)}(2) = b\\ $

$x^{(0)}(3) + a z^{(1)}(3) = b\\ $
$\cdot \cdot \cdot\\ $
$x^{(0)}(n) + a z^{(1)}(n) = b$

  引入矩陣向量記號：

$u = \begin{bmatrix} a\\ b\end{bmatrix}$ $ Y = \begin{bmatrix} x^{(0)}(2)\\ x^{(0)}(3)\\  \cdot \cdot \cdot\\  x^{(0)}(2) \end{bmatrix}$ $ B = \begin{bmatrix} -z^{(1)}(2) & 1\\  -z^{(1)}(3) & 1\\  \cdot \cdot \cdot & \cdot \cdot \cdot\\ -z^{(1)}(n) & 1 \end{bmatrix}$

  於是GM(1,1)模型可以表示為

$Y = Bu$

  那么現在的問題就是求a和b的值，我們可以用一元線性回歸，也就是最小二乘法求它們的估計值為

$u = \begin{bmatrix}a\\ b \end{bmatrix} = (B^{T}B)^{-1}B^{T}Y$

2.5 GM(1,1)的白化型

  前面我們定義的模型是在離散的情況下，那么針對連續時間$t$，則之前的$X^{(1)}$視為時間$t$函數，於是灰導數$x^{(0)}(k)$變為連續的導數$\frac{dx^{(1)}(t)}{dt}$，白化背景值$z^{(1)}(k)$對應於導數$x^{(1)}(t)$。於是GM(1,1)的灰微分方程對應的白微分方程為：

$\frac{dx^{(1)}(t)}{dt} + ax^{(1)}(t) = b$

三、GM(1,1)灰色預測的步驟

3.1 數據的檢驗與處理

  為了保證GM(1,1)建模方法的可行性，需要對已知數據做必要的檢驗處理。

  設原始數據列為$x^{(0)} = \left (x^{(0)}(1), x^{(0)}(2), \cdot \cdot \cdot ,x^{(0)}(n) \right )$，計算數列的級比

$\lambda(k) = \frac{x^{(0)}(k-1)}{x^{(0)}(k)},\ k = 2, 3, ..., n$

  如果所有的級比都落在可容覆蓋區間$X = \left(e^{\frac{-2}{n+1}},e^{\frac{2}{n+1}} \right )$內，則數列$x^{(0)}$可以建立GM(1,1)模型且可以進行灰色預測。否則，對數據做適當的變換處理，如平移變換：

$y^{(0)}(k) = x^{(0)}(k) + c, \ k = 1, 2, \cdot \cdot \cdot, n$

  取$c$使得數據列的級比都落在可容覆蓋內。

3.2 建立GM(1,1)模型

  當$x^{(0)} = \left (x^{(0)}(1), x^{(0)}(2), \cdot \cdot \cdot ,x^{(0)}(n) \right )$滿足上面的條件后，我們就可以以它為數據列建立GM(1,1)模型，即根據$Y = Bu$求解出$u = \begin{bmatrix} a\\ b\end{bmatrix}$。相應的白化模型為

$\frac{dx^{(1)}(t)}{dt} + ax^{(1)}(t) = b$

這是一個微分方程，可以求解得到：

$x^{(1)}(t) = \left(x^{(0)}(1) - \frac{b}{a} \right )e^{(-a(t-1))}+ \frac{b}{a}$

於是得到

$\hat{x}^{(1)}(k+1) = \left(x^{(0)}(1) - \frac{b}{a} \right )e^{(-ak)}+ \frac{b}{a},\ k = 1, 2, \cdot \cdot \cdot , n-1$

因此得到預測值為：

$\hat{x}^{(0)}(k) = \hat{x}^{(1)}(k) - \hat{x}^{(1)}(k-1) \ , \ k = 2, 3,\cdot \cdot \cdot, n$

3.3 檢測預測值

  模型選定后，一定要經過檢驗才能判定其是否合理，只有通過檢驗的模型才能用來作預測，灰色模型的精度檢驗一般有三種方法：相對誤差大小檢驗法、關聯度檢驗法和后驗差檢驗法。下面主要介紹相對誤差大小檢驗法和后驗差檢驗法

3.3.1 相對誤差大小檢驗法

$\varepsilon (k) = \frac{x^{(0)}(k) - \hat{x}^{(0)}(k)}{x^{(0)}(k)}, \ k = 1, 2, \cdot \cdot \cdot, n$

  如果對所有的$\left | \varepsilon (k) \right | < 0.1$，則認為到達較高的要求；若對所有的$\left | \varepsilon (k) \right | < 0.2$，則認為打到一般要求。

3.3.2 后驗差檢驗法（方差C檢驗法）

• 原始數列：$x^{(0)} = \left (x^{(0)}(1), x^{(0)}(2), \cdot \cdot \cdot ,x^{(0)}(n) \right )$
• 建模后算出的數列：$\hat{x}^{(0)} = \left (\hat{x}^{(0)}(1), \hat{x}^{(0)}(2), \cdot \cdot \cdot ,\hat{x}^{(0)}(n) \right )$
• 殘差序列為$E = \left (e(1), e(2), \cdot \cdot \cdot ,e(n) \right )$

  其中，$e(k) = x^{(0)}(k) - \hat{x}^{(0)}(k), \ k = 1, 2, \cdot \cdot \cdot , n$

  記原始序列$x^{(0)}$及殘差序列$E$的方差分別為$S_{1}^{2}$和$S_{2}^{2}$，則

$S_{1}^{2} = \frac{1}{n}\sum_{k =1}^{n}\left [ x^{(0)}(k) - \overline{x} \right ]^2$

$S_{2}^{2} = \frac{1}{n}\sum_{k =1}^{n}\left [ e(k) - \overline{e} \right ]^2$

  其中，$\overline{x} = \frac{1}{n}\sum_{k =1}^{n}x^{(0)}(k)$，$\overline{e} = \frac{1}{n}\sum_{k =1}^{n}e(k)$

  則后驗差比為：

$C = \frac {S_{2}}{S_{1}}$

  指標$C$是后驗差檢驗的重要指標。指標$C$越小越好。為什么這么說呢？這是因為，$C$越小就表明$S_{1}$大而$S_{2}$越小。

  那$S_{1}$大表明什么，它表明原始數據方差大，這也就說明原始數據的離散程度大。

  那$S_{2}$小又表明什么呢，它表明殘方差小，也就是殘差的離散程度小。

  因此，$C$小就表明盡管原始數據很離散，但模型所得的計算值與實際值之差並不太離散。這樣得到的數據可靠性就越強。這樣講是不是就懂啦。。

  下面附一張精度檢驗等級參照表

模型精度等級	均方差比值$C$
1級（好）	$C\leq 0.35$
2級（合格）	$0.35 < C \leq 0.5$
3級（勉強）	$0.5 < C \leq 0.65$
4級（不合格）	$0.65 < C$

四、GM（1,1）模型的Matlab代碼

%%構建GM模型計算出a和b

%建立符號變量a(發展系數)和b(灰作用量)
syms a b;
u = [a b]';

%原始數列X0
X0 = input('請輸入數據');
n = length(X0);

%對原始數列X做累加得到數列X1
X1 = cumsum(X0);

%對數列X1做等權鄰值生成Z1
for i = 2: n
	Z1(i) = (X1(i) + X1(i - 1)) / 2;
end
Z1(1) = [];	%去除掉第一個元素

%構造數據矩陣
B = [-Z1; ones(1, n - 1)]';
Y = [X0]'; 
Y(1, :) = [];	%去除掉第一個元素

%利用u = inv(B'*B)*B'*Y來求解
u = inv(B'*B) * B' * Y;
a = u(1, :);
b = u(2, :);


%%預測后續數據

%p是預測后續的數據量，根據情況可更改
p = 3;
%構建預測數列F
F = []; F(1) = X0(1);
for i = 2: (n + p)
　　F(i) = (X0(1) - b/a) * exp(-a*(i-1)) + b/a;
end

%對數列F累減還原，得到預測出的數據
G = []; G(1) = X0(1);
for i = 2: (n + p)
　　G(i) = F(i) - F(i-1);
end

disp('預測到的數據為： ');
G

%%模型檢驗

H = G(1: n);
%計算殘差序列
E = X0 - H;

%法一：相對殘差Q檢驗

%計算相對誤差序列
epsilon = abs(E ./ X0);
%計算相對誤差Q
disp('相對殘差Q檢驗： ');
Q = mean(epsilon)

%法二：方差比C檢驗
disp('方差比C檢驗： ');
C = std(E, 1) / std(X0, 1)

%小誤差概率P檢驗
S1 = std(X0, 1);
tmp = find(abs(E - mean(E)) < 0.6745 * S1);
disp('小誤差概率P檢驗：');
P = length(tmp) / n

%%繪制曲線圖（根據不同的題目以下的代碼會有所不同,故不貼出來，大家自己動手寫吧）