Вычислить попарный t-критерий для столбцов массива 2d numpy

Question

Я хочу выполнить попарные t-тесты для столбцов 2D-массива. Как можно получить все парные комбинации столбцов, если не с помощью itertools?

from scipy import stats
import numpy as np
a= np.random.randn(20,6)

Answer 1

Если вы хотите избежать использования для циклов или списков, вы можете реализовать t-тест с использованием numpy широковещания:

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 раз быстрее.