Śledzenie założeń wykonanych przez funkcję `ttest_ind()` SciPy

Question

Próbuję napisać własny kod w języku Python, aby obliczyć statystyki t i p-wartości dla jednego i dwóch niezależnych testów t tailed. Mogę użyć normalnego przybliżenia, ale na razie próbuję po prostu użyć rozkładu t. Nie udało mi się dopasować wyników statystyk biblioteki SciPy do moich danych testowych. Mogłabym użyć świeżej pary oczu, żeby zobaczyć, czy nie robię gdzieś głupiego błędu.

Uwaga, jest to cross-posted from Cross-Validated, ponieważ trwało to przez jakiś czas bez odpowiedzi, więc pomyślałem, że nie zaszkodzi również uzyskać opinie programistów. Próbuję zrozumieć, czy jest jakiś błąd w algorytmie, który używam, który powinien reprodukować wynik SciPy. Jest to prosty algorytm, dlatego zastanawiające jest, dlaczego nie mogę zlokalizować błędu.

Mój kod:

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 

Aktualizacja:

Po przeczytaniu nieco więcej na temat testu t-Welch, widziałem, że należy przy użyciu wzoru Welch-Satterthwaitea obliczyć stopnie wolność. Zaktualizowałem powyższy kod, aby to odzwierciedlić.

Z nowymi stopniami swobody, uzyskuję bliższy wynik. Moja dwustronna wartość p jest wyłączona o około 0,008 od wersji SciPy ... ale to nadal jest zbyt duży błąd, więc nadal muszę robić coś niepoprawnego (lub funkcje dystrybucji SciPy są bardzo złe, ale trudno w to uwierzyć są one dokładne z dokładnością do 2 miejsc po przecinku).

Druga aktualizacja:

Kontynuując spróbować rzeczy, pomyślałem, że może wersja scipy automatycznie oblicza Normalny przybliżenie do rozkładu t gdy stopnie swobody są wystarczająco wysokie (około> 30). Ponownie przetestowałem mój kod, wykorzystując dystrybucję Normalną, a obliczone wyniki są w rzeczywistości oddalone od SciPy, niż gdy używam dystrybucji t.

Bonus pytanie :) (Więcej teoria statystyczna związana; krępuj się ignorować)

Ponadto, t-statystyka jest ujemna. Właśnie zastanawiałem się, co to oznacza dla jednostronnego testu t-Studenta. Czy to zwykle oznacza, że ​​powinienem badać w kierunku ujemnym w kierunku testu? W moich danych testowych populacja 1 jest grupą kontrolną, która nie otrzymała określonego programu szkolenia zawodowego. Populacja 2 otrzymała go, a zmierzone dane są różnicami płac przed/po leczeniu.

Mam więc powody, by sądzić, że średnia dla populacji 2 będzie większa. Ale z punktu widzenia teorii statystycznej, nie wydaje się słuszne wymyślanie testu w ten sposób. Jak mógłbym sprawdzić (dla testu jednostronnego) w kierunku ujemnym bez polegania na subiektywnej wiedzy o danych? Czy jest to tylko jedna z tych częstych rzeczy, które, choć nie są filozoficznie rygorystyczne, muszą być wykonywane w praktyce?

Answer 1

Używając wbudowanej funkcji SciPy source(), mogłem zobaczyć wydruk kodu źródłowego dla funkcji ttest_ind(). W oparciu o kod źródłowy, wbudowane SciPy wykonuje test t-założenia, zakładając, że wariancje dwóch próbek są równe. Nie używa stopni swobody Welch-Satterthwaite. SciPy zakłada równe wariancje, ale nie podaje tego założenia.

Chcę tylko zwrócić uwagę, że, co najważniejsze, jest to powód, dla którego nie powinieneś po prostu zaufać funkcjom bibliotecznym. W moim przypadku faktycznie potrzebuję testu t dla populacji nierównych wariancji, a regulacja stopni swobody może mieć znaczenie dla niektórych mniejszych zestawów danych, na których to uruchomię.

Jak wspomniałem w niektórych komentarzach, rozbieżność między moim kodem a SciPy's wynosi około 0,008 dla wielkości próbek od 30 do 400, a następnie powoli idzie do zera dla większych rozmiarów próbek. Jest to efekt dodatkowego (1/n1 + 1/n2) terminu w mianowniku statystycznym o równych-wartościach wariancji. Dokładność jest bardzo ważna, szczególnie w przypadku małych próbek. Zdecydowanie potwierdza to, że muszę napisać własną funkcję. (Być może istnieją inne, lepsze biblioteki Pythona, ale przynajmniej powinno to być znane.) Szczerze mówiąc, jest to zaskakujące, że nie jest to nigdzie z przodu i centrum w dokumentacji SciPy dla ttest_ind()).

Answer 2

Wygląda na to, że zapomniałeś ** 2 do licznika twojego df. Stopnie swobody Welch-Satterthwaite.

df = (np.var(pop1)/num1 + np.var(pop2)/num2)/( (np.var(pop1)/num1)**(2.0)/(num1-1) + (np.var(pop2)/num2)**(2.0)/(num2-1)) 

powinno być:

df = (np.var(pop1)/num1 + np.var(pop2)/num2)**2/( (np.var(pop1)/num1)**(2.0)/(num1-1) + (np.var(pop2)/num2)**(2.0)/(num2-1)) 

Answer 3

Nie jesteś obliczania wariancji przykładowy, ale zamiast używasz wariancji populacji. Wariancja próbki dzieli się przez n-1, zamiast n. np.var ma opcjonalny argument o nazwie ddof z powodów podobnych do tego.

To powinno dać swój oczekiwany rezultat:

import numpy as np 
import scipy.stats as st 

def compute_t_stat(pop1,pop2): 

    num1 = pop1.shape[0] 
    num2 = pop2.shape[0]; 
    var1 = np.var(pop1, ddof=1) 
    var2 = np.var(pop2, ddof=1) 

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

    # ADDED: The Welch-Satterthwaite degrees of freedom. 
    df = ((var1/num1 + var2/num2)**(2.0))/((var1/num1)**(2.0)/(num1-1) + (var2/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 

PS: scipy jest open source i wdrażane głównie z Pythona. Mogłeś sprawdzić kod źródłowy dla ttest_ind i samemu odkryć swój błąd.

Po stronie bonusowej: Ty nie podaję numeru decydując się na stronę testu pojedynczego ogniska, patrząc na swoją wartość t. Decydujesz o tym wcześniej ze swoją hipotezą. Jeśli twoja hipoteza zerowa jest taka, że ​​środki są równe i twoja alternatywna hipoteza jest taka, że ​​druga średnia jest większa, twój ogon powinien znajdować się po lewej (ujemnej) stronie. Ponieważ wystarczająco małe (ujemne) wartości twojej wartości t wskażą, że hipoteza alternatywna jest bardziej prawdopodobna, niż hipoteza zerowa.