【數學建模】數模day13-灰色系統理論I-灰色關聯與GM(1,1)預測

接下來學習灰色系統理論。

0.

什么是灰色系統？

部分信息已知而部分信息未知的系統，我們稱之為灰色系統。相應的，知道全部信息的叫白色系統，完全未知的叫黑色系統。

為什么采用灰色系統理論？

在給定信息不多，並且無法建立客觀的物理原型，其作用原理亦不明確，內部因素難以辨識或之間關系隱蔽，人們很難准確了解這類系統的行為特征，因此對其定量描述難度較大。這時就采用“灰色系統理論”。

比如說，社會、經濟、農業、生態問題的系統中，噪聲普遍存在，一般受隨機侵蝕的系統理論立足於【概率統計】，比如回歸分析、方差分析、主成分分析等等。但是這些在小樣本（數據不足）、樣本沒有較好的統計分布規律、難以量化等問題下，都不能夠很好的勝任。尤其是，涉及到預測問題時，直接回歸方程代入數得”預測“明顯不符合客觀規律，而使用灰色預測（通常使用GM(1,1)）更可靠。

1.

關聯分析

這個方法解決的是：因素之間關聯性如何，關聯程度如何量化的問題。

討論因素之間關聯性如何，之前我們采用【回歸分析】，即因變量對自變量求回歸方程，這是基於更多樣本的量化討論。為了做【整體系統分析】，得到一個好的【直觀過程】，以及為了【定性描述】，可以考慮采用【關聯分析】。

實際上，實際應用中，我們可以【關聯分析】+【回歸分析】一起做。

關聯分析實際上是動態過程發展態勢的量化比較分析。

簡而言之，關聯分析是從整體態勢上把握兩（多）變量（每個變量的不同樣本構成數列）之間的相關程度，並且從整體上分析減少了異常點的影響。
所謂發展態勢比較，也就是系統各時期有關統計數據的幾何關系的比較。

1) 做關聯分析首先討論數據變換技術。

通過數據變換消除【量綱】，使其具有【可比性】。

2) 做關聯分析：

3) 關聯分析案例：

分析求解：

依照問題的要求，我們自然選取鉛球運動員專項成績作為參考數列，將其余的各個數列的初始化數列代入計算關聯度公式，易算出各數列的關聯度如下表：

關聯度表：

matlab程序如下：

數據：x.txt

13.6  14.01  14.54  15.64  15.69
 11.50  13.00  15.15  15.30  15.02
 13.76  16.36  16.90  16.56  17.30
 12.41  12.70  13.96  14.04  13.46
 2.48  2.49  2.56  2.64  2.59
 85  85  90  100  105
 55  65  75  80  80
 65  70  75  85  90
 12.80  15.30  16.24  16.40  17.05
 15.30  18.40  18.75  17.95  19.30
 12.71  14.50  14.66  15.88  15.70
 14.78  15.54  16.03  16.87  17.82
 7.64  7.56  7.76  7.54  7.70
 120  125  130  140  140
 80  85  90  90  95
 4.2  4.25  4.1  4.06  3.99
 13.1  13.42  12.85  12.72  12.56

程序：

% 關聯分析
load x.txt
 for i = 1:15
     x(i,:) = x(i,:)/x(i,1); %前15數列做標准化
end
 for i = 16:17
     x(i,:) = x(i,1)./x(i,:); %后兩個做標准化
end
 data = x;
 n = size(data,2); %矩陣列數，即觀測時刻的個數
ck = data(1,:); %選第一列是參考數列
bj = data(2:end,:);%其余列是比較數列
m2 = size(bj,1);%比較數列個數
for j = 1:m2
     t(j,:) = bj(j,:)-ck;
 end
 mn = min(min(abs(t'))); %最小差
mx = max(max(abs(t'))); %最大差
rho = 0.5; %分辨系數設置
ksi = (mn + rho*mx)./(abs(t)+rho*mx);%求關聯系數
r = sum(ksi')/n %關聯度
[rs,rind] = sort(r,'descend') %關聯度排序‘

matlab結果：

依次打印關聯度：

r =

  1 至 6 列

    0.5881    0.6627    0.8536    0.7763    0.8549    0.5022

  7 至 12 列

    0.6592    0.5820    0.6831    0.6958    0.8955    0.7047

  13 至 16 列

    0.9334    0.8467    0.7454    0.7261

排序后：

rs =

  1 至 6 列

    0.9334    0.8955    0.8549    0.8536    0.8467    0.7763

  7 至 12 列

    0.7454    0.7261    0.7047    0.6958    0.6831    0.6627

  13 至 16 列

    0.6592    0.5881    0.5820    0.5022

rind =

  1 至 10 列

    13    11     5     3    14     4    15    16    12    10

  11 至 16 列

     9     2     7     1     8     6

2.

優勢分析：

一個例子：假如有關聯度矩陣如下：

分析：

3.

生成數：

我們主要使用【累加生成】，其理論如下：

應用中，最常用的是 1 次累加生成。
一般地，經濟數列等實際問題的數列皆是非負數列，累加生成可使非負的擺動與非擺動的數列或任意無規律性的數列轉化為非減的數列。

有些實際問題的數列中有負數（例如溫度等），累加時略微復雜。有時，由於出現正負抵消這種信息損失的現象，數列經過累加生成后規律性非但沒得到加強，甚至可能被削弱。對於這種情形，我們可以先進行【移軸】，然后再做【累加生成】。

4.

灰色GM(1,1)模型

GM是Grey Model的簡寫。

1）GM(1,1)定義：

2）GM(1,1)的白化型：

應當注意,GM(1,1)表示模型師一階方程並且只有一個變量；推廣之，加入有m個方程，n個變量則為G(m,n)。pdf討論了G(1,N)/G(2,N)等，用到再說，對於這里，我們首先詳細使用G(1,1)。

5.

灰色預測

灰色預測是指利用 GM 模型對系統行為特征的發展變化規律進行估計【預測】，同時也可以對行為特征的異常情況發生的時刻進行估計計算(把異常時間作為數列)，以及對在特定時區內發生事件的未來時間分布情況做出研究等等。這些工作實質上是將“隨機過程”當作“灰色過程”，“隨機變量”當作“灰變量”，並主要以灰色系統理論中的 【GM(1,1)模型】來進行處理。

1) 灰色預測的方法

2) 灰色預測處理的步驟(使用)

1. 數據的檢驗和處理

2. 建立模型

參考上述GM(1,1)方法，建立GM(1,1)模型，並求解該微分方程，得到預測值：

3. 檢驗預測值

4. 預測預報

由模型 GM(1,1)所得到的指定時區內的預測值，根據實際問題的需要，給出相應的
預測預報。

另一個應用--災變預測

一個案例：

假定小於320為異常。預測下一次異常出現的時間(旱災)。

6.

灰色預測計算實例：如何使用、求解以及分析

灰色G(1,1)預測步驟：

步驟：

第一步：級比檢驗

(1)求級比，列出級比向量。

(2)級比判斷：若所有級比都落在可容覆蓋內，則通過，說明原始數據適合使用GM(1,1)

 

第二步：GM(1,1)建模

（1）原始數據做一次累加

（2）列出GM(1,1)模型

第三步：求解模型

（1）最小二乘法求解GM(1,1)

（2）求生成數列預測值以及模型還原值

（3）相應可以得到預測值-真實值比較表格

第四步：模型檢驗

（1）原始值-模型值-殘差-相對誤差0級比偏差表格

（2）根據表格作出說明

一個實例：

北方某城市 1986～1992 年道路交通噪聲平均聲級數據如下：

                                    

序號	年份	Leq
1	1986	71.1
2	1987	72.4
3	1988	72.4
4	1989	72.1
5	1990	71.4
6	1991	72.0
7	1992	71.6

求解分析：

求解MATLAB程序如下：

%此程序原著pdf上程序是有bug的，以下已經調通

clc,clear
 x0=[71.1 72.4 72.4 72.1 71.4 72.0 71.6]';%注意這里為列向量
n=length(x0);
 lamda=x0(1:n-1)./x0(2:n) %計算級比
range=minmax(lamda') %計算級比的范圍
x1=cumsum(x0); %累加運算
B=[-0.5*(x1(1:n-1)+x1(2:n)),ones(n-1,1)];
 Y=x0(2:n);
 u=B\Y %最小二乘，擬合參數
syms x(t)
 x=dsolve(diff(x)+u(1)*x==u(2),x(0)==x0(1));
 x=subs(x,{'a','b','x0'},{u(1),u(2),x1(1)});
 yuce1=subs(x,t,[0:n-1]);
 yuce1 = double(yuce1);
 %為提高預測精度，先計算預測值，再顯示微分方程的解
y=vpa(x,6) %其中的 6 表示顯示 6 位數字
yuce=[x0(1),diff(yuce1)] %差分運算，還原數據
epsilon=x0'-yuce %計算殘差
delta=abs(epsilon./x0') %計算相對誤差
rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda' %計算級比偏差值

求解結果如下：

級比：

lamda =

    0.9820
     1.0000
     1.0042
     1.0098
     0.9917
     1.0056

級比范圍：
range =

    0.9820    1.0098

模型參數：
u =

    0.0023
    72.6573

求解模型得方程：
y = 31000.0 - 30928.9*exp(-0.00234379*t)
 

模型值：

yuce =

  1 至 6 列

   71.1000   72.4057   72.2362   72.0671   71.8984   71.7301

  7 列

   71.5622

殘差：
epsilon =

  1 至 6 列

         0   -0.0057    0.1638    0.0329   -0.4984    0.2699

  7 列

    0.0378

相對誤差：
delta =

  1 至 6 列

         0    0.0001    0.0023    0.0005    0.0070    0.0037

  7 列

    0.0005

級比偏差：
rho =

    0.0203    0.0023   -0.0018   -0.0074    0.0107   -0.0032