如何更好地从 2 个维度 (m, n) 和 (n) 的数组执行 Pearson R，返回一个 (m) 大小的数组？ [Python, NumPy, SciPy]

Question

我正在尝试改进一个简单的算法以从两个数组 X(m, n) 和 Y(n) 中获取 Pearson 相关系数，返回另一个维度为 (m).
的数组 R 在这种情况下，我想知道 X 的每一行关于 Y 的行为。示例（工作）代码如下所示：

import numpy as np
from scipy.stats import pearsonr

np.random.seed(1)
m, n = 10, 5

x = 100*np.random.rand(m, n)
y = 2 + 2*x.mean(0)
r = np.empty(m)

for i in range(m):
    r[i] = pearsonr(x[i], y)[0]

对于这种特殊情况，我得到：r = array([0.95272843, -0.69134753, 0.36419159, 0.27467137, 0.76887201, 0.08823868, -0.72608421, -0.01224453, 0.58375626, 0.87442889])

对于 m 的小值（接近 10k），它运行得非常快，但我开始使用 m ~ 30k ，所以这比我预期的要长得多。我知道我可以实现 multiprocessing/multi-threading 但我相信有一种（更好的）pythonic 方法可以做到这一点。

我尝试使用 pearsonr(x, np.ones((m, n))*y)，但它 returns 仅 (nan, nan)。

Answer 1

pearsonr 内部只支持一维数组。此外，它计算此处未使用的 p-values 。因此，如果可能的话，不计算它会更有效。此外，该代码还每次都重新计算 y 向量，并且它没有有效地利用向量化的 Numpy 操作。这就是计算速度有点慢的原因。您可以在代码 here.

中查看

一种计算方法是根据 Scipy:

编写您自己的自定义实现

def multi_pearsonr(x, y):
    xmean = x.mean(axis=1)
    ymean = y.mean()
    xm = x - xmean[:,None]
    ym = y - ymean
    normxm = np.linalg.norm(xm, axis=1)
    normym = np.linalg.norm(ym)
    return np.clip(np.dot(xm/normxm[:,None], ym/normym), -1.0, 1.0)

在我的机器上 450 倍 m = 10_000。

请注意，我没有保留 Scipy 代码的检查，但如果不能保证您的输入在统计上是安全的（即计算格式正确），保留它们可能是个好主意皮尔逊测试）。

如何更好地从 2 个维度 (m, n) 和 (n) 的数组执行 Pearson R，返回一个 (m) 大小的数组？ [Python, NumPy, SciPy]

How better perform Pearson R from 2 arrays of dimensions (m, n) and (n), returning an array of (m) size? [Python, NumPy, SciPy]

python

numpy

scipy