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

Question

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

最快、最 Pythonic 的方法是什么？ (在我看来，循环 N 和 M 既不快也不像pythonic。)我期待答案涉及 numpy 和/或 scipy。现在我的数组是 numpy arrays，但我愿意将它们转换为不同的类型。

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

注意当我说“相关系数”时，我的意思是 Pearson product-moment correlation coefficient .

这里有一些注意事项:

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

Answer 1

两个二维数组之间的相关性(默认为“有效”情况):

您可以简单地使用矩阵乘法 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 的解决方案

基准测试

本节将运行时性能与建议的方法与 other answer. 中列出的基于 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