Melacak asumsi yang dibuat oleh fungsi ttest_ind () SciPy

Saya mencoba untuk menulis kode Python saya sendiri untuk menghitung t-statistik dan nilai-p untuk satu dan dua uji t independen. Saya dapat menggunakan perkiraan normal, tetapi untuk saat ini saya hanya mencoba menggunakan distribusi-t. Saya tidak berhasil dalam mencocokkan hasil perpustakaan statistik SciPy pada data pengujian saya. Saya bisa menggunakan sepasang mata baru untuk melihat apakah saya hanya membuat kesalahan bodoh di suatu tempat.

Catatan, ini bukan pertanyaan coding karena ini adalah "mengapa perhitungan ini tidak menghasilkan t-stat yang tepat?" Saya memberikan kode untuk kelengkapan, tetapi jangan mengharapkan saran perangkat lunak. Hanya membantu memahami mengapa ini tidak benar.

Kode saya:

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

Memperbarui:

Setelah membaca sedikit lebih banyak tentang uji-Welch, saya melihat bahwa saya harus menggunakan rumus Welch-Satterthwaite untuk menghitung derajat kebebasan. Saya memperbarui kode di atas untuk mencerminkan ini.

Dengan derajat kebebasan baru, saya mendapatkan hasil yang lebih dekat. Nilai dua sisi saya mati sekitar 0,008 dari versi SciPy ... tapi ini masih merupakan kesalahan yang terlalu besar sehingga saya masih harus melakukan sesuatu yang salah (atau fungsi distribusi SciPy sangat buruk, tetapi sulit untuk percaya mereka hanya akurat di 2 tempat desimal).

Pembaruan kedua:

Sambil terus mencoba berbagai hal, saya pikir mungkin versi SciPy secara otomatis menghitung perkiraan Normal ke distribusi-t ketika derajat kebebasannya cukup tinggi (kira-kira> 30). Jadi saya kembali menjalankan kode saya menggunakan distribusi Normal, dan hasil yang dihitung sebenarnya jauh dari SciPy daripada ketika saya menggunakan distribusi-t.

Dengan menggunakan sumber fungsi bawaan SciPy (), saya bisa melihat cetakan kode sumber untuk fungsi ttest_ind (). Berdasarkan kode sumber, SciPy built-in melakukan uji-t dengan asumsi bahwa varians dari dua sampel adalah sama. Itu tidak menggunakan derajat kebebasan Welch-Satterthwaite.

Saya hanya ingin menunjukkan bahwa, yang terpenting, mengapa Anda tidak hanya mempercayai fungsi perpustakaan. Dalam kasus saya, saya benar-benar membutuhkan uji-t untuk populasi dengan varian yang tidak sama, dan tingkat penyesuaian kebebasan mungkin penting untuk beberapa set data yang lebih kecil yang akan saya gunakan. SciPy mengasumsikan varian yang sama tetapi tidak menyatakan asumsi ini.

Seperti yang saya sebutkan di beberapa komentar, perbedaan antara kode saya dan SciPy adalah sekitar 0,008 untuk ukuran sampel antara 30 dan 400, dan kemudian perlahan-lahan beralih ke nol untuk ukuran sampel yang lebih besar. Ini adalah efek dari istilah ekstra (1 / n1 + 1 / n2) dalam penyebut statistik t-varian yang sama. Dari segi akurasi, ini cukup penting, terutama untuk ukuran sampel kecil. Jelas menegaskan kepada saya bahwa saya perlu menulis fungsi saya sendiri. (Mungkin ada pustaka Python lain yang lebih baik, tetapi ini setidaknya harus diketahui. Sejujurnya, ini tidak mengejutkan di mana pun di bagian depan dan tengah dalam dokumentasi SciPy untuk ttest_ind ()).