Оптимальный способ вычисления парной взаимной информации с использованием numpy

Для матрицы 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])

Конечно, должны быть лучшие/быстрые способы сделать это?

В стороне я также искал функции сопоставления по столбцам (по столбцам или по строкам) на массивах, но пока не нашел хорошего общего ответа.

Вот моя полная реализация, следуя соглашениям в странице Wiki:

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, который больше похож на декоратора)?
Существуют ли другие оптимальные реализации для этого конкретного вычисления (взаимная информация)?

Ответ 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