numpy 矩阵所有行对的相关系数和p值

bnl4lu3b  于 11个月前  发布在  其他
关注(0)|答案(8)|浏览(156)

我有一个矩阵data,有 m 行和 n 列。我曾经使用np.corrcoef计算所有行对之间的相关系数:

import numpy as np
data = np.array([[0, 1, -1], [0, -1, 1]])
np.corrcoef(data)

字符串
现在我也想看看这些系数的p值。np.corrcoef没有提供这些; scipy.stats.pearsonr有。但是,scipy.stats.pearsonr不接受输入矩阵。
有没有一种快速的方法来计算所有行对的系数和p值(例如,得到两个 m × m 矩阵,一个具有相关系数,另一个具有相应的p值),而不必手动遍历所有行对?

rlcwz9us

rlcwz9us1#

我今天遇到了同样的问题。
在谷歌上搜索了半个小时后,我在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上的性能与100 x2500矩阵相比
corrcoef需要0.06608104705810547秒循环=10
corrcoef_loop需要7.585600137710571秒loops=10

kiayqfof

kiayqfof2#

最简单的方法可能是在pandas中构建方法.corr,以获取r:

In [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.pearsonr

In [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

p4tfgftt

p4tfgftt3#

有点黑客和可能效率低下,但我认为这可能是你正在寻找的:

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's pdist是一个非常有用的函数,主要用于计算n维空间中观测值之间的成对距离。
但它允许用户自定义可调用的“距离度量”,可以利用它来执行任何类型的成对操作。结果以压缩的距离矩阵形式返回,可以使用Scipy的“squareform”函数轻松地将其更改为方阵形式。

wmtdaxz3

wmtdaxz34#

如果你不需要使用pearson correlation coefficient,你可以使用spearman correlation coefficient,因为它返回相关矩阵和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

字符串

n8ghc7c1

n8ghc7c15#

这与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

字符串

abithluo

abithluo6#

这里是@CT Zhu的答案的最小版本。我们不需要pandas,因为相关性可以直接从numpy计算,这应该更快,因为我们不需要转换为一个矩阵的步骤

import 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)

字符串

bxpogfeg

bxpogfeg7#

如果有人有类似的问题,但你的矩阵是一个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)

字符串

qlzsbp2j

qlzsbp2j8#

我有办法!我需要对大小为2000 x30,000的数组执行此操作,使用上述方法或双循环是不可行的,特别是当似乎有一个明显的解决方案我错过了。所以我研究了scipy的Pearson Correlation实现,以了解他们如何计算p值,看看我是否可以优化它为2-d数组。在笔记中,他们解释说,他们估计皮尔逊相关系数(r)的PDF,并从这个“r”计算双侧p值。
假设x和y是独立的正态分布(因此总体相关系数为0),样本相关系数r的概率密度函数为:
$$ f(r) = \frac{\left ( 1-r^2 \right )^{\frac{n}{2}-2}}{B\left ( \frac{1}{2},\frac{n}{2}-1 \right )}$$
其中n是样本数,B是beta函数。这有时被称为r的精确分布。对于具有相关系数r的给定样本,p值是从具有零相关性的总体中抽取的随机样本x'和y'的abs(r ')大于或等于abs(r)的概率。
这很容易应用于2-d数组,我很惊讶他们在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.

相关问题