本文介紹了幾個重要的變數相關性的度量，包括皮爾遜相關係數、距離相關性和最大資訊係數等，並用簡單的程式碼和示例資料展示了這些度量的適用性對比。
從訊號的角度來看，這個世界是一個嘈雜的地方。為了弄清楚所有的事情，我們必須有選擇地集中注意力到有用的資訊上。
透過數百萬年的自然選擇過程，我們人類已經變得非常擅長過濾背景訊號。我們學會將特定的訊號與特定的事件聯絡起來。
例如，假設你正在繁忙的辦公室中打乒乓球。
為了回擊對手的擊球，你需要進行大量複雜的計算和判斷，將多個相互競爭的感官訊號考慮進去。
為了預測球的運動，你的大腦必須重複取樣球的位置並估計它未來的軌跡。更厲害的球員還會將對手擊球時施加的旋轉考慮進去。
最後，為了擊球，你需要考慮對手的位置、自己的位置、球的速度，以及你打算施加的旋轉。
所有這些都涉及到了大量的潛意識微分學。一般來說，我們理所當然的認為，我們的神經系統可以自動做到這些（至少經過一些練習之後）。
同樣令人印象深刻的是，人類大腦是如何區別對待它所接收到的無數競爭訊號的重要性的。例如，球的位置被認為比你身後發生的對話或你面前開啟的門更重要。
這聽起來似乎不值得一提，但實際上這證明了可以多大程度上學習從噪聲資料中做出準確預測。
當然，一個被給予連續的視聽資料流的空白狀態機將會面臨一個困難的任務，即確定哪些訊號能夠最好地預測最佳行動方案。
幸運的是，有統計和計算方法可以用來識別帶噪聲和複雜的資料中的模式。
相關性
一般來說，當我們談到兩個變數之間的「相關性（correlation）」時，在某種意義上，我們是指它們的「關係（relatedness）」。
相關變數是包含彼此資訊的變數。兩個變數的相關性越強，其中一個變數告訴我們的關於另一個變數的資訊就越多。
你可能之前就看過：正相關、零相關、負相關
你可能已經對相關性、它的作用和它的侷限性有了一定了解。事實上，這是一個數據科學的老生常談：
「相關性不意味著因果關係」
這當然是正確的——有充分的理由說明，即使是兩個變數之間有強相關性也不保證存在因果關係。觀察到的相關性可能是由於隱藏的第三個變數的影響，或者完全是偶然的。
也就是說，相關性確實允許基於另一個變數來預測一個變數。有幾種方法可以用來估計線性和非線性資料的相關性。我們來看看它們是如何工作的。
我們將用 Python 和 R 來進行數學和程式碼實現。本文示例的程式碼可以在這裡 (https://gist.github.com/anonymous/fabecccf33f9c3feb568384f626a2c07) 找到：
GitHub 地址：https://gist.github.com/anonymous/fabecccf33f9c3feb568384f626a2c07
皮爾遜相關係數
皮爾遜相關係數是一種廣泛使用的線性相關性的度量，它通常是很多初級統計課程的第一課。從數學角度講，它被定義為「兩個向量之間的協方差，透過它們標準差的乘積來歸一化」。
兩個成對的向量之間的協方差是它們在均值上下波動趨勢的一種度量。也就是說，衡量一對向量是否傾向於在各自平均值的同側或相反。
讓我們看看在 Python 中的實現：
協方差的計算方法是從每一對變數中減去各自的均值。然後，將這兩個值相乘。
如果都高於（或都低於）均值，那麼結果將是一個正數，因為正數 × 正數 = 正數；同樣的，負數 × 負數 = 負數。
如果在均值的不同側，那麼結果將是一個負數（因為正數 × 負數 = 負數）。
一旦我們為每一對變數都計算出這些值，將它們加在一起，併除以 n-1，其中 n 是樣本大小。這就是樣本協方差。
如果這些變數都傾向於分佈在各自均值的同一側，協方差將是一個正數；反之，協方差將是一個負數。這種傾向越強，協方差的絕對值就越大。
如果不存在整體模式，那麼協方差將會接近於零。這是因為正值和負值會相互抵消。
最初，協方差似乎是兩個變數之間「關係」的充分度量。但是，請看下面的圖：
協方差 = 0.00003
看起來變數之間有很強的關係，對吧？那為什麼協方差這麼小呢（大約是 0.00003）？
這裡的關鍵是要認識到協方差是依賴於比例的。看一下 x 和 y 座標軸——幾乎所有的資料點都落在了 0.015 和 0.04 之間。協方差也將接近於零，因為它是透過從每個個體觀察值中減去平均值來計算的。
為了獲得更有意義的數字，歸一化協方差是非常重要的。方法是將其除以兩個向量標準差的乘積。
在希臘字母中ρ經常用來表示皮爾遜相關係數
在 Python 中：
這樣做的原因是因為向量的標準差是是其方差的平方根。這意味著如果兩個向量是相同的，那麼將它們的標準差相乘就等於它們的方差。
有趣的是，兩個相同向量的協方差也等於它們的方差。
因此，兩個向量之間協方差的最大值等於它們標準差的乘積（當向量完全相關時會出現這種情況）。這將相關係數限制在 -1 到 +1 之間。
箭頭指向哪個方向？
順便說一下，一個定義兩個向量的 PCC 的更酷的方法來自線性代數。
首先，我們透過從向量各自的值中減去其均值的方法來「集中」向量。
a = [1,2,3,4,5] ; b = [5,4,3,2,1]
a_centered = [i - mean(a) for i in a]
b_centered = [j - mean(b) for j in b]
現在，我們可以利用向量可以看做指向特定方向的「箭頭」的事實。
例如，在 2-D 空間中，向量 [1,3] 可以代表一個沿 x 軸 1 個單位，沿 y 軸 3 個單位的箭頭。同樣，向量 [2,1] 可以代表一個沿 x 軸 2 個單位，沿 y 軸 1 個單位的箭頭。
兩個向量 (1,3) 和 (2,1) 如箭頭所示。
類似的，我們可以將資料向量表示為 n 維空間中的箭頭（儘管當 n > 3 時不能嘗試視覺化）。
這些箭頭之間的角度 ϴ 可以使用兩個向量的點積來計算。定義為：
或者，在 Python 中：
點積也可以被定義為：
其中 ||**x**|| 是向量 **x **的大小（或「長度」）（參考勾股定理），ϴ 是箭頭向量之間的角度。
正如一個 Python 函式：
def magnitude(x):
x_sq = [i ** 2 for i in x]
return math.sqrt(sum(x_sq))
我們透過將點積除以兩個向量大小的乘積的方法得到 cos(ϴ)。
def cosTheta(x,y):
mag_x = magnitude(x)
mag_y = magnitude(y)
return dotProduct(x,y) / (mag_x * mag_y)
現在，如果你對三角學有一定了解，你可能會記得，餘弦函式產生一個在 +1 和 -1 之間震盪的圖形。
cos(ϴ) 的值將根據兩個箭頭向量之間的角度而發生變化。
當角度為零時（即兩個向量指向完全相同的方向），cos(ϴ) 等於 1。
當角度為 -180°時（兩個向量指向完全相反的方向），cos(ϴ) 等於 -1。
當角度為 90°時（兩個向量指向完全不相關的方向），cos(ϴ) 等於零。
這可能看起來很熟悉——一個介於 +1 和 -1 之間的衡量標準似乎描述了兩個向量之間的關係？那不是 Pearson』s *r *嗎？
那麼——這正是它的解釋！透過將資料視為高維空間中的箭頭向量，我們可以用它們之間的角度 ϴ 作為相似度的衡量。
A) 正相關向量； B) 負相關向量； C) 不相關向量
該角度 ϴ 的餘弦在數學上與皮爾遜相關係數相等。當被視為高維箭頭時，正相關向量將指向一個相似的方向。負相關向量將指向相反的方向。而不相關向量將指向直角。
就我個人而言，我認為這是一個理解相關性的非常直觀的方法。
統計顯著性？
正如頻率統計一樣，重要的是詢問從給定樣本計算的檢驗統計量實際上有多重要。Pearson』s r* *也不例外。
不幸的是，PCC 估計的置信區間不是完全直接的。
這是因為 Pearson』s r 被限制在 -1 和 +1 之間，因此不是正態分佈的。而估計 PCC，例如 +0.95 之上只有很少的容錯空間，但在其之下有大量的容錯空間。
幸運的是，有一個解決方案——用一個被稱為 Fisher 的 Z 變換的技巧：
1. 像平常一樣計算 Pearson』s r 的估計值。
2. 用 Fisher 的 Z 變換將 r→z，用公式 z = arctanh(r) 完成。
3. 現在計算 z 的標準差。幸運的是，這很容易計算，由 SDz = 1/sqrt(n-3) 給出，其中 n 是樣本大小。
4. 選擇顯著性閾值，alpha，並檢查與此對應的平均值有多少標準差。如果取 alpha = 0.95，用 1.96。
5. 透過計算 *z* +(1.96 × SD*z*) 找到上限，透過計算 *z -* **(1.96 × SD*z*) 找到下限。
6. 用 r = tanh(z) 將這些轉換回 r。
7. 如果上限和下限都在零的同一側，則有統計顯著性！
這裡是在 Python 中的實現：
r = Pearsons(x,y)
z = math.atanh(r)
SD_z = 1 / math.sqrt(len(x) - 3)
z_upper = z + 1.96 * SD_z
z_lower = z - 1.96 * SD_z
r_upper = math.tanh(z_upper)
r_lower = math.tanh(z_lower)
當然，當給定一個包含許多潛在相關變數的大資料集時，檢查每對的相關性可能很吸引人。這通常被稱為「資料疏浚」——在資料集中查詢變數之間的任何明顯關係。
如果確實採用這種多重比較方法，則應該用適當的更嚴格的顯著性閾值來降低發現錯誤相關性的風險（即找到純粹偶然相關的無關變數）。
一種方法是使用 Bonferroni correction。
小結
到現在為止還好。我們已經看到 Pearson』s *r* 如何用來計算兩個變數之間的相關係數，以及如何評估結果的統計顯著性。給定一組未知的資料，用於開始挖掘變數之間的重要關係是很有可能的。
但是，有一個重要的陷阱——Pearson』s r 只適用於線性資料。
看下面的圖。它們清楚地展示了一種看似非隨機的關係，但是 Pearson』s *r *非常接近於零。
原因是因為這些圖中的變數具有非線性關係。
我們通常可以將兩個變數之間的關係描繪成一個點雲，分散在一條線的兩側。點雲的分散度越大，資料越「嘈雜」，關係越弱。
然而，由於它將每個單獨的資料點與整體平均值進行比較，所以 Pearson』s r 只考慮直線。這意味著檢測非線性關係並不是很好。
在上面的圖中，Pearson』s r 並沒有顯示研究物件的相關性。
然而，這些變數之間的關係很顯然是非隨機的。幸運的是，我們有不同的相關性方法。
讓我們來看看其中幾個。
距離相關性
距離相關性與 Pearson』s *r *有一些相似之處，但是實際上是用一個相當不同的協方差概念來計算的。該方法透過用「距離」類似物替代常用的協方差和標準差（如上所定義）的概念。
類似 Pearson』s r，「距離相關性」被定義為「距離協方差」，由「距離標準差」來歸一化。
距離相關性不是根據它們與各自平均值的距離來估計兩個變數如何共同變化，而是根據與其他點的距離來估計它們是如何共同變化的，從而能更好捕捉變數之間非線性依賴關係。
深入細節
出生於 1773 年的 Robert Brown 是一名蘇格蘭植物學家。當布朗在顯微鏡下研究植物花粉時，注意到液麵上有隨機運動的有機顆粒。
他沒有想到，這一觀察竟使他名垂千古——他成為了布朗運動的（重新）發現者。
他更不會知道，近一個世紀的時間後愛因斯坦才對這種現象做出瞭解釋，從而證實了原子的存在。同年，愛因斯坦發表了關於狹義相對論的論文（E=MC²），並打開了量子理論的大門。
布朗運動是這樣一個物理過程：由於與周圍粒子的碰撞，微小粒子隨機運動。
布朗運動背後的數學原理可以被推廣為維納過程（Weiner process），維納過程在數學金融中最著名的模型 Black-Scholes 中也扮演著重要的角色。
有趣的是，Gabor Szekely 在 20 世紀中期的研究表明，布朗運動和維納過程和一個非線性關聯度量相關。
讓我們來看看如何由長度為 N 的向量 x 和 y 計算這個量。
1. 首先，我們對每個向量構建 N×N 的距離矩陣。距離矩陣和地圖中的道路距離表非常類似——每行、每列的交點顯示了相應城市間的距離。在距離矩陣中，行 i 和列 j 的交點給出了向量的第 i 個元素和第 j 個元素之間的距離。
2. 第二，矩陣是「雙中心」的。也就是說，對於每個元素，我們減去了它的行平均值和列平均值。然後，我們再加上整個矩陣的總平均值。
上述公式中，加「^」表示「雙中心」，加「-」表示「平均值」。
3. 在兩個雙中心矩陣的基礎上，將 X 中每個元素的均值乘以 Y 中相應元素的均值，則可計算出距離協方差的平方。
4. 現在，我們可以用類似的辦法找到「距離方差」。請記住，若兩個向量相同，其協方差與其方差相等。因此，距離方差可表示如下：
5. 最後，我們利用上述公式計算距離相關性。請記住，（距離）標準差與（距離）方差的平方根相等。
如果你更喜歡程式碼實現而非數學符號，那麼請看下面的 R 語言實現：
任意兩變數的距離相關性都在 0 和 1 之間。其中，0 代表兩變數相互獨立，而接近於 1 則表明變數間存在依賴關係。
如果你不想從頭開始編寫距離相關方法，你可以安裝 R 語言的 energy 包（https://cran.r-project.org/web/packages/energy/index.html），設計此方案的研究者提供了本程式碼。在該程式包中，各類可用方案呼叫的是 C 語言編寫的函式，因此有著很大的速度優勢。
置信區間？
我們可以採取「重取樣（resampling）」方法為距離相關性估計建立置信區間。一個簡單的例子是 bootstrap 重取樣。
這是一個巧妙的統計技巧，需要我們從原始資料集中隨機抽樣（替換）以「重建」資料。這個過程將重複多次（例如 1000 次），每次都計算感興趣的統計量。
這將為我們感興趣的統計量產生一系列不同的估計值。我們可以透過它們估計在給定置信水平下的上限和下限。
請看下面的 R 語言程式碼，它實現了簡單的 bootstrap 函式：
如果你想建立統計顯著性，還有另一個重取樣技巧，名為「排列檢驗（permutation test）」。
排列檢驗與上述 bootstrap 方法略有不同。在排列檢驗中，我們保持一個向量不變，並透過重取樣對另一個變數進行「洗牌」。這接近於零假設（null hypothesis）——即，在變數之間不存在依賴關係。
這個經「洗牌」打亂的變數將被用於計算它和常變數間的距離相關性。這個過程將被執行多次，然後，結果的分佈將與實際距離相關性（從未被「洗牌」的資料中獲得）相比較。
然後，大於或等於「實際」結果的經「洗牌」的結果的比例將被定為 P 值，並與給定的顯著性閾值（如 0.05）進行比較。
以下是上述過程的程式碼實現：
最大資訊係數
最大資訊係數（MIC）於 2011 年提出，它是用於檢測變數之間非線性相關性的最新方法。用於進行 MIC 計算的演算法將資訊理論和機率的概念應用於連續型資料。
由克勞德·夏農於 20 世紀中葉開創的資訊理論是數學中一個引人注目的領域。
資訊理論中的一個關鍵概念是熵——這是一個衡量給定機率分佈的不確定性的度量。機率分佈描述了與特定事件相關的一系列給定結果的機率。
機率分佈的熵是「每個可能結果的機率乘以其對數後的和」的負值。
為了理解其工作原理，讓我們比較下面兩個機率分佈：
X 軸標明瞭可能的結果；Y 軸標明瞭它們各自的機率
左側是一個常規六面骰子結果的機率分佈；而右邊的六面骰子不那麼均勻。
從直覺上來說，你認為哪個的熵更高呢？哪個骰子結果的不確定性更大？讓我們來計算它們的熵，看看答案是什麼。
不出所料，常規骰子的熵更高。這是因為每種結果的可能性都一樣，所以我們不會提前知道結果偏向哪個。但是，非常規的骰子有所不同——某些結果的發生機率遠大於其它結果——所以它的結果的不確定性也低一些。
這麼一來，我們就能明白，當每種結果的發生機率相同時，它的熵最高。而這種機率分佈也就是傳說中的「均勻」分佈。
交叉熵是熵的一個拓展概念，它引入了第二個變數的機率分佈。
crossEntropy <- function(x,y){
prX <- prop.table(table(x))
prY <- prop.table(table(y))
H <- sum(prX * log(prY,2))
return(-H)
}
兩個相同機率分佈之間的交叉熵等於其各自單獨的熵。但是對於兩個不同的機率分佈，它們的交叉熵可能跟各自單獨的熵有所不同。
這種差異，或者叫「散度」可以透過 KL 散度（Kullback-Leibler divergence）量化得出。
兩機率分佈 X 與 Y 的 KL 散度如下：
機率分佈 X 與 Y 的 KL 散度等於它們的交叉熵減去 X 的熵。
KL 散度的最小值為 0，僅當兩個分佈相同。
KL_divergence <- function(x,y){
kl <- crossEntropy(x,y) - entropy(x)
return(kl)
為了發現變數具有相關性，KL 散度的用途之一是計算兩個變數的互資訊（MI）。
互資訊可以定義為「兩個隨機變數的聯合分佈和邊緣分佈之間的 KL 散度」。如果二者相同，MI 值取 0。如若不同，MI 值就為一個正數。二者之間的差異越大，MI 值就越大。
為了加深理解，我們首先簡單回顧一些機率論的知識。
變數 X 和 Y 的聯合機率就是二者同時發生的機率。例如，如果你拋擲兩枚硬幣 X 和 Y，它們的聯合分佈將反映拋擲結果的機率。假設你拋擲硬幣 100 次，得到「正面、正面」的結果 40 次。聯合分佈將反映如下。
P(X=H, Y=H) = 40/100 = 0.4
jointDist <- function(x,y){
N <- length(x)
u <- unique(append(x,y))
joint <- c()
for(i in u){
for(j in u){
f <- x[paste0(x,y) == paste0(i,j)]
joint <- append(joint, length(f)/N)
return(joint)
邊緣分佈是指不考慮其它變數而只關注某一特定變數的機率分佈。假設兩變數獨立，二者邊緣機率的乘積即為二者同時發生的機率。仍以拋硬幣為例，假如拋擲結果是 50 次正面和 50 次反面，它們的邊緣分佈如下：
P(X=H) = 50/100 = 0.5 ; P(Y=H) = 50/100 = 0.5
P(X=H) × P(Y=H) = 0.5 × 0.5 = 0.25
marginalProduct <- function(x,y){
marginal <- c()
fX <- length(x[x == i]) / N
fY <- length(y[y == j]) / N
marginal <- append(marginal, fX * fY)
return(marginal)
現在讓我們回到拋硬幣的例子。如果兩枚硬幣相互獨立，邊緣分佈的乘積表示每個結果可能發生的機率，而聯合分佈則為實際得到的結果的機率。
如果兩硬幣完全獨立，它們的聯合機率在數值上（約）等於邊緣分佈的乘積。若只是部分獨立，此處就存在散度。
這個例子中，P(X=H,Y=H) > P(X=H) × P(Y=H)。這表明兩硬幣全為正面的機率要大於它們的邊緣分佈之積。
聯合分佈和邊緣分佈乘積之間的散度越大，兩個變數之間相關的可能性就越大。兩個變數的互資訊定義了散度的度量方式。
X 和 Y 的互資訊等於「二者邊緣分佈積和的聯合分佈的 KL 散度」
mutualInfo <- function(x,y){
joint <- jointDist(x,y)
marginal <- marginalProduct(x,y)
Hjm <- - sum(joint[marginal > 0] * log(marginal[marginal > 0],2))
Hj <- - sum(joint[joint > 0] * log(joint[joint > 0],2))
return(Hjm - Hj)
此處的一個重要假設就是機率分佈是離散的。那麼我們如何把這些概念應用到連續的機率分佈呢？
分箱演算法
其中一種方法是量化資料（使變數離散化）。這是透過分箱演算法（bining）實現的，它能將連續的資料點分配對應的離散類別。
此方法的關鍵問題是到底要使用多少「箱子（bin）」。幸運的是，首次提出 MIC 的論文給出了建議：窮舉！
也就是說，去嘗試不同的「箱子」個數並觀測哪個會在變數間取到最大的互資訊值。不過，這提出了兩個挑戰：
1. 要試多少個箱子呢？理論上你可以將變數量化到任意間距值，可以使箱子尺寸越來越小。
2. 互資訊對所用的箱子數很敏感。你如何公平比較不同箱子數目之間的 MI 值？
第一個挑戰從理論上講是不能做到的。但是，論文作者提供了一個啟發式解法（也就是說，解法不完美，但是十分接近完美解法）。他們也給出了可試箱子個數的上限。
最大可用箱子個數由樣本數 N 決定
至於如何公平比較取不同箱子數對 MI 值的影響，有一個簡單的做法……就是歸一化！這可以透過將每個 MI 值除以在特定箱子數組合上取得的理論最大值來完成。我們要採用的是產生最大歸一化 MI 總值的箱子數組合。
互資訊可以透過除以最小的箱子數的對數來歸一化
最大的歸一化互資訊就是 X 和 Y 的最大資訊係數（MIC）。我們來看看一些估算兩個連續變數的 MIC 的程式碼。
MIC <- function(x,y){
maxBins <- ceiling(N ** 0.6)
MI <- c()
for(i in 2:maxBins) {
for (j in 2:maxBins){
if(i * j > maxBins){
next
Xbins <- i; Ybins <- j
binnedX <-cut(x, breaks=Xbins, labels = 1:Xbins)
binnedY <-cut(y, breaks=Ybins, labels = 1:Ybins)
MI_estimate <- mutualInfo(binnedX,binnedY)
MI_normalized <- MI_estimate / log(min(Xbins,Ybins),2)
MI <- append(MI, MI_normalized)
return(max(MI))
x <- runif(100,-10,10)
y <- x**2 + rnorm(100,0,10)
MIC(x,y) # --> 0.751
以上程式碼是對原論文中方法的簡化。更接近原作的演算法實現可以參考 R package *minerva*（https://cran.r-project.org/web/packages/minerva/index.html）。
在 Python 中的實現請參考 minepy module（https://minepy.readthedocs.io/en/latest/）。
MIC 能夠表示各種線性和非線性的關係，並已得到廣泛應用。它的值域在 0 和 1 之間，值越高表示相關性越強。
為了建立 MIC 估計值的置信區間，你可以簡單地使用一個像我們之前介紹過的 bootstrap 函式。我們可以利用 R 語言的函數語言程式設計，透過傳遞我們想要用作引數的函式來泛化 bootstrap 函式。
bootstrap <- function(x,y,func,reps,alpha){
estimates <- c()
original <- data.frame(x,y)
N <- dim(original)[1]
for(i in 1:reps){
S <- original[sample(1:N, N, replace = TRUE),]
estimates <- append(estimates, func(S$x, S$y))
l <- alpha/2 ; u <- 1 - l
interval <- quantile(estimates, c(u, l))
return(2*(func(x,y)) - as.numeric(interval[1:2]))
bootstrap(x,y,MIC,100,0.05) # --> 0.594 to 0.88
為了總結相關性這一主題，我們來測試下各演算法在人工生成資料上的處理能力。
完整程式碼：https://gist.github.com/anonymous/fabecccf33f9c3feb568384f626a2c07
噪聲函式
set.seed(123)
# Noise
x0 <- rnorm(100,0,1)
y0 <- rnorm(100,0,1)
plot(y0~x0, pch = 18)
cor(x0,y0)
distanceCorrelation(x0,y0)
MIC(x0,y0)
皮爾遜相關係數 r = - 0.05
距離相關性 = 0.157
MIC = 0.097
簡單線性函式
# Simple linear relationship
x1 <- -20:20
y1 <- x1 + rnorm(41,0,4)
plot(y1~x1, pch =18)
cor(x1,y1)
distanceCorrelation(x1,y1)
MIC(x1,y1)
皮爾遜相關係數 r =+0.95
距離相關性 = 0.95
MIC = 0.89
簡單二次函式
# y ~ x**2
x2 <- -20:20
y2 <- x2**2 + rnorm(41,0,40)
plot(y2~x2, pch = 18)
cor(x2,y2)
distanceCorrelation(x2,y2)
MIC(x2,y2)
Pearson』s r =+0.003
距離相關性 = 0.474
MIC = 0.594
三次函式
# Cosine
x3 <- -20:20
y3 <- cos(x3/4) + rnorm(41,0,0.2)
plot(y3~x3, type="p", pch=18)
cor(x3,y3)
distanceCorrelation(x3,y3)
MIC(x3,y3)
皮爾遜相關係數 r =- 0.035
距離相關性 = 0.382
MIC = 0.484
圓函式
# Circle
n <- 50
theta <- runif (n, 0, 2*pi)
x4 <- append(cos(theta), cos(theta))
y4 <- append(sin(theta), -sin(theta))
plot(x4,y4, pch=18)
cor(x4,y4)
distanceCorrelation(x4,y4)
MIC(x4,y4)
皮爾遜相關係數 r < 0.001
距離相關性 = 0.234
MIC = 0.218
本文介紹了幾個重要的變數相關性的度量，包括皮爾遜相關係數、距離相關性和最大資訊係數等，並用簡單的程式碼和示例資料展示了這些度量的適用性對比。
從訊號的角度來看，這個世界是一個嘈雜的地方。為了弄清楚所有的事情，我們必須有選擇地集中注意力到有用的資訊上。
透過數百萬年的自然選擇過程，我們人類已經變得非常擅長過濾背景訊號。我們學會將特定的訊號與特定的事件聯絡起來。
例如，假設你正在繁忙的辦公室中打乒乓球。
為了回擊對手的擊球，你需要進行大量複雜的計算和判斷，將多個相互競爭的感官訊號考慮進去。
為了預測球的運動，你的大腦必須重複取樣球的位置並估計它未來的軌跡。更厲害的球員還會將對手擊球時施加的旋轉考慮進去。
最後，為了擊球，你需要考慮對手的位置、自己的位置、球的速度，以及你打算施加的旋轉。
所有這些都涉及到了大量的潛意識微分學。一般來說，我們理所當然的認為，我們的神經系統可以自動做到這些（至少經過一些練習之後）。
同樣令人印象深刻的是，人類大腦是如何區別對待它所接收到的無數競爭訊號的重要性的。例如，球的位置被認為比你身後發生的對話或你面前開啟的門更重要。
這聽起來似乎不值得一提，但實際上這證明了可以多大程度上學習從噪聲資料中做出準確預測。
當然，一個被給予連續的視聽資料流的空白狀態機將會面臨一個困難的任務，即確定哪些訊號能夠最好地預測最佳行動方案。
幸運的是，有統計和計算方法可以用來識別帶噪聲和複雜的資料中的模式。
相關性
一般來說，當我們談到兩個變數之間的「相關性（correlation）」時，在某種意義上，我們是指它們的「關係（relatedness）」。
相關變數是包含彼此資訊的變數。兩個變數的相關性越強，其中一個變數告訴我們的關於另一個變數的資訊就越多。
你可能之前就看過：正相關、零相關、負相關
你可能已經對相關性、它的作用和它的侷限性有了一定了解。事實上，這是一個數據科學的老生常談：
「相關性不意味著因果關係」
這當然是正確的——有充分的理由說明，即使是兩個變數之間有強相關性也不保證存在因果關係。觀察到的相關性可能是由於隱藏的第三個變數的影響，或者完全是偶然的。
也就是說，相關性確實允許基於另一個變數來預測一個變數。有幾種方法可以用來估計線性和非線性資料的相關性。我們來看看它們是如何工作的。
我們將用 Python 和 R 來進行數學和程式碼實現。本文示例的程式碼可以在這裡 (https://gist.github.com/anonymous/fabecccf33f9c3feb568384f626a2c07) 找到：
GitHub 地址：https://gist.github.com/anonymous/fabecccf33f9c3feb568384f626a2c07
皮爾遜相關係數
皮爾遜相關係數是一種廣泛使用的線性相關性的度量，它通常是很多初級統計課程的第一課。從數學角度講，它被定義為「兩個向量之間的協方差，透過它們標準差的乘積來歸一化」。
兩個成對的向量之間的協方差是它們在均值上下波動趨勢的一種度量。也就是說，衡量一對向量是否傾向於在各自平均值的同側或相反。
讓我們看看在 Python 中的實現：
協方差的計算方法是從每一對變數中減去各自的均值。然後，將這兩個值相乘。
如果都高於（或都低於）均值，那麼結果將是一個正數，因為正數 × 正數 = 正數；同樣的，負數 × 負數 = 負數。
如果在均值的不同側，那麼結果將是一個負數（因為正數 × 負數 = 負數）。
一旦我們為每一對變數都計算出這些值，將它們加在一起，併除以 n-1，其中 n 是樣本大小。這就是樣本協方差。
如果這些變數都傾向於分佈在各自均值的同一側，協方差將是一個正數；反之，協方差將是一個負數。這種傾向越強，協方差的絕對值就越大。
如果不存在整體模式，那麼協方差將會接近於零。這是因為正值和負值會相互抵消。
最初，協方差似乎是兩個變數之間「關係」的充分度量。但是，請看下面的圖：
協方差 = 0.00003
看起來變數之間有很強的關係，對吧？那為什麼協方差這麼小呢（大約是 0.00003）？
這裡的關鍵是要認識到協方差是依賴於比例的。看一下 x 和 y 座標軸——幾乎所有的資料點都落在了 0.015 和 0.04 之間。協方差也將接近於零，因為它是透過從每個個體觀察值中減去平均值來計算的。
為了獲得更有意義的數字，歸一化協方差是非常重要的。方法是將其除以兩個向量標準差的乘積。
在希臘字母中ρ經常用來表示皮爾遜相關係數
在 Python 中：
這樣做的原因是因為向量的標準差是是其方差的平方根。這意味著如果兩個向量是相同的，那麼將它們的標準差相乘就等於它們的方差。
有趣的是，兩個相同向量的協方差也等於它們的方差。
因此，兩個向量之間協方差的最大值等於它們標準差的乘積（當向量完全相關時會出現這種情況）。這將相關係數限制在 -1 到 +1 之間。
箭頭指向哪個方向？
順便說一下，一個定義兩個向量的 PCC 的更酷的方法來自線性代數。
首先，我們透過從向量各自的值中減去其均值的方法來「集中」向量。
a = [1,2,3,4,5] ; b = [5,4,3,2,1]
a_centered = [i - mean(a) for i in a]
b_centered = [j - mean(b) for j in b]
現在，我們可以利用向量可以看做指向特定方向的「箭頭」的事實。
例如，在 2-D 空間中，向量 [1,3] 可以代表一個沿 x 軸 1 個單位，沿 y 軸 3 個單位的箭頭。同樣，向量 [2,1] 可以代表一個沿 x 軸 2 個單位，沿 y 軸 1 個單位的箭頭。
兩個向量 (1,3) 和 (2,1) 如箭頭所示。
類似的，我們可以將資料向量表示為 n 維空間中的箭頭（儘管當 n > 3 時不能嘗試視覺化）。
這些箭頭之間的角度 ϴ 可以使用兩個向量的點積來計算。定義為：
或者，在 Python 中：
點積也可以被定義為：
其中 ||**x**|| 是向量 **x **的大小（或「長度」）（參考勾股定理），ϴ 是箭頭向量之間的角度。
正如一個 Python 函式：
def magnitude(x):
x_sq = [i ** 2 for i in x]
return math.sqrt(sum(x_sq))
我們透過將點積除以兩個向量大小的乘積的方法得到 cos(ϴ)。
def cosTheta(x,y):
mag_x = magnitude(x)
mag_y = magnitude(y)
return dotProduct(x,y) / (mag_x * mag_y)
現在，如果你對三角學有一定了解，你可能會記得，餘弦函式產生一個在 +1 和 -1 之間震盪的圖形。
cos(ϴ) 的值將根據兩個箭頭向量之間的角度而發生變化。
當角度為零時（即兩個向量指向完全相同的方向），cos(ϴ) 等於 1。
當角度為 -180°時（兩個向量指向完全相反的方向），cos(ϴ) 等於 -1。
當角度為 90°時（兩個向量指向完全不相關的方向），cos(ϴ) 等於零。
這可能看起來很熟悉——一個介於 +1 和 -1 之間的衡量標準似乎描述了兩個向量之間的關係？那不是 Pearson』s *r *嗎？
那麼——這正是它的解釋！透過將資料視為高維空間中的箭頭向量，我們可以用它們之間的角度 ϴ 作為相似度的衡量。
A) 正相關向量； B) 負相關向量； C) 不相關向量
該角度 ϴ 的餘弦在數學上與皮爾遜相關係數相等。當被視為高維箭頭時，正相關向量將指向一個相似的方向。負相關向量將指向相反的方向。而不相關向量將指向直角。
就我個人而言，我認為這是一個理解相關性的非常直觀的方法。
統計顯著性？
正如頻率統計一樣，重要的是詢問從給定樣本計算的檢驗統計量實際上有多重要。Pearson』s r* *也不例外。
不幸的是，PCC 估計的置信區間不是完全直接的。
這是因為 Pearson』s r 被限制在 -1 和 +1 之間，因此不是正態分佈的。而估計 PCC，例如 +0.95 之上只有很少的容錯空間，但在其之下有大量的容錯空間。
幸運的是，有一個解決方案——用一個被稱為 Fisher 的 Z 變換的技巧：
1. 像平常一樣計算 Pearson』s r 的估計值。
2. 用 Fisher 的 Z 變換將 r→z，用公式 z = arctanh(r) 完成。
3. 現在計算 z 的標準差。幸運的是，這很容易計算，由 SDz = 1/sqrt(n-3) 給出，其中 n 是樣本大小。
4. 選擇顯著性閾值，alpha，並檢查與此對應的平均值有多少標準差。如果取 alpha = 0.95，用 1.96。
5. 透過計算 *z* +(1.96 × SD*z*) 找到上限，透過計算 *z -* **(1.96 × SD*z*) 找到下限。
6. 用 r = tanh(z) 將這些轉換回 r。
7. 如果上限和下限都在零的同一側，則有統計顯著性！
這裡是在 Python 中的實現：
r = Pearsons(x,y)
z = math.atanh(r)
SD_z = 1 / math.sqrt(len(x) - 3)
z_upper = z + 1.96 * SD_z
z_lower = z - 1.96 * SD_z
r_upper = math.tanh(z_upper)
r_lower = math.tanh(z_lower)
當然，當給定一個包含許多潛在相關變數的大資料集時，檢查每對的相關性可能很吸引人。這通常被稱為「資料疏浚」——在資料集中查詢變數之間的任何明顯關係。
如果確實採用這種多重比較方法，則應該用適當的更嚴格的顯著性閾值來降低發現錯誤相關性的風險（即找到純粹偶然相關的無關變數）。
一種方法是使用 Bonferroni correction。
小結
到現在為止還好。我們已經看到 Pearson』s *r* 如何用來計算兩個變數之間的相關係數，以及如何評估結果的統計顯著性。給定一組未知的資料，用於開始挖掘變數之間的重要關係是很有可能的。
但是，有一個重要的陷阱——Pearson』s r 只適用於線性資料。
看下面的圖。它們清楚地展示了一種看似非隨機的關係，但是 Pearson』s *r *非常接近於零。
原因是因為這些圖中的變數具有非線性關係。
我們通常可以將兩個變數之間的關係描繪成一個點雲，分散在一條線的兩側。點雲的分散度越大，資料越「嘈雜」，關係越弱。
然而，由於它將每個單獨的資料點與整體平均值進行比較，所以 Pearson』s r 只考慮直線。這意味著檢測非線性關係並不是很好。
在上面的圖中，Pearson』s r 並沒有顯示研究物件的相關性。
然而，這些變數之間的關係很顯然是非隨機的。幸運的是，我們有不同的相關性方法。
讓我們來看看其中幾個。
距離相關性
距離相關性與 Pearson』s *r *有一些相似之處，但是實際上是用一個相當不同的協方差概念來計算的。該方法透過用「距離」類似物替代常用的協方差和標準差（如上所定義）的概念。
類似 Pearson』s r，「距離相關性」被定義為「距離協方差」，由「距離標準差」來歸一化。
距離相關性不是根據它們與各自平均值的距離來估計兩個變數如何共同變化，而是根據與其他點的距離來估計它們是如何共同變化的，從而能更好捕捉變數之間非線性依賴關係。
深入細節
出生於 1773 年的 Robert Brown 是一名蘇格蘭植物學家。當布朗在顯微鏡下研究植物花粉時，注意到液麵上有隨機運動的有機顆粒。
他沒有想到，這一觀察竟使他名垂千古——他成為了布朗運動的（重新）發現者。
他更不會知道，近一個世紀的時間後愛因斯坦才對這種現象做出瞭解釋，從而證實了原子的存在。同年，愛因斯坦發表了關於狹義相對論的論文（E=MC²），並打開了量子理論的大門。
布朗運動是這樣一個物理過程：由於與周圍粒子的碰撞，微小粒子隨機運動。
布朗運動背後的數學原理可以被推廣為維納過程（Weiner process），維納過程在數學金融中最著名的模型 Black-Scholes 中也扮演著重要的角色。
有趣的是，Gabor Szekely 在 20 世紀中期的研究表明，布朗運動和維納過程和一個非線性關聯度量相關。
讓我們來看看如何由長度為 N 的向量 x 和 y 計算這個量。
1. 首先，我們對每個向量構建 N×N 的距離矩陣。距離矩陣和地圖中的道路距離表非常類似——每行、每列的交點顯示了相應城市間的距離。在距離矩陣中，行 i 和列 j 的交點給出了向量的第 i 個元素和第 j 個元素之間的距離。
2. 第二，矩陣是「雙中心」的。也就是說，對於每個元素，我們減去了它的行平均值和列平均值。然後，我們再加上整個矩陣的總平均值。
上述公式中，加「^」表示「雙中心」，加「-」表示「平均值」。
3. 在兩個雙中心矩陣的基礎上，將 X 中每個元素的均值乘以 Y 中相應元素的均值，則可計算出距離協方差的平方。
4. 現在，我們可以用類似的辦法找到「距離方差」。請記住，若兩個向量相同，其協方差與其方差相等。因此，距離方差可表示如下：
5. 最後，我們利用上述公式計算距離相關性。請記住，（距離）標準差與（距離）方差的平方根相等。
如果你更喜歡程式碼實現而非數學符號，那麼請看下面的 R 語言實現：
任意兩變數的距離相關性都在 0 和 1 之間。其中，0 代表兩變數相互獨立，而接近於 1 則表明變數間存在依賴關係。
如果你不想從頭開始編寫距離相關方法，你可以安裝 R 語言的 energy 包（https://cran.r-project.org/web/packages/energy/index.html），設計此方案的研究者提供了本程式碼。在該程式包中，各類可用方案呼叫的是 C 語言編寫的函式，因此有著很大的速度優勢。
置信區間？
我們可以採取「重取樣（resampling）」方法為距離相關性估計建立置信區間。一個簡單的例子是 bootstrap 重取樣。
這是一個巧妙的統計技巧，需要我們從原始資料集中隨機抽樣（替換）以「重建」資料。這個過程將重複多次（例如 1000 次），每次都計算感興趣的統計量。
這將為我們感興趣的統計量產生一系列不同的估計值。我們可以透過它們估計在給定置信水平下的上限和下限。
請看下面的 R 語言程式碼，它實現了簡單的 bootstrap 函式：
如果你想建立統計顯著性，還有另一個重取樣技巧，名為「排列檢驗（permutation test）」。
排列檢驗與上述 bootstrap 方法略有不同。在排列檢驗中，我們保持一個向量不變，並透過重取樣對另一個變數進行「洗牌」。這接近於零假設（null hypothesis）——即，在變數之間不存在依賴關係。
這個經「洗牌」打亂的變數將被用於計算它和常變數間的距離相關性。這個過程將被執行多次，然後，結果的分佈將與實際距離相關性（從未被「洗牌」的資料中獲得）相比較。
然後，大於或等於「實際」結果的經「洗牌」的結果的比例將被定為 P 值，並與給定的顯著性閾值（如 0.05）進行比較。
以下是上述過程的程式碼實現：
最大資訊係數
最大資訊係數（MIC）於 2011 年提出，它是用於檢測變數之間非線性相關性的最新方法。用於進行 MIC 計算的演算法將資訊理論和機率的概念應用於連續型資料。
深入細節
由克勞德·夏農於 20 世紀中葉開創的資訊理論是數學中一個引人注目的領域。
資訊理論中的一個關鍵概念是熵——這是一個衡量給定機率分佈的不確定性的度量。機率分佈描述了與特定事件相關的一系列給定結果的機率。
機率分佈的熵是「每個可能結果的機率乘以其對數後的和」的負值。
為了理解其工作原理，讓我們比較下面兩個機率分佈：
X 軸標明瞭可能的結果；Y 軸標明瞭它們各自的機率
左側是一個常規六面骰子結果的機率分佈；而右邊的六面骰子不那麼均勻。
從直覺上來說，你認為哪個的熵更高呢？哪個骰子結果的不確定性更大？讓我們來計算它們的熵，看看答案是什麼。
不出所料，常規骰子的熵更高。這是因為每種結果的可能性都一樣，所以我們不會提前知道結果偏向哪個。但是，非常規的骰子有所不同——某些結果的發生機率遠大於其它結果——所以它的結果的不確定性也低一些。
這麼一來，我們就能明白，當每種結果的發生機率相同時，它的熵最高。而這種機率分佈也就是傳說中的「均勻」分佈。
交叉熵是熵的一個拓展概念，它引入了第二個變數的機率分佈。
crossEntropy <- function(x,y){
prX <- prop.table(table(x))
prY <- prop.table(table(y))
H <- sum(prX * log(prY,2))
return(-H)
}
兩個相同機率分佈之間的交叉熵等於其各自單獨的熵。但是對於兩個不同的機率分佈，它們的交叉熵可能跟各自單獨的熵有所不同。
這種差異，或者叫「散度」可以透過 KL 散度（Kullback-Leibler divergence）量化得出。
兩機率分佈 X 與 Y 的 KL 散度如下：
機率分佈 X 與 Y 的 KL 散度等於它們的交叉熵減去 X 的熵。
KL 散度的最小值為 0，僅當兩個分佈相同。
KL_divergence <- function(x,y){
kl <- crossEntropy(x,y) - entropy(x)
return(kl)
}
為了發現變數具有相關性，KL 散度的用途之一是計算兩個變數的互資訊（MI）。
互資訊可以定義為「兩個隨機變數的聯合分佈和邊緣分佈之間的 KL 散度」。如果二者相同，MI 值取 0。如若不同，MI 值就為一個正數。二者之間的差異越大，MI 值就越大。
為了加深理解，我們首先簡單回顧一些機率論的知識。
變數 X 和 Y 的聯合機率就是二者同時發生的機率。例如，如果你拋擲兩枚硬幣 X 和 Y，它們的聯合分佈將反映拋擲結果的機率。假設你拋擲硬幣 100 次，得到「正面、正面」的結果 40 次。聯合分佈將反映如下。
P(X=H, Y=H) = 40/100 = 0.4
jointDist <- function(x,y){
N <- length(x)
u <- unique(append(x,y))
joint <- c()
for(i in u){
for(j in u){
f <- x[paste0(x,y) == paste0(i,j)]
joint <- append(joint, length(f)/N)
}
}
return(joint)
}
邊緣分佈是指不考慮其它變數而只關注某一特定變數的機率分佈。假設兩變數獨立，二者邊緣機率的乘積即為二者同時發生的機率。仍以拋硬幣為例，假如拋擲結果是 50 次正面和 50 次反面，它們的邊緣分佈如下：
P(X=H) = 50/100 = 0.5 ; P(Y=H) = 50/100 = 0.5
P(X=H) × P(Y=H) = 0.5 × 0.5 = 0.25
marginalProduct <- function(x,y){
N <- length(x)
u <- unique(append(x,y))
marginal <- c()
for(i in u){
for(j in u){
fX <- length(x[x == i]) / N
fY <- length(y[y == j]) / N
marginal <- append(marginal, fX * fY)
}
}
return(marginal)
}
現在讓我們回到拋硬幣的例子。如果兩枚硬幣相互獨立，邊緣分佈的乘積表示每個結果可能發生的機率，而聯合分佈則為實際得到的結果的機率。
如果兩硬幣完全獨立，它們的聯合機率在數值上（約）等於邊緣分佈的乘積。若只是部分獨立，此處就存在散度。
這個例子中，P(X=H,Y=H) > P(X=H) × P(Y=H)。這表明兩硬幣全為正面的機率要大於它們的邊緣分佈之積。
聯合分佈和邊緣分佈乘積之間的散度越大，兩個變數之間相關的可能性就越大。兩個變數的互資訊定義了散度的度量方式。
X 和 Y 的互資訊等於「二者邊緣分佈積和的聯合分佈的 KL 散度」
mutualInfo <- function(x,y){
joint <- jointDist(x,y)
marginal <- marginalProduct(x,y)
Hjm <- - sum(joint[marginal > 0] * log(marginal[marginal > 0],2))
Hj <- - sum(joint[joint > 0] * log(joint[joint > 0],2))
return(Hjm - Hj)
}
此處的一個重要假設就是機率分佈是離散的。那麼我們如何把這些概念應用到連續的機率分佈呢？
分箱演算法
其中一種方法是量化資料（使變數離散化）。這是透過分箱演算法（bining）實現的，它能將連續的資料點分配對應的離散類別。
此方法的關鍵問題是到底要使用多少「箱子（bin）」。幸運的是，首次提出 MIC 的論文給出了建議：窮舉！
也就是說，去嘗試不同的「箱子」個數並觀測哪個會在變數間取到最大的互資訊值。不過，這提出了兩個挑戰：
1. 要試多少個箱子呢？理論上你可以將變數量化到任意間距值，可以使箱子尺寸越來越小。
2. 互資訊對所用的箱子數很敏感。你如何公平比較不同箱子數目之間的 MI 值？
第一個挑戰從理論上講是不能做到的。但是，論文作者提供了一個啟發式解法（也就是說，解法不完美，但是十分接近完美解法）。他們也給出了可試箱子個數的上限。
最大可用箱子個數由樣本數 N 決定
至於如何公平比較取不同箱子數對 MI 值的影響，有一個簡單的做法……就是歸一化！這可以透過將每個 MI 值除以在特定箱子數組合上取得的理論最大值來完成。我們要採用的是產生最大歸一化 MI 總值的箱子數組合。
互資訊可以透過除以最小的箱子數的對數來歸一化
最大的歸一化互資訊就是 X 和 Y 的最大資訊係數（MIC）。我們來看看一些估算兩個連續變數的 MIC 的程式碼。
MIC <- function(x,y){
N <- length(x)
maxBins <- ceiling(N ** 0.6)
MI <- c()
for(i in 2:maxBins) {
for (j in 2:maxBins){
if(i * j > maxBins){
next
}
Xbins <- i; Ybins <- j
binnedX <-cut(x, breaks=Xbins, labels = 1:Xbins)
binnedY <-cut(y, breaks=Ybins, labels = 1:Ybins)
MI_estimate <- mutualInfo(binnedX,binnedY)
MI_normalized <- MI_estimate / log(min(Xbins,Ybins),2)
MI <- append(MI, MI_normalized)
}
}
return(max(MI))
}
x <- runif(100,-10,10)
y <- x**2 + rnorm(100,0,10)
MIC(x,y) # --> 0.751
以上程式碼是對原論文中方法的簡化。更接近原作的演算法實現可以參考 R package *minerva*（https://cran.r-project.org/web/packages/minerva/index.html）。
在 Python 中的實現請參考 minepy module（https://minepy.readthedocs.io/en/latest/）。
MIC 能夠表示各種線性和非線性的關係，並已得到廣泛應用。它的值域在 0 和 1 之間，值越高表示相關性越強。
置信區間？
為了建立 MIC 估計值的置信區間，你可以簡單地使用一個像我們之前介紹過的 bootstrap 函式。我們可以利用 R 語言的函數語言程式設計，透過傳遞我們想要用作引數的函式來泛化 bootstrap 函式。
bootstrap <- function(x,y,func,reps,alpha){
estimates <- c()
original <- data.frame(x,y)
N <- dim(original)[1]
for(i in 1:reps){
S <- original[sample(1:N, N, replace = TRUE),]
estimates <- append(estimates, func(S$x, S$y))
}
l <- alpha/2 ; u <- 1 - l
interval <- quantile(estimates, c(u, l))
return(2*(func(x,y)) - as.numeric(interval[1:2]))
}
bootstrap(x,y,MIC,100,0.05) # --> 0.594 to 0.88
小結
為了總結相關性這一主題，我們來測試下各演算法在人工生成資料上的處理能力。
完整程式碼：https://gist.github.com/anonymous/fabecccf33f9c3feb568384f626a2c07
噪聲函式
set.seed(123)
# Noise
x0 <- rnorm(100,0,1)
y0 <- rnorm(100,0,1)
plot(y0~x0, pch = 18)
cor(x0,y0)
distanceCorrelation(x0,y0)
MIC(x0,y0)
皮爾遜相關係數 r = - 0.05
距離相關性 = 0.157
MIC = 0.097
簡單線性函式
# Simple linear relationship
x1 <- -20:20
y1 <- x1 + rnorm(41,0,4)
plot(y1~x1, pch =18)
cor(x1,y1)
distanceCorrelation(x1,y1)
MIC(x1,y1)
皮爾遜相關係數 r =+0.95
距離相關性 = 0.95
MIC = 0.89
簡單二次函式
# y ~ x**2
x2 <- -20:20
y2 <- x2**2 + rnorm(41,0,40)
plot(y2~x2, pch = 18)
cor(x2,y2)
distanceCorrelation(x2,y2)
MIC(x2,y2)
Pearson』s r =+0.003
距離相關性 = 0.474
MIC = 0.594
三次函式
# Cosine
x3 <- -20:20
y3 <- cos(x3/4) + rnorm(41,0,0.2)
plot(y3~x3, type="p", pch=18)
cor(x3,y3)
distanceCorrelation(x3,y3)
MIC(x3,y3)
皮爾遜相關係數 r =- 0.035
距離相關性 = 0.382
MIC = 0.484
圓函式
# Circle
n <- 50
theta <- runif (n, 0, 2*pi)
x4 <- append(cos(theta), cos(theta))
y4 <- append(sin(theta), -sin(theta))
plot(x4,y4, pch=18)
cor(x4,y4)
distanceCorrelation(x4,y4)
MIC(x4,y4)
皮爾遜相關係數 r < 0.001
距離相關性 = 0.234
MIC = 0.218