Entender mejor la prueba T usando un experimento

Question

Considere lo siguiente simple python código:

import numpy as np
import scipy.stats as stats
np.random.seed(1000)#(seed)
for i in range(10):
    X = np.random.normal(0,1,2000000)
    Y = np.random.normal(0,1,2000000)
    print(stats.ttest_ind(X,Y))

Ya que ambos $X$ y $Y$ se muestrean a partir de la misma distribución normal (misma media y misma varianza) Esperaba que los valores p devueltos fueran bastante altos. Sin embargo, en realidad, en mi máquina, obtengo la siguiente salida:

(-0.63504239247902361, 0.52540080331521977)
(0.95766169154615455, 0.33823343778490922)
(0.77495702809122613, 0.43836509803793244)
(1.0889385249318819, 0.27618106157897093)
(0.97293351367380154, 0.33058640601636757)
(1.3556340433654583, 0.17521571107906755)
(0.58921373947134548, 0.55571793482831677)
(0.91391361607063426, 0.36076229350915234)
(0.62066063893057222, 0.5348229679159715)
(2.3510055029223147, 0.018722805219522423)

Fíjese que para $ \alpha =0.05$ el último par de muestras sugeriría que debería rechazar la hipótesis nula, a saber que $ \mu_X = \mu_Y $ .

Como un control de cordura corrí,

stats.ttest_ind(X,X)

que da lugar a (0.0, 1.0) como se esperaba. Además, también he hecho lo siguiente:

import numpy as np
import scipy.stats as stats
np.random.seed(1000)#(seed)
for i in range(10):
    X = np.random.normal(0,1,2000000)
    Y = np.random.normal(10,1,2000000) # Different mean for Y
    print(stats.ttest_ind(X,Y))

que siempre dio un valor p cero, como esperaba.

La pregunta: ¿Qué está pasando aquí? ¿Por qué es este el comportamiento y qué me estoy perdiendo?

Answer 1

Buena pregunta. Supongo que su verdadera pregunta es: ¿por qué obtengo errores de tipo I cuando utilizo 2 millones de muestras? Si reflexiona sobre ello, su observación depende realmente sólo de $\alpha$ (tamaño de la prueba). Para $\alpha = 0.05$ si muestrea su estadística de prueba 100 veces, tendría un error de tipo I 5 veces de media. En su caso, para 10 muestras ha observado un error de tipo I una vez, lo que es posible (con probabilidad $10*0.95^9*0.05 = 0.315$ ). Un experimento más relevante podría ser:

import scipy.stats as stats
np.random.seed(1000)#(seed)
alpha = 0.05
n = 1000
count = 0
for i in range(n):
    X = np.random.normal(0,1,2000)
    Y = np.random.normal(0,1,2000)
    count += 1 if stats.ttest_ind(X,Y)[1] <= alpha else 0
count/float(n)

Salida: .055 (en línea con la elegida $\alpha$ )