Compute pairwise t-test on 2d numpy array columns

user_13038520 0
Viewed 743

I want to perform pairwise t-tests on 2D array columns. What is a way to get all the pairwise combinations of columns if not with itertools?

from scipy import stats
import numpy as np

a = np.random.randn(20,6)
python-3.x numpy arrays python
1 Answers

If you want to avoid using for loops or list comprehensions you can implement the t-test using numpy broadcasting:

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)

Checking performance with respect to naive list comprehension implementation:

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)

Showing that the numpy implementation is about 20x faster.

user_124115360
Related