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

Question

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

什么是最快,最pythonic的方式来做到这一点？(循环过来N,M在我看来既不快也不是pythonic.)我期待答案涉及numpy和/或scipy.现在我的数组是numpy arrays,但我愿意将它们转换为不同的类型.

我期待我的输出是一个形状的数组N X M.

NB当我说"相关系数"时,我的意思是Pearson乘积矩相关系数.

以下是一些需要注意的事项:

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

Answer 1

两个2D数组之间的相关性(默认"有效"情况):

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

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

与"valid"两个输入数组的每个成对行组合(row1,row2)之间的默认情况的相关性将对应于每个(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循环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)