首页 文章

Pearson使用scipy.stats.pearsonr与python 3.5中的numpy.corrcoeff对比来自两个2D数组的所有行对之间的相关系数

提问于
浏览 次
0

我试图计算两个2D阵列中每对行之间的Pearson相关系数 . 然后，基于对角线元素对相关矩阵的行/列进行排序 . 首先，在以下代码中从一个随机矩阵（即'randmtx'）计算相关系数矩阵（即'ccmtx'）：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import pearsonr

def correlation_map(x, y):
    n_row_x = x.shape[0]
    n_row_y = x.shape[0]
    ccmtx_xy = np.empty((n_row_x, n_row_y))
    for n in range(n_row_x):
        for m in range(n_row_y):
            ccmtx_xy[n, m] = pearsonr(x[n, :], y[m, :])[0]

    return ccmtx_xy

randmtx = np.random.randn(100, 1000) # generating random matrix
#ccmtx = np.corrcoef(randmtx, randmtx) # cc matrix based on numpy.corrcoef
ccmtx = correlation_map(randmtx, randmtx) # cc matrix based on scipy pearsonr
#
ccmtx_diag = np.diagonal(ccmtx)
#
ids, vals = np.argsort(ccmtx_diag, kind = 'mergesort'), np.sort(ccmtx_diag, kind = 'mergesort')
#ids, vals = np.argsort(ccmtx_diag, kind = 'quicksort'), np.sort(ccmtx_diag, kind = 'quicksort')

plt.plot(ids)
plt.show()

plt.plot(ccmtx_diag[ids])
plt.show()

vals[0]

这里的问题是当使用'pearsonr'时，'ccmtx'的对角线元素恰好是1.0，这是有道理的 . 然而，使用'corrcoef'，'ccmtrix'的对角线元素不是一个（对于某些对角线略小于1），似乎是由于浮点数的精确误差 .

我发现令人讨厌的是单个矩阵的自相关矩阵具有不是1.0的诊断元素，因为当基于对角元素对矩阵进行排序时，这导致相关矩阵的行/列的混洗 .

我的问题是：

[1]当我坚持使用'pearsonr'函数时，有什么好方法可以加快计算时间吗？ （例如，矢量化的皮尔逊？）

[2]在numpy中使用'corrcoef'时，是否有任何好的方法/做法可以防止这种精确错误？ （例如np.around中的'decimals'选项？）

我搜索了两个矩阵的所有行或列对之间的相关系数计算 . 然而，由于算法包含某种“cov / variance”操作，因此这种精确问题似乎总是存在 .

小点：'mergesort'选项似乎提供了比'quicksort'更可靠的结果，因为quicksort改组了1d数组，正好是1到随机顺序 .

任何想法/意见将不胜感激！

python arrays python-3.x numpy scipy

1 回答

  • 0

    对于问题1，矢量化的皮尔逊看到问题的评论 .

    我将只回答问题2：如何提高np.corrcoef的精度 .

    根据协方差矩阵 C 计算相关矩阵 R

    enter image description here
    .

    该实现针对性能和内存使用进行了优化 . 它计算协方差矩阵，然后通过 sqrt(C_ii) 和 sqrt(Cjj) 执行两个除法 . 这个单独的平方根是不精确的来源 . 例如：

    np.sqrt(3 * 3) - 3 == 0.0

np.sqrt(3) * np.sqrt(3) - 3  == -4.4408920985006262e-16

    我们可以通过实现我们自己的简单 corrcoef 例程来解决这个问题：

    def corrcoef(a, b):
    c = np.cov(a, b)
    d = np.diag(c)
    return c / np.sqrt(d[:, None] * d[None, :])

    请注意，此实现需要比numpy实现更多的内存，因为它需要存储大小为n * n的临时矩阵，并且它稍微慢一点，因为它需要执行n ^ 2平方根而不是仅2 n .

    回复于

相关问题