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

7gcisfzg  于 2022-11-10  发布在  Python
关注(0)|答案(1)|浏览(158)

我试图改进一个简单的算法,从两个数组X(m,n)Y(n)中获得皮尔逊相关系数,返回另一个维度为(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值(接近10 k),这运行得相当快,但是我开始使用m ~ 30 k,所以这花费的时间比我预期的要长得多。我知道我可以实现多处理/多线程,但是我相信有一种(更好的)pythonic 方法来实现这一点。
我尝试使用pearsonr(x, np.ones((m, n))*y),但它只返回(nan, nan)

8ljdwjyq

8ljdwjyq1#

pearsonr只在内部支持1D数组。此外,它计算p值,而这里没有使用。因此,如果可能的话,不计算它会更有效。此外,代码每次都重新计算y向量,它不能有效地使用向量化的Numpy运算。这就是计算有点慢的原因。您可以在这里的代码中检查这一点。
计算此值的一种方法是基于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)

在我的机器上,m = 10_000时,速度快了450倍
请注意,我没有保留Scipy代码的检查,但如果您的输入不能保证统计安全(即为皮尔逊检验的计算格式良好),保留它们可能是一个好主意。

相关问题