Вычисление коэффициента корреляции между двумя многомерными массивами

У меня есть два массива, которые имеют формы N X T и M X T. Я хотел бы вычислить коэффициент корреляции через T между всеми возможными парами строк n и m (из n и m соответственно).

Какой самый быстрый, самый питонический способ сделать это? (Циклирование над n и m показалось бы мне не быстрым и не pythonic.) Я ожидаю, что ответ будет включать numpy и/или scipy. Прямо сейчас мои массивы numpy array s, но я открыт для их преобразования в другой тип.

Я ожидаю, что мой вывод будет массивом с формой N X M.

N.B. Когда я говорю "коэффициент корреляции", я имею в виду коэффициент корреляции момента продукта Пирсона.

Вот что нужно отметить:

Функция numpy correlate требует, чтобы массивы ввода были одномерными.
Функция numpy corrcoef принимает двумерные массивы, но они должны иметь одинаковую форму.
Функция scipy.stats pearsonr требует, чтобы массивы ввода были одномерными.

Ответ 1

Корреляция (по умолчанию "действительный" случай) между двумя 2D-массивами:

Вы можете просто использовать матричное умножение np.dot, как это -

out = np.dot(arr_one,arr_two.T)

Корреляция со случаем по умолчанию "valid" между каждой комбинацией рядов строк (строка1, строка2) двух входных массивов будет соответствовать результату умножения в каждой позиции (row1, row2).

Расчет коэффициента корреляции по шкале для двух 2D-массивов:

def corr2_coeff(A,B):
    # Rowwise mean of input arrays & subtract from input arrays themeselves
    A_mA = A - A.mean(1)[:,None]
    B_mB = B - B.mean(1)[:,None]

    # Sum of squares across rows
    ssA = (A_mA**2).sum(1);
    ssB = (B_mB**2).sum(1);

    # Finally get corr coeff
    return np.dot(A_mA,B_mB.T)/np.sqrt(np.dot(ssA[:,None],ssB[None]))

Это основано на этом решении How to apply corr2 functions in Multidimentional arrays in MATLAB

Бенчмаркинг

В этом разделе сравнивается производительность выполнения с предлагаемым подходом к подходу generate_correlation_map и loopy pearsonr, указанному в другом ответе. (взято из функции test_generate_correlation_map() без кода проверки правильности значения в конце его). Обратите внимание, что тайминги предлагаемого подхода также включают проверку в начале проверки одинакового количества столбцов в двух входных массивах, также как и в этом другом ответе. Время выполнения указано ниже.

Случай №1:

In [106]: A = np.random.rand(1000,100)

In [107]: B = np.random.rand(1000,100)

In [108]: %timeit corr2_coeff(A,B)
100 loops, best of 3: 15 ms per loop

In [109]: %timeit generate_correlation_map(A, B)
100 loops, best of 3: 19.6 ms per loop

Случай № 2:

In [110]: A = np.random.rand(5000,100)

In [111]: B = np.random.rand(5000,100)

In [112]: %timeit corr2_coeff(A,B)
1 loops, best of 3: 368 ms per loop

In [113]: %timeit generate_correlation_map(A, B)
1 loops, best of 3: 493 ms per loop

Случай №3:

In [114]: A = np.random.rand(10000,10)

In [115]: B = np.random.rand(10000,10)

In [116]: %timeit corr2_coeff(A,B)
1 loops, best of 3: 1.29 s per loop

In [117]: %timeit generate_correlation_map(A, B)
1 loops, best of 3: 1.83 s per loop

Другой метод loopy pearsonr based казался слишком медленным, но вот время выполнения для одного небольшого размера данных -

In [118]: A = np.random.rand(1000,100)

In [119]: B = np.random.rand(1000,100)

In [120]: %timeit corr2_coeff(A,B)
100 loops, best of 3: 15.3 ms per loop

In [121]: %timeit generate_correlation_map(A, B)
100 loops, best of 3: 19.7 ms per loop

In [122]: %timeit pearsonr_based(A,B)
1 loops, best of 3: 33 s per loop

Ответ 2

@Divakar предоставляет отличный вариант для вычисления немасштабированной корреляции, о чем я и просил.

Для вычисления коэффициента корреляции требуется немного больше:

import numpy as np

def generate_correlation_map(x, y):
    """Correlate each n with each m.

    Parameters
    ----------
    x : np.array
      Shape N X T.

    y : np.array
      Shape M X T.

    Returns
    -------
    np.array
      N X M array in which each element is a correlation coefficient.

    """
    mu_x = x.mean(1)
    mu_y = y.mean(1)
    n = x.shape[1]
    if n != y.shape[1]:
        raise ValueError('x and y must ' +
                         'have the same number of timepoints.')
    s_x = x.std(1, ddof=n - 1)
    s_y = y.std(1, ddof=n - 1)
    cov = np.dot(x,
                 y.T) - n * np.dot(mu_x[:, np.newaxis],
                                  mu_y[np.newaxis, :])
    return cov / np.dot(s_x[:, np.newaxis], s_y[np.newaxis, :])

Здесь проверяется эта функция, которая проходит:

from scipy.stats import pearsonr

def test_generate_correlation_map():
    x = np.random.rand(10, 10)
    y = np.random.rand(20, 10)
    desired = np.empty((10, 20))
    for n in range(x.shape[0]):
        for m in range(y.shape[0]):
            desired[n, m] = pearsonr(x[n, :], y[m, :])[0]
    actual = generate_correlation_map(x, y)
    np.testing.assert_array_almost_equal(actual, desired)