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

Posted 2023-02-23

【中文标题】使用 numpy 计算成对互信息的最佳方法

【英文标题】：Optimal way to compute pairwise mutual information using numpy

【发布时间】：2013-12-27 18:24:11

【问题描述】：

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

mutual information，我的意思是：

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])

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

顺便说一句，我还在数组的列（按列或按行操作）上寻找映射函数，但还没有找到一个好的通用答案。

这是我的完整实现，遵循the Wiki page 中的约定：

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 的所有列（以成对计算它们互信息）？

我也想知道：

是否有有效的方法来映射函数以对np.arrays 的列（或行）进行操作（可能像np.vectorize，看起来更像一个装饰器）？

对于这个特定的计算（互信息）是否还有其他优化的实现方式？

【问题讨论】：

您能否扩展您的示例代码以包含calc_MI 和A 的示例输入？制作它，以便我们可以复制、粘贴和运行。将极大地帮助任何试图回答您的问题的人。

请阅读此sscce.org 并更新您的示例代码以包含calc_MI 和A 的示例输入。

我之前的评论是在我打算回应建议时无意中输入的。感谢您指向 sscce.org。

对于您当前的方法，这是一个准确的自包含示例吗？ pastebin.com/kbzyvA6K.

如果您的矩阵大小为(n, m)，则没有简单的方法可以仅对您所追求的n * (n - 1) / 2 唯一值的计算进行矢量化，尽管进行矢量化计算通常更快完整笛卡尔积中的n * n 值，即使有重复项。这样做的问题是它需要一次创建所有中间计算对象。使用上述方法，您将不得不想出一种创建 4D histogramdd 的方法...我认为它不适用于您的庞大数据集。我会研究 Cython 或 C 扩展...

【参考方案1】：

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

在 scipy 0.13 中，lambda_ 参数被添加到 scipy.stats.chi2_contingency 此参数控制函数计算的统计量。如果 你使用lambda_="log-likelihood"（或lambda_=0），对数似然比 被退回。这通常也称为 G 或 G² 统计量。以外 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”而不是“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

【讨论】：

不错的代码！ bin 数量的合理默认值是多少？

@felbo 这是一个很好的问题，而且不是一个容易回答的问题。如果您在stats.stackexchange.com

询问，您可能会得到一些想法

如果某些计数为零，则此方法不起作用。为什么你建议这种方法超过密度估计？另外，我赞成您的回答，因为它确实提供了一种在某些情况下计算 MI 的有效方法。

“你为什么推荐这种方法而不是密度估计？”我没有。我只建议了问题中给出的代码的几个替代实现。

建议的两种方法因连续性校正而异。更改为 chi2_contingency(correction = False) 消除了这种不一致。