MI & MIC
互信息(Mutual Information, MI)是一个统计量,用来衡量两个随机变量之间的依赖程度或信息共享的多少。简单来说,就是它告诉我们知道一个变量后,能减少多少对另一个变量的不确定性。
最大互信息系数(Maximal Information Coefficient, MIC)是一个更先进的统计量,用来检测两个变量之间的关系。不像互信息只关注信息共享,MIC可以捕捉到各种复杂关系,包括线性和非线性的。MIC的主要特点是它能在不同情况下识别出不同的关系类型,并且不需要事先假设关系的形式。简单来说,它能灵活地找到最能描述两个变量之间关系的模式。和R Square有些类似地,MIC的取值为[0, 1],1表示完全相关,0则表示毫无关系。通过这种方式,我们能够更有效地挖掘数据中的隐含关系。
计算互信息
最简单的方式就是调用sklearn的方法。需要注意的是,如果待分析的数组包含连续变量,则需要先离散化,否则计算结果会不准确。
import numpy as np
from sklearn.metrics import mutual_info_score
from sklearn.preprocessing import KBinsDiscretizer
'''
传入以下函数的x和y均为numpy数组,且等长
'''
# 离散化方法1
def calc_mi(x, y, bins=10): # bins是离散化时的分箱数
# 使用直方图来估计联合概率分布
c_xy = np.histogram2d(x, y, bins)[0]
# 计算互信息
mi = mutual_info_score(None, None, contingency=c_xy)
return mi
# 离散化方法2,更精细一些
def calc_mi_2(x, y, bins=10):
'''
KBinsDiscretizer():用于将连续数据转换成离散数据
params:
n_bins: 分箱数
encode: 编码为整数值('ordinal')或者独热编码('onehot')等
strategy: 等距离分箱('uniform')、分位数分箱('quantile')或者聚类分箱('kmeans')等
'''
est = KBinsDiscretizer(n_bins=bins, encode='ordinal', strategy='quantile')
x_binned = est.fit_transform(x.reshape(-1, 1)).astype(int).flatten()
y_binned = est.fit_transform(y.reshape(-1, 1)).astype(int).flatten()
# 计算互信息
mi = mutual_info_score(x_binned, y_binned)
return x_binned, y_binned, mi
互信息的随机化检验
随机化检验(Permutation Test)是一种非参数统计方法,用来评估统计量(比如互信息)的显著性。具体来说,它帮助我们确定观测到的互信息是否显著高于随机情况下的值,从而判断两个变量之间的关系是否非随机。
步骤通常如下:
- 计算原始互信息: 首先计算实际数据 ( x ) 和 ( y ) 之间的互信息值。
- 生成随机化数据: 通过打乱其中一个变量的数据顺序,生成许多随机化的数据集。在每次随机化时,我们假定两个变量之间没有关系。
- 计算随机互信息分布: 对每个随机化数据集计算互信息,得到一个互信息的分布。
- 评估显著性: 将原始互信息值与随机化互信息分布进行比较,计算原始互信息值在随机化分布中的百分位数(p-value)。如果 p-value 很小(例如小于0.05),我们就可以认为原始互信息显著高于随机情况下的值,从而认为两变量之间存在显著关系。
import numpy as np
from sklearn.metrics import mutual_info_score
'''
传入以下函数的位置参数均由上述calc_mi_2()函数生成
'''
def permutation_test(x_binned, y_binned, mi, num_permutations=1000): # num_permutations是重复测试次数
permuted_mi = np.zeros(num_permutations)
for _ in range(num_permutations):
np.random.shuffle(y_binned)
permuted_mi[_] = mutual_info_score(x_binned, y_binned)
p_value = np.mean(permuted_mi >= mi)
return p_value
计算最大互信息系数
全网推荐最多的是使用minepy库的方法。但这个库的安装条件极其苛刻,因此本文弃用该思路。但为了保持方法的通用性,默认参数沿用minepy库的定义法。
其中,参数alpha
控制网格分辨率的范围,具体来说,它决定了用于计算互信息的网格的最大尺寸(对数据点的平方根进行缩放);参数c
控制每个网格中的数据点对的数量上限,用于防止过拟合。
import numpy as np
def mic(x, y, alpha=0.6, c=15):
def mutual_info(x_bins, y_bins): # 和前面不同,这里计算MI是纯numpy的方法
# 计算联合概率分布
joint_prob, _, _ = np.histogram2d(x, y, bins=(x_bins, y_bins))
joint_prob /= joint_prob.sum()
# 计算边缘概率分布
x_prob = joint_prob.sum(axis=1)
y_prob = joint_prob.sum(axis=0)
# 计算互信息
mi = 0
for i in range(len(x_prob)):
for j in range(len(y_prob)):
if joint_prob[i, j] > 0:
mi += joint_prob[i, j] * np.log(joint_prob[i, j] / (x_prob[i] * y_prob[j]))
return mi
def entropy(prob):
return -np.sum([p * np.log(p) for p in prob if p > 0])
n = len(x)
mic = 0
# 网格分辨率范围
for nx in range(2, int(alpha * n**0.5) + 2):
for ny in range(2, int(alpha * n**0.5) + 2):
# 限制数据点对的数量
if nx * ny <= n / c:
mi = mutual_info(nx, ny)
h = min(entropy(np.histogram(x, bins=nx)[0] / n), entropy(np.histogram(y, bins=ny)[0] / n))
mic = max(mic, mi / h if h > 0 else 0)
return mic