计算两个多维数组之间的相关系数

Question

我有两个数组，形状分别为 N X T 和 M X T。我想计算每对可能的行 n 和 m 之间 T 的相关系数（分别来自 N 和 M）。

最快、最符合 Python 风格的方法是什么？（循环遍历 N 和 M 在我看来既不快也不符合 pythonic。）我期待答案涉及 numpy and/or scipy。现在我的数组是 numpy arrays，但我愿意将它们转换为不同的类型。

我希望我的输出是一个数组，形状为 N X M。

N.B。当我说 "correlation coefficient," 时，我的意思是 Pearson product-moment correlation coefficient.

这里有一些注意事项：

numpy 函数 correlate 要求输入数组是一维的。
numpy函数corrcoef接受二维数组，但它们必须具有相同的形状。
scipy.stats 函数 pearsonr 要求输入数组是一维的。

Answer 1

两个二维数组之间的相关性（默认 'valid' 情况）：

你可以像这样简单地使用矩阵乘法 np.dot -

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

与两个输入数组的每个成对行组合 (row1,row2) 之间的默认 "valid" 情况的相关性将对应于每个 (row1,row2) 位置的乘法结果。

两个二维数组的行相关系数计算：

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 和循环 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

另一种循环 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

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

Answer 3

对于那些对计算一维和二维数组之间的皮尔逊相关系数感兴趣的人，我编写了以下函数，其中 x 是一维数组，y 是二维数组。

def pearsonr_2D(x, y):
    """computes pearson correlation coefficient
       where x is a 1D and y a 2D array"""

    upper = np.sum((x - np.mean(x)) * (y - np.mean(y, axis=1)[:,None]), axis=1)
    lower = np.sqrt(np.sum(np.power(x - np.mean(x), 2)) * np.sum(np.power(y - np.mean(y, axis=1)[:,None], 2), axis=1))
    
    rho = upper / lower
    
    return rho

示例运行：

>>> x
Out[1]: array([1, 2, 3])

>>> y
Out[2]: array([[ 1,  2,  3],
               [ 6,  7, 12],
               [ 9,  3,  1]])

>>> pearsonr_2D(x, y)
Out[3]: array([ 1.        ,  0.93325653, -0.96076892])

计算两个多维数组之间的相关系数

Computing the correlation coefficient between two multi-dimensional arrays

python

arrays

numpy

scipy

correlation