计算两个multidimensional array之间的相关系数
我有两个NXT
和MXT
形状的arrays。 我想计算在每一对可能的行n
和m
(分别从N
和M
)之间的T
间的相关系数。
什么是最快,最pythonic方式做到这一点? (在N
和M
循环似乎对我来说既不快也不pythonic)。我期待的答案涉及numpy
和/或scipy
。 现在我的数组是numpy
array
,但我打算将它们转换为不同的types。
我期待我的输出是NXM
形状的arrays。
NB当我说“相关系数”时,我的意思是Pearson乘积矩相关系数 。
这里有一些事情要注意:
-
numpy
函数correlate
要求input数组是一维的。 -
numpy
函数corrcoef
接受二维数组,但它们必须具有相同的形状。 -
scipy.stats
函数pearsonr
要求input数组是一维的。
两个二维数组之间的相关性(默认“有效”情况):
你可以简单地使用像这样的matrix乘法np.dot
–
out = np.dot(arr_one,arr_two.T)
与两个input数组的每个成对行组合(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]))
这是基于这个解决schemeHow to apply corr2 functions in Multidimentional arrays in MATLAB
标杆
本节将运行时性能与针对其他答案中列出的基于generate_correlation_map
&loopy pearsonr
的方法的build议方法进行比较。 (从functiontest_generate_correlation_map()
而在它的末尾没有值正确性validation码)。 请注意,所提出的方法的时间安排还包括在开始时检查两个input数组中相同的列数,如在其他答案中所做的那样。 运行时间在下面列出。
情况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
@Divakar为计算未缩放的关联提供了一个很好的select,这正是我最初要求的。
为了计算相关系数,需要多一点:
import numpy as np def generate_correlation_map(x, y): """Correlate each n with each m. Parameters ---------- x : np.array Shape NX T. y : np.array Shape MX T. Returns ------- np.array NXM 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, yT) - n * np.dot(mu_x[:, np.newaxis], mu_y[np.newaxis, :]) return cov / np.dot(s_x[:, np.newaxis], s_y[np.newaxis, :])
这是对这个函数的testing,它通过:
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)