8

片側および両側の独立したt検定のt統計量とp値を計算する独自のPythonコードを作成しようとしています。正規近似を使用できますが、今のところ、t分布を使用しようとしています。SciPyの統計ライブラリの結果をテストデータと照合することに失敗しました。私は新鮮な目を使って、どこかでばかげた間違いをしているのかどうかを確認することができました。

注：これは、しばらくの間応答がないため、Cross-Validatedからクロスポストされているため、ソフトウェア開発者の意見を得るのも悪くないと思いました。使用しているアルゴリズムにエラーがあるかどうかを理解しようとしています。これにより、SciPyの結果が再現されます。これは単純なアルゴリズムなので、間違いを見つけられない理由は不可解です。

私のコード：

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

アップデート：

ウェルチのt検定についてもう少し読んだ後、自由度を計算するためにウェルチ-サタスウェイトの公式を使用する必要があることがわかりました。これを反映するために上記のコードを更新しました。

新しい自由度で、私はより近い結果を得ることができます。私の両側p値はSciPyバージョンから約0.008ずれています...しかし、これはまだ非常に大きなエラーなので、私はまだ何か間違ったことをしている必要があります（またはSciPy分布関数は非常に悪いですが、信じがたいです小数点以下2桁までしか正確ではありません）。

2番目の更新：

試してみながら、SciPyのバージョンでは、自由度が十分に高い場合（約> 30）に、t分布の正規近似が自動的に計算されるのではないかと思いました。そのため、代わりに正規分布を使用してコードを再実行しました。計算結果は、実際には、t分布を使用した場合よりもSciPyから遠く離れています。

ボーナスの質問: )（より統計理論に関連しています;無視してください）

また、t統計量は負です。これが片側t検定にとって何を意味するのか疑問に思っていました。これは通常、テストのために負の軸方向を見る必要があることを意味しますか？私のテストデータでは、母集団1は、特定の雇用訓練プログラムを受けていない対照群です。人口2はそれを受け取りました、そして、測定されたデータは治療の前後の賃金の違いです。

したがって、母集団2の平均が大きくなると考える理由がいくつかあります。しかし、統計理論の観点からは、この方法でテストを作成することは正しくないようです。データに関する主観的な知識に頼らずに、（片側テストの場合）負の方向にチェックすることをどのようにして知ることができましたか？それとも、これは、哲学的に厳密ではありませんが、実際に行う必要がある頻度主義的なことの1つにすぎませんか？

4

3 に答える 3

9

SciPyの組み込み関数を使用することsource()で、関数のソースコードのプリントアウトを見ることができましたttest_ind()。ソースコードに基づいて、SciPyビルトインは、2つのサンプルの分散が等しいと仮定してt検定を実行しています。ウェルチ・サタスウェイト自由度を使用していません。SciPyは等しい分散を想定していますが、この想定を述べていません。

重要なのは、これがライブラリ関数を信頼するだけではいけない理由です。私の場合、分散が等しくない母集団のt検定が実際に必要であり、これを実行する小さなデータセットのいくつかでは自由度の調整が重要になる可能性があります。

いくつかのコメントで述べたように、私のコードとSciPyの不一致は、サンプルサイズが30〜400の場合は約0.008であり、サンプルサイズが大きい場合はゆっくりとゼロになります。これは、等分散t統計分母の余分な（1 / n1 + 1 / n2）項の効果です。精度に関しては、これは特にサンプルサイズが小さい場合に非常に重要です。自分で関数を書く必要があることは間違いありません。（おそらく他のより良いPythonライブラリがありますが、これは少なくとも知っておく必要があります。率直に言って、これがSciPyドキュメントの中心にないのは驚くべきことですttest_ind()）。

于 2012-04-12T21:07:40.920 に答える
2

サンプル分散を計算していませんが、代わりに母分散を使用しています。サンプルの分散はn-1、ではなく、で除算されnます。これと同様の理由でnp.var呼び出されるオプションの引数があります。ddof

これにより、期待どおりの結果が得られるはずです。

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はオープンソースであり、ほとんどがPythonで実装されています。ソースコードをチェックしttest_indて、自分で間違いを見つけたかもしれません。

ボーナス側の場合：t値を見て、片側検定の側を決定することはありません。あなたはあなたの仮説でそれを前もって決定します。帰無仮説が平均が等しいというものであり、対立仮説が2番目の平均が大きいというものである場合、裾は左側（負）にある必要があります。t値の値が十分に小さい（負の）場合は、帰無仮説ではなく対立仮説が真である可能性が高いことを示しているためです。

于 2012-04-06T04:33:10.787 に答える
0

dfの分子に**2を忘れたようです。ウェルチ-サタスウェイト自由度。

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) )

する必要があります：

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) )
于 2012-04-06T03:43:35.123 に答える