2013-02-22 61 views
46

正規分布を仮定して、信頼区間を計算したいサンプルデータがあります。サンプルデータからの信頼区間を計算する

私はnumpyパッケージとscipyパッケージを見つけてインストールしました。平均と標準偏差(numpy.mean(data))を返すようnumpyになっています。サンプルの信頼区間を得るためのアドバイスは、非常に高く評価されます。

答えて

82
import numpy as np 
import scipy as sp 
import scipy.stats 

def mean_confidence_interval(data, confidence=0.95): 
    a = 1.0*np.array(data) 
    n = len(a) 
    m, se = np.mean(a), scipy.stats.sem(a) 
    h = se * sp.stats.t._ppf((1+confidence)/2., n-1) 
    return m, m-h, m+h 

このように計算できます。

ここ
+1

sp.stats.stderrは推奨されていません。私はsp.stats.semを置き換えました。 – Bmayer0122

+1

'scipy'をインポートしても必ずしもすべてのサブパッケージが自動的にインポートされるわけではありません。サブパッケージ 'scipyをインポートする方がよいでしょう。統計情報を明示的に表示する。 – Vikram

+22

'sp.stats.t._ppf'の「私的な」使用には注意が必要です。私はそれ以上の説明なしでそこにそれで快適ではない。あなたが何をしているのか分からない限り、 'sp.stats.t.ppf'を直接使うのが良いでしょう。 [ソース](https://github.com/scipy/scipy/blob/v0.13.0/scipy/stats/distributions.py#L1474)の素早い検査では、 '_ppf'でスキップされたかなりの量のコードがあります。恐らく良性ですが、おそらく危険な最適化の試みでしょうか? – Russ

6

look-up tableから、あなたの希望する信頼区間のためにz-valueを検索することから始めます。信頼区間はmean +/- z*sigmaで、sigmaはサンプル平均の推定標準偏差で、sigma = s/sqrt(n)で与えられます。sはサンプルデータから計算された標準偏差で、nはサンプルサイズです。

+20

( 'scipy.stats.norm.intervalここでは、適切なオプションが(基本的に)同じ信頼区間を与える例

信頼、loc =平均、スケール=シグマ) ' – Jaime

+0

私はその機能を見ていませんでした。ありがとう! – bogatron

+3

元の質問者は、正規分布が仮定されることを示しましたが、小さな標本集団(N <100程度)の場合、[Student tの分布](http: //en.wikipedia.org/wiki/Student%27s_t-distribution)ではなく、[正規分布](http://en.wikipedia.org/wiki/Standard_normal_table)の代わりに使用されます。シャサンの答えはすでにこれをしている。 – Russ

45

shasanのコードの短縮版、配列aの平均の95%信頼区間を計算する:

import numpy as np, scipy.stats as st 

st.t.interval(0.95, len(a)-1, loc=np.mean(a), scale=st.sem(a)) 

しかしStatsModels' tconfint_meanを使用すると、間違いなくでもよりよいです:

import statsmodels.stats.api as sms 

sms.DescrStatsW(a).tconfint_mean() 

両方のための根本的な仮定は、サンプル(配列a)が未知の標準偏差を有する正規分布から独立して描かれたことである(MathWorldまたはWikipedia)。

サンプルサイズが大きい場合、標本平均は通常分布し、st.norm.interval()(Jaimeのコメントで示唆されているように)を使用して信頼区間を計算できます。しかし、上記の解は、小さいnについても正しい。ここで、st.norm.interval()は、信頼区間が狭すぎる(すなわち、「偽の信頼」)。詳細については、私のanswerと同様の質問を参照してください(Russのコメントの1つ)。

In [9]: a = range(10,14) 

In [10]: mean_confidence_interval(a) 
Out[10]: (11.5, 9.4457397432391215, 13.554260256760879) 

In [11]: st.t.interval(0.95, len(a)-1, loc=np.mean(a), scale=st.sem(a)) 
Out[11]: (9.4457397432391215, 13.554260256760879) 

In [12]: sms.DescrStatsW(a).tconfint_mean() 
Out[12]: (9.4457397432391197, 13.55426025676088) 

そしてst.norm.interval()を使用して最終的には、誤った結果:

In [13]: st.norm.interval(0.95, loc=np.mean(a), scale=st.sem(a)) 
Out[13]: (10.23484868811834, 12.76515131188166) 
+0

95%信頼区間を得るために' st.t.interval(0.05) 'を呼び出すべきだと思います。 – Scimonster

+1

いいえ、 'st.t.interval(0.95)'は95%信頼区間に対して正しいです。[docs](http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats。 t.html)に 'scipy.stats.t'を追加しました。 SciPyの引数alphaの命名は理想的ではないようだ。 –

関連する問題