利用数值计算成对互信息的最优方法

对于 M x n矩阵，计算所有列对(N x n)的互信息的最佳(最快)方法是什么？

我说的 共同信息是指:

I (X，Y) = H (X) + H (Y)-H (X，Y)

其中 H (X)指的是 X的香农熵。

目前我正在使用 np.histogram2d和 np.histogram来计算 (X，Y)和单个 (X 或 Y)的联合计数。对于给定的矩阵 A(例如一个250000X1000的浮点矩阵) ，我正在做一个嵌套的 for循环,

    n = A.shape[1]
for ix = arange(n)
for jx = arange(ix+1,n):
matMI[ix,jx]= calc_MI(A[:,ix],A[:,jx])

肯定有更好/更快的方法来做到这一点吧？

另外，我还在数组的列上寻找映射函数(列操作或行操作) ，但是还没有找到一个很好的通用答案。

以下是我的完整实现，遵循 维基页面的惯例:

import numpy as np


def calc_MI(X,Y,bins):


c_XY = np.histogram2d(X,Y,bins)[0]
c_X = np.histogram(X,bins)[0]
c_Y = np.histogram(Y,bins)[0]


H_X = shan_entropy(c_X)
H_Y = shan_entropy(c_Y)
H_XY = shan_entropy(c_XY)


MI = H_X + H_Y - H_XY
return MI


def shan_entropy(c):
c_normalized = c / float(np.sum(c))
c_normalized = c_normalized[np.nonzero(c_normalized)]
H = -sum(c_normalized* np.log2(c_normalized))
return H


A = np.array([[ 2.0,  140.0,  128.23, -150.5, -5.4  ],
[ 2.4,  153.11, 130.34, -130.1, -9.5  ],
[ 1.2,  156.9,  120.11, -110.45,-1.12 ]])


bins = 5 # ?
n = A.shape[1]
matMI = np.zeros((n, n))


for ix in np.arange(n):
for jx in np.arange(ix+1,n):
matMI[ix,jx] = calc_MI(A[:,ix], A[:,jx], bins)

虽然我的嵌套 for循环的工作版本可以以合理的速度完成，但是我想知道是否有一种更优化的方法可以在 A的所有列上应用 calc_MI(计算它们成对的相互信息) ？

我还想知道:

  1. 是否有有效的方法来映射函数来操作 np.arrays的列(或行)(也许像 np.vectorize，它看起来更像一个装饰器) ？

  2. 是否还有其他针对此特定计算的最佳实现(互信息) ？

57016 次浏览
小开

我不能建议对 n * (n-1)/2上的外循环进行更快的计算 向量，但是可以简化 calc_MI(x, y, bins)的实现 如果您可以使用 scipy 版本0.13或 Scikit-learn。

在 scipy 0.13中，将 lambda_参数添加到 < a href = “ http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chi2 _ contency.html”rel = “ noReferrer”> scipy.stats.chi2_contingency 此参数控制由函数计算的统计信息。如果 你使用 lambda_="log-likelihood"(或 lambda_=0) ，对数似然比 返回。这也通常称为 G 或 G2统计量。除了 A 因子为2 * n (其中 n 是偶然性中的样本总数 表) ，这个 是的相互信息。因此，你可以实现 calc_MI 作为:

from scipy.stats import chi2_contingency


def calc_MI(x, y, bins):
c_xy = np.histogram2d(x, y, bins)[0]
g, p, dof, expected = chi2_contingency(c_xy, lambda_="log-likelihood")
mi = 0.5 * g / c_xy.sum()
return mi

它与您的实现之间的唯一区别是 实现使用自然对数代替以2为基数的对数 (因此它用“ nats”而不是“ bit”表示信息) 你真的喜欢位，只要除以日志(2)的 mi。

如果您有(或可以安装) sklearn(即 scikit-learn) ，则可以使用 sklearn.metrics.mutual_info_score ，并按照以下方式实施 calc_MI:

from sklearn.metrics import mutual_info_score


def calc_MI(x, y, bins):
c_xy = np.histogram2d(x, y, bins)[0]
mi = mutual_info_score(None, None, contingency=c_xy)
return mi

提交答案