python 相关性矩阵_矩阵所有行对的相关系数和p值

我今天也遇到同样的问题。

经过半个小时的搜索，我在numpy/scipy库中找不到任何代码可以帮助我做到这一点。

所以我写了自己的版本import numpy as np

from scipy.stats import pearsonr, betai

def corrcoef(matrix):

r = np.corrcoef(matrix)

rf = r[np.triu_indices(r.shape[0], 1)]

df = matrix.shape[1] - 2

ts = rf * rf * (df / (1 - rf * rf))

pf = betai(0.5 * df, 0.5, df / (df + ts))

p = np.zeros(shape=r.shape)

p[np.triu_indices(p.shape[0], 1)] = pf

p[np.tril_indices(p.shape[0], -1)] = p.T[np.tril_indices(p.shape[0], -1)]

p[np.diag_indices(p.shape[0])] = np.ones(p.shape[0])

return r, p

def corrcoef_loop(matrix):

rows, cols = matrix.shape[0], matrix.shape[1]

r = np.ones(shape=(rows, rows))

p = np.ones(shape=(rows, rows))

for i in range(rows):

for j in range(i+1, rows):

r_, p_ = pearsonr(matrix[i], matrix[j])

r[i, j] = r[j, i] = r_

p[i, j] = p[j, i] = p_

return r, p

第一个版本使用np.corrcoef的结果，然后基于corrcoef矩阵的三角形上数值计算p值。

第二个循环版本只是在行上迭代，手动执行pearsonr。def test_corrcoef():

a = np.array([

[1, 2, 3, 4],

[1, 3, 1, 4],

[8, 3, 8, 5],

[2, 3, 2, 1]])

r1, p1 = corrcoef(a)

r2, p2 = corrcoef_loop(a)

assert np.allclose(r1, r2)

assert np.allclose(p1, p2)

考试通过了，他们是一样的。def test_timing():

import time

a = np.random.randn(100, 2500)

def timing(func, *args, **kwargs):

t0 = time.time()

loops = 10

for _ in range(loops):

func(*args, **kwargs)

print('{} takes {} seconds loops={}'.format(

func.__name__, time.time() - t0, loops))

timing(corrcoef, a)

timing(corrcoef_loop, a)

if __name__ == '__main__':

test_corrcoef()

test_timing()

我的Macbook对100x2500矩阵的性能corrcoef takes 0.06608104705810547 seconds loops=10

corrcoef_loop takes 7.585600137710571 seconds loops=10