Retrouver les hypothèses de la fonction ttest_ind () de SciPy

8

J'essaie d'écrire mon propre code Python pour calculer les statistiques t et les valeurs p pour des tests t indépendants à une et deux queues. Je peux utiliser l'approximation normale, mais pour le moment, j'essaie simplement d'utiliser la distribution t. Je n'ai pas réussi à faire correspondre les résultats de la bibliothèque de statistiques de SciPy sur mes données de test. Je pourrais utiliser une nouvelle paire d'yeux pour voir si je fais juste une erreur stupide quelque part.

Remarque, ce n'est pas tant une question de codage que c'est un "pourquoi ce calcul ne donne-t-il pas le bon t-stat?" Je donne le code pour l'exhaustivité, mais n'attendez aucun conseil logiciel. Aidez simplement à comprendre pourquoi ce n'est pas bien.

Mon code:

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

Mise à jour:

Après avoir lu un peu plus sur le test t de Welch, j'ai vu que je devrais utiliser la formule de Welch-Satterthwaite pour calculer les degrés de liberté. J'ai mis à jour le code ci-dessus pour refléter cela.

Avec les nouveaux degrés de liberté, j'obtiens un résultat plus proche. Ma valeur p bilatérale est désactivée d'environ 0,008 par rapport à la version SciPy ... mais c'est toujours une erreur beaucoup trop importante donc je dois encore faire quelque chose de incorrect (ou les fonctions de distribution SciPy sont très mauvaises, mais c'est difficile à croire ils ne sont précis qu'à 2 décimales près).

Deuxième mise à jour:

Tout en continuant à essayer, je pensais que la version de SciPy calcule automatiquement l'approximation normale de la distribution t lorsque les degrés de liberté sont suffisamment élevés (environ> 30). J'ai donc réexécuté mon code en utilisant la distribution normale à la place, et les résultats calculés sont en fait plus éloignés de SciPy que lorsque j'utilise la distribution en t.

ely
la source
Peut-être que SciPy calcule le test t de Welch - La documentation de SciPy ne spécifie pas ...
Cyan
La formule que j'utilise dans mon calcul est la même que la statistique t de Welch. À ma connaissance, c'est la chose «standard» à faire lorsque les tailles d'échantillon et les variances de population peuvent être différentes, n'est-ce pas?
le
4
N'avez-vous pas besoin de prendre le carré du numérateur (actuel) dans le calcul des degrés de liberté? De plus, avec pratiquement aucun changement de code, il existe des moyens beaucoup plus sûrs de calculer les valeurs . La façon dont il est actuellement implémenté est extrêmement susceptible d' erreurs massives en raison de l'annulation . p
Cardinal
4
( 1 ) Consultez la documentation de numpy.var. La version que j'ai vue semble indiquer que l'estimation MLE est calculée par défaut au lieu de l'estimation non biaisée. Pour obtenir l'estimation impartiale, il faut l'appeler avec l'option ddof=1. ( 2 ) Pour la queue supérieure -valeur, utiliser la symétrie de la -distribution, à savoir, et ( 3 ) pour les deux-tailed -valeur, faire quelque chose de similaire: . ptone_tailed_p_value = st.t.cdf(-t_stat,df)ptwo_tailed_p_value = 2*st.t.cdf(-np.abs(t_stat),df)
Cardinal
2
Je ne pense pas que ce soit si insignifiant, dans le sens où il y a souvent un écart considérable entre avoir une formule mathématique pour quelque chose à portée de main et connaître un moyen sûr et efficace de le calculer. C'est l'une de ces choses où il est agréable d'avoir un large corpus de connaissances déjà disponible, car il faudrait une éternité virtuelle pour apprendre de telles astuces, une par une, tout seul. :)
Cardinal

Réponses:

4

En utilisant la fonction intégrée SciPy source (), j'ai pu voir une impression du code source de la fonction ttest_ind (). Sur la base du code source, le SciPy intégré effectue le test t en supposant que les variances des deux échantillons sont égales. Il n'utilise pas les degrés de liberté Welch-Satterthwaite.

Je veux juste souligner que, c'est surtout pourquoi vous ne devriez pas simplement faire confiance aux fonctions de bibliothèque. Dans mon cas, j'ai en fait besoin du test t pour les populations de variances inégales, et les degrés d'ajustement de liberté pourraient avoir de l'importance pour certains des plus petits ensembles de données sur lesquels je vais exécuter cela. SciPy suppose des variances égales mais n'énonce pas cette hypothèse.

Comme je l'ai mentionné dans certains commentaires, l'écart entre mon code et SciPy est d'environ 0,008 pour les tailles d'échantillon comprises entre 30 et 400, puis passe lentement à zéro pour les tailles d'échantillon plus grandes. Il s'agit d'un effet du terme supplémentaire (1 / n1 + 1 / n2) dans le dénominateur de la statistique t à variances égales. En termes de précision, cela est assez important, en particulier pour les petits échantillons. Cela me confirme définitivement que j'ai besoin d'écrire ma propre fonction. (Il existe peut-être d'autres meilleures bibliothèques Python, mais cela devrait au moins être connu. Franchement, il est surprenant que ce ne soit nulle part au centre de la documentation SciPy pour ttest_ind ()).

ely
la source
3
Il semble que cela soit désormais correctement implémenté à partir de Scipy 0.11.0 via un paramètre facultatif pour spécifier le test t de Welch: docs.scipy.org/doc/scipy/reference/generated/…
Abhijit Rao