我试图改进一个简单的算法,从两个数组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)
。
1条答案
按热度按时间8ljdwjyq1#
pearsonr
只在内部支持1D数组。此外,它计算p值,而这里没有使用。因此,如果可能的话,不计算它会更有效。此外,代码每次都重新计算y
向量,它不能有效地使用向量化的Numpy运算。这就是计算有点慢的原因。您可以在这里的代码中检查这一点。计算此值的一种方法是基于Scipy之一编写您自己的自定义实现:
在我的机器上,m = 10_000时,速度快了450倍。
请注意,我没有保留Scipy代码的检查,但如果您的输入不能保证统计安全(即为皮尔逊检验的计算格式良好),保留它们可能是一个好主意。