Pythonプログラミング(ステップ7・平均値の検定とt分布)

このページでは、t検定を行う際に必要となるt統計量とその分布関数について考える。

平均値の検定

平均 μ と分散 σ2 は不明であるものの、正規分布に従うことはわかっている母集団 X から、n個のデータを取得したとしよう。 i番目のデータを Xi と表記することにする。 このサンプルから、母平均がある値(例えば0)であると見なしてよいかどうか、統計的に吟味したいとしよう。

言うまでもなく、母平均はサンプルの平均 X¯=1ni=1nXi によって推定できるが、nが十分に大きくない限り、X¯の値は無視できない程度の確率的な変動を伴うはずだ。 そして、その変動は正規分布 X¯N(μ,σ2/n) に従う。しかしながら、母平均 μ と母分散 σ2 は不知である。

このような設定の下で、得られたサンプル平均がどれくらい予想される値として理解可能なのかを吟味する(平均値の検定を行う)には、どのようにすればよいだろうか。

t統計量

Yを正規分布に従う確率変数とする。 観測によってその不偏分散がU2と推定されているとき、 (1)T=YmU で定義される量を t統計量と呼ぶ。 このt統計量は、Uで規格化されているため、母分散の大小には依らず、Ymからの系統的な「ずれ」に対する統一的な指標となることが期待できる。

さて、正規分布に従うXから、n回の観測によって X1,X2,,Xn のサンプルが得られたとすると、 Xの不偏分散は UX2=1n1i=1n(XiX¯)2 で与えられる。

ここで X の標本平均 X¯=1ni=1nXi を確率変数と見なせば、標本平均の不変分散 U2 は Xの不偏分散 UX2を使って U2=UX2/n となる。

そこで、(1)式のt統計量の定義で Y=X¯ とおけば、平均についてのt統計量 (2)T=YmU=X¯mUX/n を構成できる。

(1)式では X¯ の分散の大小をちょうど相殺するようにこれをUX/nで割った格好になっている。 ここで、分子のX¯も、分母のUXも、共にサンプル毎に確率的に変動するので、全体としてTは正規分布とはならない (ただしnが十分大きければ、不偏分散はσ2/nに近づくから、tは平均0分散1の正規分布に漸近する)。

このt統計量を使うと、比較的サンプル数が少ない場合でも、『母集団の平均が m である』という帰無仮説を設定して、 分散の大小には影響されることなく検定することができる。

このt統計量は t分布と呼ばれる分布関数に従うことが知られている: (3)f(t)=Γ((ν+1)/2)νπΓ(ν/2)(1+t2/ν)(ν+1)/2 ここで、ν=n1は分布の自由度と呼ばれる。Γ()ガンマ関数である。

自由度1のt分布

計算に興味の無い者 は、このセクションは読み飛ばして構わない。 一方で、計算力の腕試しをしたい諸君は検算してみること。

デルタ関数は、本来の意味の「函数」ではなくて、 δ(ts)f(s)ds=f(t) となるような性質を持った分布関数のことである。

もう少し詳しくは こちらのページを参照。

デルタ関数を使ってt分布関数を表現する

では、t統計量はどのような確率分布をするのか、その様子を見るために、最も自由度の小さな場合(ν=1)について考えよう。 平均は0、分散σ2のガウス分布に従う独立な確率変数 XY から(ここでのYは前節の説明とは関係ない)、(2)式に従ってt統計量を「合成」すると T=(X+Y)/2(X(X+Y)2)2+(Y(X+Y)2)2/2 となる。 このTが従う分布関数 f(t) を具体的に調べてみよう。

デルタ関数を使うと、f(t)f(t)=1(2πσ)2++δ(t(x+y)/2(x(x+y)2)2+(y(x+y)2)2/2)exp[x2+y22σ2]dxdy と表現することができるので、この式から出発して考えてみる。

積分を極座標に変換し、動径方向と角度方向の積分に分ける

このままでは見通しが悪いので、積分を極座標 x=rcosθ,y=rsinθ に変数変換すると f(t)=1(2πσ)202π0+δ(tcosθ+sinθ|cosθsinθ|)exp[r22σ2]rdrdθ =12πσ20+rexp[r22σ2]dr02πδ(tcosθ+sinθ|cosθsinθ|)dθ となる。デルタ関数の中のrは分子と分母でキャンセルして消えることに注意。ここで、rについての積分は 0+rexp[r22σ2]dr=σ2 となる。

θについての積分は、絶対値の中の符号に気をつけながら書き下すと 02πδ(tcosθ+sinθ|cosθsinθ|)dθ=0π/4δ(tcosθ+sinθcosθsinθ)dθ+π/43π/4δ(tcosθ+sinθsinθcosθ)dθ+3π/42πδ(tcosθ+sinθcosθsinθ)dθ となる。ここで、変数変換 s=tcosθ+sinθ|cosθsinθ| を施すと dθds=12(1tanθ)21+tan2θ tanθ(t,s)=1+(ts)±1+(ts) であるから(複号同順)、デルタ関数の性質、および(1tanθ(t,0))21+tan2θ(t,0)=21+t2を使って 02πδ(tcosθ+sinθ|cosθsinθ|)dθ=12t1tδ(s)(1tanθ(t,s))21+tan2θ(t,s)ds+12tt+δ(s)(1tanθ(t,s))21+tan2θ(t,s)ds12t+t1δ(s)(1tanθ(t,s))21+tan2θ(t,s)ds =21+t2 を得る(この部分は実質的に積分計算は不要)。 計算の過程で、あるtに対して、sの積分区間が0をよぎるのは、上記の3つ項のうちの2つであること、 そして、全体にマイナス符号が付いている積分は、積分範囲が「大きいほうから小さいほう」の向きである点に注意。

結果をまとめて1自由度のt分布関数を得る

以上をまとめると、2つの確率変数に対するt分布関数(自由度1)は f(t)=1π(1+t2) であることがわかった。 このように、自由度が小さい場合、t分布と正規分布はかなり異なることに注意が必要である。

式(3)にν=1を代入すると、Γ(1)=1,Γ(1/2)=π であるから、 11ππ(1+t2/1)(1+1)/2=1π(1+t2) で、確かに上記の計算結果と一致することが確かめられる。

t統計量の分布は母分散に依らない

変数がn個(自由度がn1)の場合でも同様に、極座標で積分すると(動径成分をr、角度成分をΘ、角度変数のヤコビ行列式をJ(Θ)と表記すると) f(t)=1(2πσ)n0+rn1exp[r22σ2]drδ(tG(Θ)))J(Θ)dΘ のようにrΘの積分を分けることができて、 0+rn1exp[r22σ2]dr=2n21Γ(n/2)σn であることを使うと f(t)=Γ(n/2)2πnδ(tG(Θ)))J(Θ)dΘ となって、分布関数はσに依存しないことがわかる (t統計量は標本平均を不偏分散の平方根で規格化した量なので、これは当然である)。

標本平均と分散を確率変数と見なして計算する

母平均をμ、母分散をσ2とすると n個のデータから得られる標本平均は平均μ、分散σ2/nの正規分布に従う。

不偏分散をU2とすると、(n1)U2/σ2は自由度n1のカイ二乗分布に従う。

このように、ガウス分布するランダムな変数から出発してt分布関数を求めるのは結構な骨が折れるので、 通常は、標本平均X¯と不偏分散U2をそれぞれ独立なランダムな変数と考えて計算を行う。例えば上記の自由度1の場合は、

f(t)=++δ(twu/2)δ(u(xw)2+(yw)2)δ(w(x+y)/2)ex2+y22σ2(2πσ)2dxdy=+dw0+duδ(twu/2)g(u)h(w) のように予め式を変形しておけば、標本平均の分布関数h(w)と不偏分散の平方根の分布関数g(u)から計算を出発することができる。 詳しくは専門書にあたること。

平均値のt検定

30cmのものさしを作っている工場があったとする。 製品の中からランダムに5個のサンプルを抽出し、実際の長さを検査してみたとしよう。 ここで、帰無仮説として「製品(ものさし)の平均の長さは30cmである」を設け、 これを検定してみよう。

製品の誤差が正規分布すると考えるのは自然であるから、t統計量を使った検定を行うための前提を満たしている。

サンプルから得た長さi[cm]から、標本平均 w=15(1+2+3+4+5) と、不偏分散 u2=14[(1w)2+(2w)2+(3w)2+(4w)2+(5w)2] が得られるから、それらを用いて t統計量 t=w30u/5 を求める。得られた t が生じるであろう確率が有意水準よりも小さければ、帰無仮説は棄却される(製品の平均長が30cmである可能性は確率的に小さい)。

自由度νのt分布関数をfν(t)とすると、両側検定の場合、 ±|t|よりも分布の端を占める割合 Fν(t)=|t|fν(s)ds+|t|fν(s)ds=2|t|fν(s)ds を求め、有意水準(例えば0.05)と比較すればよい(t分布関数は中心対称なので、片側分だけ計算して2倍すればよい)。

icon-pc 練習:ものさしチェッカー

ものさしの実測値を5件入力して、「製品の平均の長さが30cmである」かどうかを有意水準5%で判定するコードを作成しなさい。

icon-hint ヒント

コードの先頭部分に

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)

と記述することで、自由度νのt分布関数を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()

平均値の信頼区間

正規分布に従う母集団 からのn個のサンプル X1,X2,,Xn が得られているとき、 母平均は「どの範囲」にあると考えられるだろうか。 大数の法則から、標本平均 X¯=1niXin で 母平均 μ に漸近することが保証されているものの、 nがそれほど大きくなければ X¯μ の周りで確率的に揺らぐ。

その様子を推定するには、 t統計量 T=X¯mU/n を考え、m の範囲を色々と変えてみて、確率密度(t分布関数)の山の中央付近をとればよさそうである。 t検定の帰無仮説を設定する際のmの値を沢山試してみて、帰無仮説が棄却されるようなmの範囲だけ除くわけである。

平均値の例えば95%信頼区間を見積もるには、t分布密度の両端の面積が 0.025 に相当する部分を除いた、中央の範囲を取ればよい。 その境い目となる T の値を、それぞれ T, +T とすれば、 TX¯mU/n+T であり、これから m の信頼区間は X¯TU/nmX¯+TU/n で与えられる。

上の議論では、t分布は0を中心に対称であることを使っている。

自由度が大きい場合はt分布は正規分布に漸近する。その場合、95%の信頼区間を与えるTは、1.96に近づく。

データを与えて、その平均値の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%の信頼区間を青、その両端が赤でプロットされる。