灰色關聯度分析（Grey Relation Analysis，GRA）原理詳解

釋名

灰色關聯度分析（Grey Relation Analysis，GRA），是一種多因素統計分析的方法。簡單來講，就是在一個灰色系統中，我們想要了解其中某個我們所關注的某個項目受其他的因素影響的相對強弱，再直白一點，就是說：我們假設以及知道某一個指標可能是與其他的某幾個因素相關的，那么我們想知道這個指標與其他哪個因素相對來說更有關系，而哪個因素相對關系弱一點，依次類推，把這些因素排個序，得到一個分析結果，我們就可以知道我們關注的這個指標，與因素中的哪些更相關。

（ note ： 灰色系統這個概念的提出是相對於白色系統和黑色系統而言的。這個概念最初是由控制科學與工程（hhh熟悉的一級學科）的教授鄧聚龍提出的。按照控制論的慣例，顏色一般代表的是對於一個系統我們已知的信息的多少，白色就代表信息充足，比如一個力學系統，元素之間的關系都是能夠確定的，這就是一個白色系統；而黑色系統代表我們對於其中的結構並不清楚的系統，通常叫做黑箱或黑盒的就是這類系統。灰色介於兩者之間，表示我們只對該系統有部分了解。）

舉例

為了說明灰色關聯度分析的應用場景，我們利用下圖進行說明：

該圖來源於參考文獻1。這篇文獻研究的內容是旅游業發展的影響因子，看該表格，第一行為五年的旅游總收入，代表着旅游業發展的程度，而下面的這些要素就是我們需要分析的因子，比如在校大學生數，旅行社數，星級飯店數，A級景區數等等。最終目的是要得到一個排序，從而說明這些因子對旅游總收入的關聯性的程度。

操作步驟與原理詳解

（1） 確立母序列

（參考序列，在上面栗子中就是1998~2002年的旅游總收入序列）和子序列（比較序列，也就是需要確立順序的因素序列，上栗中的除了第一行以外的所有因素都可以作為參考序列）

為了后面的表述方便，這里統一一下notation：

我們用x_i(k)表示第i個因素的第k個數值
用上面的栗子來說
比如第一個因素是在校大學生人數，那么x_1(1)就表示在校大學生人數在1998年的取值，也就是341，x_1(2)就是1999年的取值，
而x_2(1)就是表示旅游從業人數在1998年的數值。以此類推。
我們用x_0(k)表示母序列，i≥1的表示子序列，也就是要分析的要素的序列。
如果不寫括號，比如x_i ，就代表這個元素的整個序列，也就是向量 x_i = [x_i(1), x_i(2), ... , x_i(n)] 
n為每個向量的維度，也就是每個元素的特征的數量，在上栗中，n就是5，因為有五年的數據，代表五維向量。
以下所有表述都用該notation表示。

詳解： 這個就是我們任務的目的（找到子序列和參考序列的關聯程度），所以不需要再解釋了吧～

（2）歸一化，或者叫 無量綱化。

詳解： 因為我們的這些要素是不同質的東西的指標，因此可能會有的數字很大有的數字很小，但是這並不是由於它們內稟的性質決定的，而只是由於量綱不同導致的，因此我們需要對它們進行無量綱化。這個操作一般在數據處理領域叫做歸一化（normalization），也就是減少數據的絕對數值的差異，將它們統一到近似的范圍內，然后重點關注其變化和趨勢。
如下圖所示，這是上面表格中前3個元素隨年的變化曲線，以及作為母序列的旅游總收入：

可以看到，有兩個曲線絕對數值很大，而另外兩個很小，如果不做處理必然導致大的數值的影響會”淹沒“掉小數值的變量的影響。

所以我們要對數據進行歸一化處理，主要方法有如下幾個：

（1） 初值化： 顧名思義，就是把這一個序列的數據統一除以最開始的值，由於同一個因素的序列的量級差別不大，所以通過除以初值就能將這些值都整理到1這個量級附近。

公式： x_i(k)' = x_i(k) / x_i(1)   i = 1,...,m, k = 1,...,n 
(m為因素個數，n為每個因素的數據維度，仍如上栗，n=5，m=3（我們只看前三個因素，就是曲線圖里畫的這三種，和旅游總收入的關聯，數據維度為5，即五年）)

（2） 均值化： 顧名思義，就是把這個序列的數據除以均值，由於數量級大的序列均值比較大，所以除掉以后就能歸一化到1的量級附近。

公式： x_i(k)' = x_i(k) / ( mean(x_i) )        （除以均值）
其中 ： mean(x_i) = (1/n) sum_k=1^n (x_i(k))     （求第i個因素序列的均值）

其余還有如區間化，即把序列的值規范到一個區間，比如[0, 1]，之間。這個方法實際上在數據處理中應用比較多，但是在GRA中似乎常用均值化或者初值化，所以在此不介紹。

這里我們按照參考文獻中的采用的方法，用初值化進行歸一化，得到的結果如下圖：

可以看到，歸一化以后的數據，量級差別變小了，這是為了后面提供鋪墊，因為我們關注的實際上是曲線的形狀的差異，而不希望絕對數值對后面的計算有影響。

（3）計算灰色關聯系數

先放上公式：

詳解：
首先，我們把i看做固定值，也就是說對於某一個因素，其中的每個維度進行計算，得到一個新的序列，這個序列中的每個點就代表着該子序列與母序列對應維度上的關聯性（數字越大，代表關聯性越強）。
仔細觀察這個公式，rho是一個可調節的系數，取值為(0,1)，大於零小於一，這一項的目的是為了調節輸出結果的差距大小，我們放在后面講。我們先假設把rho取成0，那么，這個式子就變成了

pseudo_zeta_i(k) = min min |x_0(k) - x_i(k)| / |x_0(k) - x_i(k)| = constant / |x_0(k) - x_i(k)|

我們看上面這個式子，可以發現，分子上這個數值，對於所有子序列來說都是一樣的（，分子上這個數實際上就是所有因素的所有維度中，與母序列（參考序列，即我們要比較的序列）距離最近的維度上的距離。為什么要這樣做呢？這樣來想，假如我們沒有進行歸一化，或者不是用的初值化，而是用的均值化或者其他方法，可能會導致曲線之間，也就是母序列和各個子序列之間仍然有一段距離，那么這個距離最小值與下面的每個維度的距離相除，實際上也可以看成是一種取消量綱的手段。對於所有子序列，這個分子是相同的，所以實際上，這個系數pseudo_zeta是與第k個維度上，子序列與母序列的距離（差的絕對值，通常叫做l1范數（l1-norm））成反比，也就是說，這兩個數距離越遠，我們認為越不相關，這是符合直覺的。

當然，如果用了初值化歸一化數據，如上面的圖2所示，min min |x_0(k) - x_i(k)| 對每個i都會變成0，這樣就不好了，因為這樣一來，所有的zeta_i(k)都成了0，是無意義的。所以這時候我們就看到后面的 rho max max這一項的作用了。這一項對於每個i來說也是一個不變的常數constant，所以可以理解為給上面那個式子的分子分母同時加上某個數值，如下所示：

zeta_i(k) = (aconstant + bconstant) / (|x_0(k) - x_i(k)| + bconstant)

這樣做的目的是什么呢？ 我們舉個栗子： 對於兩個分數： 1/5 和 1/4 ，它們的分子一樣，分母相差為1，這時候他們的值相差1/20，也就是0.05，這就是沒有＋rho max max那一項的情況，分子相同，分母的差代表着與參考序列的距離。 如果我們給他們分子分母同時加上20，那么就是21/25和21/24，它們相差為0.035，可以看到，加入這一項會導致同樣的距離的點的系數差，會因為計算而變小。很顯然地，rho取得越大，不同zeta系數的差距就越小。

另外，由於分子上是min min，也就是距離的全局最小值，這就導致下面的分母必然大於分子（不考慮 rho max max 項），而且，如果分母非常大，曲線距離非常遠，那么，zeta接近0； 相反，如果x_i和x_0在所有維度上的差完全一樣，那么分數的值就是1。這樣zeta取值范圍就是0~1之間，0表示不相關，1表示強關聯性。這也符合認知。考慮上rho max max 項之后，我們知道對於一個真分數，分子分母都加一個同樣的值，仍然是真分數（實際上是一個添加溶質的溶液的問題）。也就是說，仍然是0到1。

總結來說，rho是控制zeta系數區分度的一個系數，rho取值0到1， rho越小，區分度越大，一般取值0.5較為合適。zeta關聯系數取值落在0到1之間。

接上栗，我們對上述三個子序列做出關聯系數zeta的序列，結果如下：

其實從這個圖中已經可以看出，大學生這一因素對旅游也的相關性普遍要高一些，從業人員相對影響少一些。星級飯店的數量居中。

(4) 計算關聯系數均值，形成關聯序

根據上圖其實已經可以看出大概的趨勢，但是這只是因為這個恰好所有維度上的趨勢比較一致，實際上，我們得到zeta關聯系數的值以后，應該對每個因素在不同維度上的值求取均值，換句話說，也就是對於上面那些zeta 的曲線，同一個顏色的求取均值。結果如下：

>> mean(zeta_1)
ans =
   0.7505
>> mean(zeta_2)
ans =
    0.5848
>> mean(zeta_3)
ans =
   0.7154

可以看到，根據關聯系數大小，排序結果為：

大學生人數 ＞ 星級飯店數量 ＞ 從業人員人數

這和參考文獻論文中的結論是一致的：

（由於論文中用了所有的因子，導致max max 這個全局最大值不同，所以計算出的關聯度數值與復現計算的結果數值不一樣，但是這三個因子的關聯度的排序是一致的，說明關聯度是一個相對的指標，它反映的是不同因子與參考內容的關聯程度）

總結

GRA算法本質上來講就是提供了一種度量兩個向量之間距離的方法，對於有時間性的因子，向量可以看成一條時間曲線，而GRA算法就是度量兩條曲線的形態和走勢是否相近。為了避免其他干擾，凸出形態特征的影響，GRA先做了歸一化，將所有向量矯正到同一個尺度和位置，然后計算每個點的距離。最后，通過min min 和max max 的矯正，使得最終輸出的結果落在0到1之間，從而符合系數的一般定義。rho調節不同關聯系數之間的差異，換句話說，就是輸出的分布，使其可以變得更加稀疏或者緊密。以數學角度要言之，該算法即度量已歸一化的子向量與母向量的每一維度的l1-norm距離的倒數之和，並將其映射到0～1區間內，作為子母向量的關聯性之度量的一種策略。

附錄：MATLAB代碼

下面是該博文中為舉栗子復現的參考文獻1中的計算過程MATLAB代碼：

% Grey relation analysis

clear all
close all
clc

zongshouru = [3439, 4002, 4519, 4995, 5566];
daxuesheng = [341, 409, 556, 719, 903];
congyerenyuan = [183, 196, 564, 598, 613];
xingjifandian = [3248, 3856, 6029, 7358, 8880];

% define comparative and reference
x0 = zongshouru;
x1 = daxuesheng;
x2 = congyerenyuan;
x3 = xingjifandian;

% normalization
x0 = x0 ./ x0(1);
x1 = x1 ./ x1(1);
x2 = x2 ./ x2(1);
x3 = x3 ./ x3(1);

% global min and max
global_min = min(min(abs([x1; x2; x3] - repmat(x0, [3, 1]))));
global_max = max(max(abs([x1; x2; x3] - repmat(x0, [3, 1]))));

% set rho
rho = 0.5;

% calculate zeta relation coefficients
zeta_1 = (global_min + rho * global_max) ./ (abs(x0 - x1) + rho * global_max);
zeta_2 = (global_min + rho * global_max) ./ (abs(x0 - x2) + rho * global_max);
zeta_3 = (global_min + rho * global_max) ./ (abs(x0 - x3) + rho * global_max);

% show
figure;
plot(x0, 'ko-' )
hold on
plot(x1, 'b*-')
hold on
plot(x2, 'g*-')
hold on
plot(x3, 'r*-')
legend('zongshouru', 'daxuesheng', 'congyerenyuan', 'xingjifandian')

figure;
plot(zeta_1, 'b*-')
hold on
plot(zeta_2, 'g*-')
hold on
plot(zeta_3, 'r*-')
title('Relation zeta')
legend('daxuesheng', 'congyerenyuan', 'xingjifandian')

參考文獻：

1. 馬曉龍. 旅游業發展影響因子灰色關聯分析[J]. 人文地理, 2006, 21(2):37-40.
2. 譚學瑞, 鄧聚龍. 灰色關聯分析：多因素統計分析新方法[J]. 統計研究, 1995, 12(3):46-48.
3. 劉思峰, 蔡華, 楊英傑, et al. 灰色關聯分析模型研究進展[J]. 系統工程理論與實踐, 2013, 33(8):2041-2046.

2018年12月29日02:06:16
to Rita ～