Rastreando las suposiciones hechas por la función ttest_ind() de SciPy

Question

Estoy intentando escribir mi propio código en Python para calcular las estadísticas t y los valores p para pruebas t independientes de una y dos colas. Puedo utilizar la aproximación normal, pero por el momento estoy tratando de utilizar sólo la distribución t. No he tenido éxito en hacer coincidir los resultados de la biblioteca de estadísticas de SciPy en mis datos de prueba. Me vendría bien un par de ojos frescos para ver si estoy cometiendo un error tonto en alguna parte.

Nota, esto no es tanto una pregunta de codificación como un "¿por qué este cálculo no está dando la t-stat correcta?" Doy el código para completar, pero no esperes ningún consejo de software. Sólo ayuda para entender por qué no es correcto.

Mi código:

import numpy as np
import scipy.stats as st

def compute_t_stat(pop1,pop2):

    num1 = pop1.shape[0]; num2 = pop2.shape[0];

    # The formula for t-stat when population variances differ.
    t_stat = (np.mean(pop1) - np.mean(pop2))/np.sqrt( np.var(pop1)/num1 + np.var(pop2)/num2 )

    # ADDED: The Welch-Satterthwaite degrees of freedom.
    df = ((np.var(pop1)/num1 + np.var(pop2)/num2)**(2.0))/(   (np.var(pop1)/num1)**(2.0)/(num1-1) +  (np.var(pop2)/num2)**(2.0)/(num2-1) ) 

    # Am I computing this wrong?
    # It should just come from the CDF like this, right?
    # The extra parameter is the degrees of freedom.

    one_tailed_p_value = 1.0 - st.t.cdf(t_stat,df)
    two_tailed_p_value = 1.0 - ( st.t.cdf(np.abs(t_stat),df) - st.t.cdf(-np.abs(t_stat),df) )    

    # Computing with SciPy's built-ins
    # My results don't match theirs.
    t_ind, p_ind = st.ttest_ind(pop1, pop2)

    return t_stat, one_tailed_p_value, two_tailed_p_value, t_ind, p_ind

Actualización:

Después de leer un poco más sobre la prueba t de Welch, vi que debería utilizar la fórmula Welch-Satterthwaite para calcular los grados de libertad. He actualizado el código anterior para reflejar esto.

Con los nuevos grados de libertad, obtengo un resultado más cercano. Mi valor p de dos caras se aleja aproximadamente 0,008 del de la versión de SciPy... pero sigue siendo un error demasiado grande, así que debo estar haciendo algo incorrecto (o las funciones de distribución de SciPy son muy malas, pero es difícil creer que sólo sean precisas con dos decimales).

Segunda actualización:

Mientras seguía probando cosas, pensé que tal vez la versión de SciPy calcula automáticamente la aproximación Normal a la distribución t cuando los grados de libertad son lo suficientemente altos (aproximadamente > 30). Así que volví a ejecutar mi código usando la distribución Normal en su lugar, y los resultados calculados están realmente más lejos de los de SciPy que cuando uso la distribución t.

Answer 1

Utilizando la función incorporada de SciPy source(), pude ver una impresión del código fuente de la función ttest_ind(). Basado en el código fuente, la función incorporada de SciPy está realizando la prueba t asumiendo que las varianzas de las dos muestras son iguales. No está utilizando los grados de libertad de Welch-Satterthwaite.

Sólo quiero señalar que, fundamentalmente, esta es la razón por la que no se debe confiar únicamente en las funciones de la biblioteca. En mi caso, realmente necesito la prueba t para poblaciones de varianzas desiguales, y el ajuste de los grados de libertad podría ser importante para algunos de los conjuntos de datos más pequeños en los que voy a ejecutar esto. SciPy asume que las varianzas son iguales, pero no indica esta suposición.

Como mencioné en algunos comentarios, la discrepancia entre mi código y el de SciPy es de aproximadamente 0,008 para tamaños de muestra entre 30 y 400, y luego va lentamente a cero para tamaños de muestra mayores. Esto es un efecto del término extra (1/n1 + 1/n2) en el denominador del estadístico t de varianzas iguales. Desde el punto de vista de la precisión, esto es bastante importante, especialmente para tamaños de muestra pequeños. Definitivamente me confirma que tengo que escribir mi propia función. (Posiblemente haya otras bibliotecas de Python mejores, pero esto al menos debería saberse. Francamente, es sorprendente que esto no esté en ningún lugar en la documentación de SciPy para ttest_ind()).