https://blog.csdn.net/qtlyx/article/details/50780400
最后的效果就是這樣的。很明顯可以看到，左下角那個有點像三角函數的關系，Pearson系數(就是線性相關系數)為0，而MIC則有0.8。
摘自：http://tech.ifeng.com/a/20180323/44917506_0.shtml
最大信息系數
最大信息系數（MIC）於 2011 年提出，它是用於檢測變量之間非線性相關性的最新方法。用於進行 MIC 計算的算法將信息論和概率的概念應用於連續型數據。
深入細節
由克勞德·香農於 20 世紀中葉開創的信息論是數學中一個引人注目的領域。
信息論中的一個關鍵概念是熵——這是一個衡量給定概率分布的不確定性的度量。概率分布描述了與特定事件相關的一系列給定結果的概率。
概率分布的熵是「每個可能結果的概率乘以其對數后的和」的負值
為了理解其工作原理，讓我們比較下面兩個概率分布：
X 軸標明了可能的結果；Y 軸標明了它們各自的概率
左側是一個常規六面骰子結果的概率分布；而右邊的六面骰子不那么均勻。
從直覺上來說，你認為哪個的熵更高呢？哪個骰子結果的不確定性更大？讓我們來計算它們的熵，看看答案是什么。
entropy <- function(x){
pr <- prop.table(table(x))
H <- sum(pr * log(pr,2))
return(-H)
}
dice1 <- 1:6
dice2 <- c(1,1,1,1,2:6)
entropy(dice1) # --> 2.585
entropy(dice2) # --> 2.281不出所料，常規骰子的熵更高。這是因為每種結果的可能性都一樣，所以我們不會提前知道結果偏向哪個。但是，非常規的骰子有所不同——某些結果的發生概率遠大於其它結果——所以它的結果的不確定性也低一些。
這么一來，我們就能明白，當每種結果的發生概率相同時，它的熵最高。而這種概率分布也就是傳說中的「均勻」分布。
交叉熵是熵的一個拓展概念，它引入了第二個變量的概率分布。
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 的論文給出了建議：窮舉！
也就是說，去嘗試不同的「箱子」個數並觀測哪個會在變量間取到最大的互信息值。不過，這提出了兩個挑戰：
-
要試多少個箱子呢？理論上你可以將變量量化到任意間距值，可以使箱子尺寸越來越小。
-
互信息對所用的箱子數很敏感。你如何公平比較不同箱子數目之間的 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 之間，值越高表示相關性越強。