有更快的方法来获得相同的结果吗？

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 987654326 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 STEP可以包含一些布尔屏蔽。

（我可以在需要时将stats.pearsonr代码添加到我的答案中）。

---

此版本将采用任何维度a，b（只要它们是相同的），并沿指定轴进行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)