最大信息系数java_相关性分析之最大互信息系数

最大互信息系数简介

互信息Mutual Information(MI)是用来评价一个事件的出现对于另一个事件的出现所贡献的信息量。在先前聚类算法的评估指标中有过简单的介绍。抛开公式，通俗的理解：原来我对X有些不确定(不确定性为H(X))，告诉我Y后我对X不确定性变为H(X|Y)，这个不确定性的减少量就是X,Y之间的互信息I(X;Y)=H(X)-H(X|Y)。互信息指的是两个随机变量之间的关联程度：即给定一个随机变量后，另一个随机变量不确定性的削弱程度，因而互信息取值最小为0，意味着给定一个随机变量对确定另一个随机变量没有关系，最大取值为随机变量的熵，意味着给定一个随机变量，能完全消除另一个随机变量的不确定性。

把互信息直接用于特征选择存在一些问题：

它不属于度量方式，也没有办法归一化，在不同数据及上的结果无法做比较

对于连续变量的计算不是很方便(X和Y都是集合, $x_i$, y 都是离散的取值)，通常变量需要先离散化，而互信息的结果对离散化的方式很敏感。

MIC 即：Maximal Information Coefficient 最大互信息系数克服了这两个问题。使用MIC来衡量两个基因之间的关联程度，线性或非线性关系，相较于Mutual Information(MI)互信息而言有更高的准确度。它首先寻找一种最优的离散化方式，然后把互信息取值转换成一种度量方式，取值区间在 [0,1]。

MIC的优越性

根据 MIC 的性质，MIC 具有普适性、公平性和对称性。所谓普适性，是指在样本量足够大(包含了样本的大部分信息)时，能够捕获各种各样的有趣的关联，而不限定于特定的函数类型(如线性函数、指数函数或周期函数)，或者说能均衡覆盖所有的函数关系。一般变量之间的复杂关系不仅仅是通过单独一个函数就能够建模的，而是需要叠加函数来表现。所谓公平性，是指在样本量足够大时能为不同类型单噪声程度相似的相关关系给出相近的系数。例如，对于一个充满相同噪声的线性关系和一个正弦关系，一个好的评价算法应该给出相同或相近的相关系数。

算法对比:

最大互信息系数原理

MIC基本原理会利用到互信息概念，互信息的概念使用以下方程来说明：

$$I(x;y)=\int p(x,y) log_2 \frac{p(x,y)}{p(x)p(y)}\mathrm{d}x\mathrm{d}y$$

一般情况下联合概率计算相对来说比较麻烦，MIC的想法是针对两个变量之间的关系离散在二维空间中，并且使用散点图来表示，将当前二维空间在 x,y 方向分别划分为一定的区间数，然后查看当前的散点在各个方格中落入的情况，这样就解决了在互信息中的联合概率难求的问题。下面的公式给出MIC的计算公式：

$$mic(x;y)=\max_{a*b

上式中 a,b 是在 x,y 方向上的划分格子的个数，本质上就是网格分布，B 是变量，在原作者的论文当中提到 B 的大小设置是数据量的 0.6 次方左右。

MIC计算分为三个步骤：

给定i、j，对XY构成的散点图进行i列j行网格化，并求出最大的互信息值

对最大的互信息值进行归一化

选择不同尺度下互信息的最大值作为MIC值

互信息的计算划分方式示例：

那么，给定了某个网格化方案后，如何计算其对应的互信息值呢？这里以上图中红色的网格化方案为例进行说明。红色网格化方案将所有数据点分为四个区域：左上，右上，左下，右下。每个区域对应的数据点数量为1，4，4，1。将数据点数归一化得到四个区域的数据点频率，分别为0.1，0.4，0.4，0.1。也就是说，此时，X有两种取值：左和右，Y有两种取值：上和下。P(X=左,Y=上)=0.1，P(X=右,Y=上)=0.4，P(X=左,Y=下)=0.4，P(X=右,Y=下)=0.1。并且，P(X=左)=0.5，P(X=右)=0.5，P(Y=上)=0.5，P(Y=下)=0.5。根据互信息计算公式，得到X和Y在这种分区下的互信息为：

$$\begin{aligned}&p(X=\text{left}, Y=u p) \log (\frac{p(X=\text{left}, Y=u p)}{p(X=\text{left}) p(Y=u p)})\\ &+p(X=r i g h t, Y=u p) \log (\frac{p(X=\text {right, } Y=u p)}{p(X=\text {right}) p(Y=u p)})\\ &+p(X=\text{left}, Y=\text{down}) \log (\frac{p(X=\text{lef} i, Y=\text {down})}{p(X=\text {left}) p(Y=\text {down})})\\&+p(X=\text {right, } Y=\text {down}) \log (\frac{p(X=\text {right, } Y=\text {down})}{p(X=\text {right}) p(Y=\text {down})})\\ &=0.1 * \log (\frac{0.1}{0.5 * 0.5})+0.4 * \log (\frac{0.4}{0.5 * 0.5})+0.4 * \log (\frac{0.4}{0.5 * 0.5})+0.1 * \log (\frac{0.1}{0.5 * 0.5})\end{aligned}$$

以此类推，算出哪种方案得到的互信息值最大，最大的互信息值是多少。将得到的最大互信息除以log(min(X,Y))，即为归一化。

上面讲述了给定i和j的情况下M(X,Y,D,i,j)的计算方法。这一步就是给定很多(i,j)值，计算每一种情况下M(X,Y,D,i,j)的值，将所有M(X,Y,D,i,j)中的最大那个值作为MIC值。注意的是，这里的(i,j)是有条件的，要满足，n表示数据集D的数据量。当然，B(n)这个值可以自己定，这里是别人做实验认为效果最好的值。

最大互信息系数Python下的使用

在Python中的minepy类库中实现了MIC算法，具体使用如下：

import numpy as np

from minepy import MINE

x = np.linspace(0, 1, 1000)

y = np.sin(10 * np.pi * x) + x

mine = MINE(alpha=0.6, c=15)

mine.compute_score(x, y)

print("Without noise:")

print("MIC", mine.mic())

print()

np.random.seed(0)

y += np.random.uniform(-1, 1, x.shape[0]) # add some noise

mine.compute_score(x, y)

print("With noise:")

print("MIC", mine.mic())

上述代码中的参数解释：

alpha：(float数据类型，取值范围为(0，1] 或 >= 4)，如果alpha的取值范围在(0,1]之内，那么B的取值范围为$(n^{\alpha },4)$，其中N是样本的数目。如果alpha的取值范围是是>= 4。 alpha直接定义B参数。如果alpha高于样本数n，则它将被限制为n，因此B的取值实际上是个分段函数，具体公式为：B = min(alpha，n)。

c(float,取值范围为大于0))，确定比每个分区中的列多个块。默认值为15，这意味着当尝试在x轴上绘制x网格线时，算法将以最多$15 \times x$个团块开始。

结合sklearn的使用：

from minepy import MINE

from sklearn.feature_selection import SelectKBest

from sklearn.datasets import load_iris

import numpy as np

iris = load_iris()

# 由于MINE的设计不是函数式的，定义mic方法将其为函数式的，返回一个二元组，二元组的第2项设置成固定的P值0.5

def mic(x, y):

m = MINE()

m.compute_score(x, y)

return (m.mic(), 0.5)

X_new = SelectKBest(lambda X, Y: tuple(map(tuple, np.array(list(map(lambda x: mic(x, Y), X.T))).T)),

k=2).fit_transform(iris.data, iris.target)

print(X_new) # 输出具体特征列

最大互信息系数可视化

在具体的使用中，有时候我们还需要进行可视化来进行数据探索等等诸多任务。我们使用UCI的红酒质量数据集。然后利用minepy.MINE计算不同特征之间的MIC，然后利用searbon进行矩阵可视化。

import pandas as pd

import numpy as np

from minepy import MINE

import seaborn as sns

import matplotlib.pyplot as plt

# 从硬盘读取数据进入内存

wine = pd.read_csv("data/winequality-red.csv", sep=&＃39;;&＃39;)

def MIC_matirx(dataframe, mine):

data = np.array(dataframe)

n = len(data[0, :])

result = np.zeros([n, n])

for i in range(n):

for j in range(n):

mine.compute_score(data[:, i], data[:, j])

result[i, j] = mine.mic()

result[j, i] = mine.mic()

RT = pd.DataFrame(result)

return RT

def ShowHeatMap(DataFrame):

colormap = plt.cm.RdBu

ylabels = DataFrame.columns.values.tolist()

f, ax = plt.subplots(figsize=(8, 8))

ax.set_title(&＃39;GRA HeatMap&＃39;)

sns.heatmap(DataFrame.astype(float),

cmap=colormap,

ax=ax,

annot=True,

yticklabels=ylabels,

xticklabels=ylabels)

plt.show()

mine = MINE(alpha=0.6, c=15)

data_wine_mic = MIC_matirx(wine, mine)

print(data_wine_mic)

ShowHeatMap(data_wine_mic)

参考链接：