在2D的NumPy数组列上计算成对的t检验。

Question

在2D的NumPy数组列上计算成对的t检验。

4

我想对2D数组列执行成对的t检验。如果不使用itertools，有什么方法可以获得所有成对列的组合？

from scipy import stats
import numpy as np

a = np.random.randn(20,6)

- hippocampus

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- FBruzzesi · Accepted Answer

如果你希望避免使用for循环或列表推导式，可以使用numpy广播实现t检验：

a = np.random.randn(20,6)
n1, n2 = a.shape
# Columns mean and squared std
m = np.mean(a,axis=0)
s2 = np.square(np.std(a,axis=0, ddof=1))

# Compute the test statistic
t = (m[:,np.newaxis] - m[np.newaxis,:])/np.sqrt((s2[:,np.newaxis] + s2[np.newaxis,:])/n1)

# Compute the 2-sided p-value
df = 2*n1 - 2
p = 2 - 2*stats.t.cdf(t,df=df)

相对于朴素列表推导式实现来检查性能：

def t_test(a):
    n1, n2 = a.shape
    m = np.mean(a,axis=0)
    s2 = np.square(np.std(a,axis=0, ddof=1))

    t = (m[:,np.newaxis] - m[np.newaxis,:])/np.sqrt((s2[:,np.newaxis] + s2[np.newaxis,:])/n1)

    df = 2*n1 - 2
    p = 2 - 2*stats.t.cdf(t,df=df)

    return t,p


%timeit t_test(a)
213 µs ± 13 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%timeit [[ (i,j, stats.ttest_ind(a[:,i], a[:,j])) for i in range(n2) if i <j] for j in range(n2)]
4.36 ms ± 139 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

展示了numpy实现的速度大约快了20倍。