Pythonプログラミング(ステップ7・平均値の検定とt分布)
このページでは、t検定を行う際に必要となるt統計量とその分布関数について考える。
平均値の検定
平均
言うまでもなく、母平均はサンプルの平均
このような設定の下で、得られたサンプル平均がどれくらい予想される値として理解可能なのかを吟味する(平均値の検定を行う)には、どのようにすればよいだろうか。
t統計量
さて、正規分布に従うXから、
ここで X の標本平均
そこで、(1)式のt統計量の定義で
(1)式では
このt統計量を使うと、比較的サンプル数が少ない場合でも、『母集団の平均が
このt統計量は t分布と呼ばれる分布関数に従うことが知られている:
自由度1のt分布
計算に興味の無い者 は、このセクションは読み飛ばして構わない。 一方で、計算力の腕試しをしたい諸君は検算してみること。
デルタ関数は、本来の意味の「函数」ではなくて、
もう少し詳しくは こちらのページを参照。
デルタ関数を使ってt分布関数を表現する
では、t統計量はどのような確率分布をするのか、その様子を見るために、最も自由度の小さな場合(
デルタ関数を使うと、
積分を極座標に変換し、動径方向と角度方向の積分に分ける
このままでは見通しが悪いので、積分を極座標
結果をまとめて1自由度のt分布関数を得る
以上をまとめると、2つの確率変数に対するt分布関数(自由度1)は
式(3)に
t統計量の分布は母分散に依らない
変数が
標本平均と分散を確率変数と見なして計算する
母平均を
不偏分散を
このように、ガウス分布するランダムな変数から出発してt分布関数を求めるのは結構な骨が折れるので、
通常は、標本平均
平均値のt検定
30cmのものさしを作っている工場があったとする。
製品の中からランダムに
製品の誤差が正規分布すると考えるのは自然であるから、t統計量を使った検定を行うための前提を満たしている。
サンプルから得た長さ
自由度
練習:ものさしチェッカー
ものさしの実測値を5件入力して、「製品の平均の長さが30cmである」かどうかを有意水準5%で判定するコードを作成しなさい。
ヒント
コードの先頭部分に
import math def tdist(t,nu): return math.gamma((nu+1)/2)/(math.sqrt(nu*math.pi)*math.gamma(nu/2)) * (1+t**2/nu)**(-(nu+1)/2)
と記述することで、自由度tdist(t,nu)
によって呼び出せるようになる。
積分の計算には台形法を用いよ。 DE法(半無限区間)を使うと、さらに精度を高めることができるはずである。
あるいは、SciPyライブラリのscipy.stats.tを使えば、t分布の累積分布を簡単に求めることができる。 以下は、t分布の累積密度関数をプロットするコードの例である:
# coding: utf-8 import math import numpy as np from scipy.stats import t import matplotlib.pyplot as plt nu = 3 xaxis = np.arange(-5, 5, 0.01) tcdf = t.cdf(xaxis,nu) plt.plot(xaxis,tcdf,'-', linewidth=1.0, label="nu="+str(nu)) plt.title("CUMULATIVE T-DISTRIBUTION") plt.xlim(-5,5) plt.ylim(0,1) plt.legend() plt.xlabel('x') plt.ylabel('Cumulative probability') plt.grid(True) plt.show()
平均値の信頼区間
正規分布に従う母集団 からの
その様子を推定するには、
t統計量
平均値の例えば95%信頼区間を見積もるには、t分布密度の両端の面積が 0.025 に相当する部分を除いた、中央の範囲を取ればよい。
その境い目となる
上の議論では、t分布は0を中心に対称であることを使っている。
自由度が大きい場合はt分布は正規分布に漸近する。その場合、95%の信頼区間を与える
データを与えて、その平均値の95%の信頼区間を出力するPythonコードの例を以下に示す。
t分布の計算にはSciPyライブラリを用いている。
確率からt値を逆算するために ppf()
関数が用意されている。
# coding: utf-8 import math import numpy as np from scipy.stats import t import matplotlib.pyplot as plt X=[29.978, 29.933, 29.9245, 30.0944, 30.0981] n = len(X) nu = n-1 x = np.array(X) s = np.mean(x) u2 = np.var(x,ddof=1) u = np.sqrt(u2) t_lower = t.ppf(0.025, nu) t_upper = t.ppf(1-0.025, nu) m_lower = - t_upper * u/n + s m_upper= - t_lower * u/n + s print(m_lower, m_upper) xaxis_low = np.arange(-5, t_lower, 0.01) xaxis_mid = np.arange(t_lower, t_upper, 0.01) xaxis_upper = np.arange(t_upper, 5, 0.01) tdist_low = t.pdf(xaxis_low,nu) tdist_mid = t.pdf(xaxis_mid,nu) tdist_upper = t.pdf(xaxis_upper,nu) fig = plt.figure() ax = fig.add_subplot(111) ax.grid() ax.fill_between(xaxis_low,tdist_low,'-', color='red') ax.plot(xaxis_mid,tdist_mid,'-', color='blue', linewidth=1.0, label="nu="+str(nu)) ax.fill_between(xaxis_upper,tdist_upper,'-', color='red', linewidth=1.0) ax.set_title("T-DISTRIBUTION") ax.set_xlim(-5,5) ax.set_yticks([0, 0.1, 0.2, 0.3, 0.4]) ax.set_yticklabels(["0", "0.1", "0.2", "0.3", "0.4"]) ax.legend() ax.set_xlabel('x') ax.set_ylabel('p(x)') plt.show()
コード実行すると、t分布の95%の信頼区間を青、その両端が赤でプロットされる。