我有一个矩阵
datanp.corrcoefimport numpy as np
data = np.array([[0, 1, -1], [0, -1, 1]])
np.corrcoef(data)
现在我还想看看这些系数的 p 值。
np.corrcoefscipy.stats.pearsonrscipy.stats.pearsonr是否有一种快速方法可以计算所有行对的系数和 p 值（例如，到达两个 m × m 矩阵，一个具有相关系数，另一个具有相应的 p 值），而无需必须手动检查所有对？
我今天也遇到了同样的问题。
经过半个小时的谷歌搜索，我在 numpy/scipy 库中找不到任何代码可以帮助我做到这一点。
所以我写了我自己的版本 corrcoef
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 需要 0.06608104705810547 秒循环=10
corrcoef_loop 需要 7.585600137710571 秒循环=10
最简洁的方法可能是
.corrpandasIn [79]:
import pandas as pd
m=np.random.random((6,6))
df=pd.DataFrame(m)
print df.corr()
0 1 2 3 4 5
0 1.000000 -0.282780 0.455210 -0.377936 -0.850840 0.190545
1 -0.282780 1.000000 -0.747979 -0.461637 0.270770 0.008815
2 0.455210 -0.747979 1.000000 -0.137078 -0.683991 0.557390
3 -0.377936 -0.461637 -0.137078 1.000000 0.511070 -0.801614
4 -0.850840 0.270770 -0.683991 0.511070 1.000000 -0.499247
5 0.190545 0.008815 0.557390 -0.801614 -0.499247 1.000000
使用 t 检验获取 p 值：
In [84]:
n=6
r=df.corr()
t=r*np.sqrt((n-2)/(1-r*r))
import scipy.stats as ss
ss.t.cdf(t, n-2)
Out[84]:
array([[ 1. , 0.2935682 , 0.817826 , 0.23004382, 0.01585695,
0.64117917],
[ 0.2935682 , 1. , 0.04363408, 0.17836685, 0.69811422,
0.50661121],
[ 0.817826 , 0.04363408, 1. , 0.39783538, 0.06700715,
0.8747497 ],
[ 0.23004382, 0.17836685, 0.39783538, 1. , 0.84993082,
0.02756579],
[ 0.01585695, 0.69811422, 0.06700715, 0.84993082, 1. ,
0.15667393],
[ 0.64117917, 0.50661121, 0.8747497 , 0.02756579, 0.15667393,
1. ]])
In [85]:
ss.pearsonr(m[:,0], m[:,1])
Out[85]:
(-0.28277983892175751, 0.58713640696703184)
In [86]:
#be careful about the difference of 1-tail test and 2-tail test:
0.58713640696703184/2
Out[86]:
0.2935682034835159 #the value in ss.t.cdf(t, n-2) [0,1] cell
你也可以使用你在OP中提到的
scipy.stats.pearsonrIn [95]:
#returns a list of tuples of (r, p, index1, index2)
import itertools
[ss.pearsonr(m[:,i],m[:,j])+(i, j) for i, j in itertools.product(range(n), range(n))]
Out[95]:
[(1.0, 0.0, 0, 0),
(-0.28277983892175751, 0.58713640696703184, 0, 1),
(0.45521036266021014, 0.36434799921123057, 0, 2),
(-0.3779357902414715, 0.46008763115463419, 0, 3),
(-0.85083961671703368, 0.031713908656676448, 0, 4),
(0.19054495489542525, 0.71764166168348287, 0, 5),
(-0.28277983892175751, 0.58713640696703184, 1, 0),
(1.0, 0.0, 1, 1),
#etc, etc
有点黑客，可能效率低下，但我认为这可能就是您正在寻找的：
import scipy.spatial.distance as dist
import scipy.stats as ss
# Pearson's correlation coefficients
print dist.squareform(dist.pdist(data, lambda x, y: ss.pearsonr(x, y)[0]))
# p-values
print dist.squareform(dist.pdist(data, lambda x, y: ss.pearsonr(x, y)[1]))
Scipy 的 pdist 是一个非常有用的函数，主要用于查找 n 维空间中观测值之间的成对距离。
但它允许用户定义可调用的“距离度量”，可以利用它来执行任何类型的成对操作。结果以压缩距离矩阵形式返回，可以使用 Scipy 的“squareform”函数轻松将其更改为方阵形式。
如果您不必使用 pearson 相关系数，您可以使用 spearman 相关系数，因为它同时返回相关矩阵和 p 值（请注意，前者要求您的数据呈正态分布，而斯皮尔曼相关性是一种非参数测量，因此不假设数据呈正态分布）。示例代码：
from scipy import stats
import numpy as np
data = np.array([[0, 1, -1], [0, -1, 1], [0, 1, -1]])
print 'np.corrcoef:', np.corrcoef(data)
cor, pval = stats.spearmanr(data.T)
print 'stats.spearmanr - cor:\n', cor
print 'stats.spearmanr - pval\n', pval
这与 MATLAB 中的 corrcoef 的性能完全相同：
要使用此功能，您需要安装 pandas 和 scipy。
# Compute correlation correfficients matrix and p-value matrix
# Similar function as corrcoef in MATLAB
# dframe: pandas dataframe
def corrcoef(dframe):
fmatrix = dframe.values
rows, cols = fmatrix.shape
r = np.ones((cols, cols), dtype=float)
p = np.ones((cols, cols), dtype=float)
for i in range(cols):
for j in range(cols):
if i == j:
r_, p_ = 1., 1.
else:
r_, p_ = pearsonr(fmatrix[:,i], fmatrix[:,j])
r[j][i] = r_
p[j][i] = p_
return r, p
这是@CT Zhu 的回答的最小版本。我们不需要
pandasnumpyimport numpy as np
import scipy.stats as ss
def corr_significance_two_sided(cc, nData):
# We will divide by 0 if correlation is exactly 1, but that is no problem
# We would simply set the test statistic to be infinity if it evaluates to NAN
with np.errstate(divide='ignore'):
t = -np.abs(cc) * np.sqrt((nData - 2) / (1 - cc**2))
t[t == np.nan] = np.inf
return ss.t.cdf(t, nData - 2) * 2 # multiply by two to get two-sided p-value
x = np.random.uniform(0, 1, (8, 1000))
cc = np.corrcoef(x)
pVal = corr_significance_two_sided(cc, 1000)
如果有人有类似的问题，但你的矩阵是 pd.DataFrame 对象，我编写了以下代码：
from scipy.stats import pearsonr
def corr_pval(df):
corr_pval_df = pd.DataFrame(index=df.columns, columns=df.columns)
for i in range(len(corr_pval_df.index)):
for c in range(len(corr_pval_df.columns)):
corr_pval_df.iloc[i, c] = pearsonr(df[corr_pval_df.index[i]], df[corr_pval_df.columns[c]])
return corr_pval_df
corr_pval(corr_df)
我有解决办法！ 我需要对大小为 2000x30,000 的数组执行此操作，并且使用上述方法或双循环是不可行的，特别是当我似乎缺少一个明显的解决方案时。因此，我研究了 scipy 的 Pearson Correlation 实现，以弄清楚他们如何计算 p 值，并看看是否可以针对二维数组对其进行优化。在注释中，他们解释说，他们估计 Pearson 相关系数 (r) 的 PDF，并根据该“r”计算两侧 p 值。
假设x和y取自独立正态分布（因此总体相关系数为0），则样本相关系数r的概率密度函数为：
$$ f(r) = rac{\左 ( 1-r^2 右 )^{ rac{n}{2}-2}}{B\left ( rac{1}{2}, rac{n}{2}-1 右）}$$
其中n是样本数，B是beta函数。这有时被称为 r 的精确分布。 对于相关系数为 r 的给定样本，p 值是从相关性为零的总体中抽取的随机样本 x' 和 y' 的 abs(r') 大于或等于 abs(r)
的概率
这很容易应用于二维数组，我很惊讶他们在 np.corrcoef 中没有这个功能。
import numpy as np
from scipy import stats
N, L = 100, 2500
signals = np.random.random((N, L)).astype(np.float64)
R = np.corrcoef(signals)
以下内容直接摘自Scipy的pearsonr
的笔记dist = stats.beta(L/2 - 1, L/2 - 1, loc=-1, scale=2)
P = 2*dist.cdf(-abs(R))
此方法几乎是即时的，大约与 np.corrcoef 一样快。 我也通过比较双循环方式来检查这些值是否正确，并得到了这个。
testR = np.empty_like(R)
testP = np.empty_like(P)
for i, sigA in enumerate(signals):
for j, sigB in enumerate(signals):
testR[i,j], testP[i,j] = pearsonr(sigA, sigB)
print(np.abs(testP - P).mean(), np.abs(testR - R).mean(), np.abs(R - R.T).mean())
#1.3860685982216431e-14 2.8794441450755465e-17 8.459507779608882e-19
#There's probably some float rounding error that causes the discrepancies in the 14th/17th decimal place.