python - 有没有更快的方法来获得相同的结果？

Question

我有两个给定的数组: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 是否有更好更快的方法。

Answer 1

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

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 代码添加到我的答案中)。

---

此版本将采用任何维度 a、b(只要它们相同)，并沿着 pearsonr 计算指定轴。无需 reshape 。

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)