我有两个给定的数组:x 和 y。我想计算两个数组之间的相关系数如下:

import numpy as np
from scipy.stats import pearsonr

x = np.array([[[1,2,3,4],
               [5,6,7,8]],
              [[11,22,23,24],
               [25,26,27,28]]])


i,j,k = x.shape

y = np.array([[[31,32,33,34],
               [35,36,37,38]],
              [[41,42,43,44],
               [45,46,47,48]]])



xx = np.row_stack(np.dstack(x))
yy = np.row_stack(np.dstack(y))

results = []

for a, b in zip(xx,yy):
    r_sq, p_val = pearsonr(a, b)
    results.append(r_sq)

results = np.array(results).reshape(j,k)

print results

[[ 1.  1.  1.  1.]
 [ 1.  1.  1.  1.]]

答案是正确的。但是,想知道是否有更好更快的方法来使用 numpy 和/或 scipy。

最佳答案

另一种方法(不一定更好)是:

xx = x.reshape(2,-1).T  # faster, minor issue though
yy = y.reshape(2,-1).T
results = [pearsonr(a,b)[0] for a,b in zip(xx,yy)]
results = np.array(results).reshape(x.shape[1:])

当前另一个线程正在讨论使用列表推导式迭代数组的值:Confusion about numpy's apply along axis and list comprehensions

正如那里所讨论的,另一种方法是初始化 results ,并在迭代期间填充值。对于非常大的情况,这可能会更快,但对于适度的情况,这
np.array([... for .. in ...])

是合理的。

更深层次的问题是 pearsonr 或其他替代方法是否可以计算许多对的相关性,而不仅仅是一对。这可能需要研究 pearsonr 的内部结构,或 stats 中的其他函数。

这是向量化 stats.pearsonr 的第一个切入点:
def pearsonr2(a,b):
    # stats.pearsonr adapted for
    # x and y are (N,2) arrays
    n = x.shape[1]
    mx = x.mean(1)
    my = y.mean(1)
    xm, ym = x-mx[:,None], y-my[:,None]
    r_num = np.add.reduce(xm * ym, 1)
    r_den = np.sqrt(stats.ss(xm,1) * stats.ss(ym,1))
    r = r_num / r_den
    r = np.clip(r, -1.0, 1.0)
    return r

print pearsonr2(xx,yy)

它符合您的情况,尽管这些测试值并没有真正行使该功能。我只是采用了 pearsonr 代码,在大多数行中添加了 axis=1 参数,并确保一切正常。 prob 步骤可以包含在一些 bool 掩码中。

(如果需要,我可以将 stats.pearsonr 代码添加到我的答案中)。

此版本将采用任何维度 ab (只要它们相同),并沿指定的轴进行 pearsonr calc。无需重塑。
def pearsonr_flex(a,b, axis=1):
    # stats.pearsonr adapted for
    # x and y are (N,2) arrays
    n = x.shape[axis]
    mx = x.mean(axis, keepdims=True)
    my = y.mean(axis, keepdims=True)
    xm, ym = x-mx, y-my
    r_num = np.add.reduce(xm * ym, axis)
    r_den = np.sqrt(stats.ss(xm, axis) * stats.ss(ym, axis))
    r = r_num / r_den
    r = np.clip(r, -1.0, 1.0)
    return r

pearsonr_flex(xx, yy, 1)
preasonr_flex(x, y, 0)

关于python - 有没有更快的方法来获得相同的结果?,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/27649503/

10-12 17:17