使用numpy计算成对互信息的最佳方法

Question

对于一个 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（例如，一个由浮点数组成的 250000 X 1000 矩阵），我正在进行嵌套的 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循环版本可以以合理的速度完成，但我想知道是否有更优化的方法来应用calc_MI到A的所有列（计算它们的成对互信息）？

我还想知道：

1. 是否有有效的方法将函数映射到np.arrays的列（或行）上进行操作（也许像np.vectorize这样的装饰器）？

2. 是否有其他针对此特定计算（互信息）的最佳实现方式？

- nahsivar

---

你能否扩展你的示例代码，包括 calc_MI 和 A 的示例输入？让它可以被复制、粘贴和运行。这将极大地帮助任何试图回答你问题的人。 - YXD

如果您的矩阵大小为 (n, m)，则没有简单的方法可以将您想要计算的 n * (n - 1) / 2 个唯一值向量化，尽管使用完整笛卡尔积的 n * n 个值进行矢量化计算通常更快，即使存在重复值。问题在于，这需要一次性创建所有中间计算对象。使用上述方法，您需要找出一种创建 4D 的 'histogramdd' 的方法... 我不认为它适用于您的大型数据集。我建议探索Cython或C扩展。 - Jaime

在 shan_entropy 函数中，看起来 H = -sum(...) 应该改为 H = -np.sum(...)。 - Warren Weckesser

是的，E先生。感谢您提供的SSCCE。我只是进行了一些小改动并在这里更新了：http://pastebin.com/2bJM6uSi - nahsivar

是的@M4rtini，已经更正了。感谢您发现了它。不是完全填充零，而是在列之间没有MI的地方填充零。 - nahsivar

Answer 1

针对外循环遍历n*(n-1)/2个向量的快速计算，我无法提供更快的建议。但是，如果您可以使用scipy版本0.13或scikit-learn，则可以简化calc_MI(x, y, bins)的实现。

在scipy 0.13中，scipy.stats.chi2_contingency添加了lambda_参数。该参数控制函数计算的统计量。如果使用lambda_="log-likelihood"（或lambda_=0），则返回对数似然比。这通常也称为G或G²统计量。除了样本表中的总样本数n的因子2*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”而不是“bits”）。如果您真的喜欢“bits”，只需将mi除以log(2)即可。

如果您有（或可以安装）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

- Warren Weckesser

---

好代码！对于箱子数量，合理的默认值是多少？ - pir

@felbo 这是个好问题，但并不容易回答。如果你在http://stats.stackexchange.com/上提问的话，或许会得到一些灵感。 - Warren Weckesser

两种建议的方法在连续校正上有所不同。将 chi2_contingency(correction=False) 改为此，即可消除这种不一致性。 - shouldsee

小心离散化（即箱子大小）。本博客建议使用“Jack Knifed Estimate”来克服这个问题，或者我也希望密度估计技术能够提供帮助。 - Josh.F

这个实现和 https://scikit-learn.org/stable/modules/generated/sklearn.feature_selection.mutual_info_classif.html?highlight=mutual_info_classif#sklearn.feature_selection.mutual_info_classif 有什么区别？ - dada

Answer 2

为了摆脱外部循环（在某种程度上），一种方法是重写calc_MI以调用在整个数组c_XY上使用的向量化函数，这些函数用于构建matMI。

shan_entropy使用可以处理任意大小数组的函数，c_X和c_Y是np.histogram2d输出的第一维和第二维上的边际总数。同样，对matMI的循环赋值可以使用np.triu_indices进行向量化处理。

这里唯一无法避免的循环是对A的列对进行循环，以调用np.histogram2d，但这部分可以使用joblib.Parallel进行并行化处理。

def shan_entropy(c):
    return - (c * np.log2(c, out=np.zeros(c.shape), where=(c!=0))).sum(axis=1)

def calc_pairwise_mutual_info(A, bins):
    
    m, n = A.shape
    matMI = np.zeros((n, n))

    c_XYs = np.array([np.histogram2d(A[:,ix], A[:,jx], bins)[0] for ix in range(n-1) for jx in range(ix+1, n)]) / m
    c_Xs = c_XYs.sum(axis=2)
    c_Ys = c_XYs.sum(axis=1)
    c_XYs = c_XYs.reshape(len(c_XYs), -1)

    H_X, H_Y, H_XY = map(shan_entropy, (c_Xs, c_Ys, c_XYs))

    MI = H_X + H_Y - H_XY
    matMI[np.triu_indices(n, k=1)] = MI
    return matMI


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

MI_arr = calc_pairwise_mutual_info(A, bins)

您可以验证MI_arr与OP中计算的matMI相同（np.allclose返回True）。

- cottontail

Answer 3

你还可以使用scipy.stat.entropy。

- Mohammad Reza Rezaie

---

你的回答可以通过提供更多支持性信息来改进。请编辑以添加进一步的细节，例如引用或文档，以便他人能够确认你的回答是否正确。你可以在帮助中心找到关于如何撰写好回答的更多信息。 - Community

虽然这个链接可能回答了问题，但最好在这里包含答案的关键部分，并提供链接作为参考。仅有链接的答案如果链接页面发生变化，就可能失效。- 来自评论 - Chenmunka