Python: SciPyを使った仮説検定

Python 統計学

前回はSciPyを使った推定をとりあげましたが、今回は正規分布に関する仮説検定をSciPyを使って行います。

サンプルには前回同様、Irisデータセットを使いますので、UCIのWebサイトからダウンロードして品種ごと(setosa, virginica, versicolor)に分けておきます。それぞれの標本数は50です。

import pandas as pd
import math
from scipy import stats

iris = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data', 
                 header=None, names=['sepal_length', 'sepal_width', 'petal_length', 'petal_width', 'class'])

setosa = iris[iris['class'] == 'Iris-setosa'] # Iris setosa(ヒオウギアヤメ)
virginica = iris[iris['class'] == 'Iris-virginica'] # Iris virginica
versicolor = iris[iris['class'] == 'Iris-versicolor'] # Iris versicolor

また、赤本の練習問題からも一部抜粋しております。

統計学入門 (基礎統計学?)

統計学入門 (基礎統計学?)

仮説検定のプロセス

最初に仮説検定のプロセスを確認します。

  1. 帰無仮説と対立仮説を立てる
    例えば、setosaのがく片の長さ(sepal_width)は母平均は3.5であるという仮設(帰無仮説)を立てます。ここでは H_0: \mu_0=3.5 とします。
    対して、平均3.5ではないという仮設(対立仮説)を立てます。同じく H_1: \mu_0 \neq 3.5 とします。2つの仮設で大小は関係ない(ただ異なっていることを示す)ため、両側検定となります。

  2. 有意水準を決める
    帰無仮説を棄却するに足るための水準を決めます。
    例えば、両側検定で有意水準95%とした場合、[0.025, 0.975]が H_0 を支持する区間(採択域)となり、この区間から外れた場合に H_0 が棄却(つまり H_1 が支持)される。

  3. 統計量を計算する
    例えば、正規分布の平均であればz値を計算することで、標準正規分布上にプロットできるようになる。
    そうすると「何%の確率で \mu_0 が得られるのか」というのが確率密度が計算できる。

  4. 有意水準と照らして判定する
    得られた確率が採択域から外れていれば H_0 を棄却( H_1 を支持)、採択域に入っていれば H_0 を棄却しないと判定する。

母平均の仮説検定

大標本(n>30)の場合は正規分布、小標本(n≦30)の場合はt分布に従うと仮定します。ただし、nが大きくなればt分布は正規分布に近づくため、nに関わらず小標本の求め方でOKだと思います。

大標本 (n>30)

setosaのがく片の長さ(sepal_width)の母平均について、帰無仮説 H_0: \mu_0=3.5 、対立仮説 H_1: \mu_0 \neq 3.5 とし、有意水準95%で仮説検定します。

z統計量を計算します。未知の母分散 \sigma^2 については、大標本のため標本不偏分散 s^2 で推定しています。

$$ z=\frac{\bar{X} - \mu_0}{\sqrt{\frac{\sigma^2}{n}}}=\frac{\bar{X} - \mu_0}{\sqrt{\frac{s^2}{n}}} $$

z値として-1.52が得られましたので、標準正規分布における累積確率([-∞, -1.52]の区間に入る確率)を累積分布関数cdfで計算します。
結果は0.0640でしたので、95%有意水準の採択域[0.025, 0.975]に入っているため、 H_0 は棄却できない、と判定されます。

mu_0 = 3.5
z = (setosa['sepal_width'].mean() - mu_0) / math.sqrt(setosa['sepal_width'].var() / len(setosa))
print(z) # -1.5217596660094042

stats.norm.cdf(z)
# 0.06403465560380407

小標本 (n≦30)

小標本の場合、正規分布ではなくt分布で推定される。
そのため、まずは以下でt統計量を計算する。そして、自由度n-1のt分布上で累積確率を計算する。

$$ t=\frac{\bar{X} - \mu_0}{\sqrt{\frac{s^2}{n}}} $$

n=10として、同じくsetosaのがく片の長さの母平均について、帰無仮説 H_0: \mu_0=3.5 、対立仮説 H_1: \mu_0 \neq 3.5 において、有意水準95%で仮説検定した実装が以下です。

SciPyではttest_1samp()に、データと帰無仮説(popmean)を渡すと、1発で計算できます。
返り値はタプルで、1つ目がt統計量、2つ目がp値となっています。p値は、帰無仮説が成り立つと仮定した場合に、t統計量以上の値(ここでは|-1.96|=1.96以上の値)が得られる確率を表します。p値は0.082となっていますので、95%有意水準に照らして H_0 は棄却できないという判定になります。

n = 10
mu_0 = 3.5

stats.ttest_1samp(setosa.iloc[:n]['sepal_width'], popmean=mu_0)
# Ttest_1sampResult(statistic=-1.9562349357055535, pvalue=0.08213989920574083)

母平均の差の仮説検定

2つの標本XとYの平均の差について仮説検定します。等分散を仮定できるかどうかで、求め方が異なるので注意が必要。

ここでは赤本の練習問題12.2のデータを使います。同年代・同学歴・同職種の男女各10人の賃金(万円)を表すを架空の例です。
帰無仮説 H_0: \mu_X-\mu_Y=0 、対立仮説 H_1: \mu_X-\mu_Y \neq 0 とし、有意水準95%で検定します。

payment_male = [15.4, 18.3, 16.5, 17.4, 18.9, 17.2, 15.0, 15.7, 17.9, 16.5]
payment_female = [14.2, 15.9, 16.0, 14.0, 17.0, 13.8, 15.2, 14.5, 15.0, 14.4]

等分散が仮定できる

2つの標本間で母分散が等しいと仮定できる場合、2つの標本の平均差 \bar{X}-\bar{Y} の分散を以下で計算して、2標本のt統計量を求めて仮説検定します。

$$ s^2=\frac{(n_X-1)s_X^2+(n_Y-1)s_Y^2}{n_X+n_Y-2} $$

$$ t=\frac{\bar{X}-\bar{Y}}{s\sqrt{\frac{1}{n_X}+\frac{1}{n_Y}}} $$

ナイーブに上の式でも計算できますが、ttest_ind()を使うと、2標本のデータを渡すだけでt統計量とp値を計算できます。

t統計量は3.61、p値は0.0020で有意水準(0.05/2=0.025)を大きく下回るため、帰無仮説 H_0 は棄却されます。

stats.ttest_ind(payment_male, payment_female)
# Ttest_indResult(statistic=3.606503775109937, pvalue=0.002017792820428875)

等分散を仮定できない

母分散が等しいと仮定できない場合、ウェルチの近似法でt統計量を近似します。

等分散でない平均の差は、以下の計算で得られた \upsilon に最も近い整数 \upsilon' を自由度とするt分布に従うことを利用して検定します(ウェルチのt検定と言います)。

$$ \upsilon=\frac{(\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2})^2}{\frac{(\frac{s_1^2}{n_1})^2}{n_1-1}+\frac{(\frac{s_2^2}{n_2})^2}{n_2-1}} $$

ttest_ind()のオプション引数equal_varをFalseとすることで、ウェルチのt検定が行われます(デフォルトではTrueで、これは上の通り等分散を仮定したt検定となります)。
賃金のデータに適用すると、t統計量は変わらず3.61でした。実際に \upsilon を計算してみると17.17で、 \upsilon' は17となり、上で求めた自由度18のt統計量と異なるはずですが、そうなってないです(四捨五入ではなく切り上げで補正されているのでしょうか)。
一方で、t統計量は同じにもかかわらず、p値が少し異なっているのも気になります(何か誤解しているかもしれませんので、注意してください)。

stats.ttest_ind(payment_1, payment_2, equal_var=False)
# Ttest_indResult(statistic=3.606503775109937, pvalue=0.0021488395200926777)

標本間で対応がある場合の平均差の検定

同じ個体で比較する場合、t統計量を以下で求めます。例えば、高血圧の人にある薬を投薬してその前後の血圧の平均の差を検定する、などです。このt統計量を使った検定を、対応ありのt検定と言います。
d_i=X_i-Y_i 、 \bar{d}=\frac{1}{n}\sum_i^n d_i とすると以下で計算されます。

$$ s^2=\frac{1}{n-1}\sum_i^n d_i^2 $$

$$ t=\frac{\bar{d}-\mu_0}{\sqrt{\frac{s^2}{n}}} $$

ここでは降圧剤を想定した架空の例として、投薬前の血圧値(before)と投薬後の血圧値(after)をデータに対して検定してみます。

SciPyでは、ttest_ind()の代わりに、ttest_rel()を使います。
t統計量1.56、p値0.14が得られました。

before = [141, 139, 140, 148, 138, 133, 142, 134, 139, 132, 140, 138]
after = [143, 134, 136, 140, 141, 133, 145, 128, 137, 133, 140, 134]

stats.ttest_rel(before, after)
# Ttest_relResult(statistic=1.560009076442849, pvalue=0.147047493905328)

母分散の仮説検定

母分散 \sigma^2 の仮説検定は、 \chi^2 検定統計量が自由度(n-1)の \chi^2 分布に従うことを利用して行います。

帰無仮説 H_0:\sigma^2=\sigma_0^2 としたときの \chi^2 は以下で計算されます。

$$ \chi^2=\frac{(n-1)s^2}{\sigma_0^2} $$

有意水準95%(両側検定)の場合、この \chi^2 統計量が \chi^2 分布の[ \chi^2_{0.975}(n-1), \chi^2_{0.025}(n-1) ]に含まれるかどうかは計算し、含まれない場合は帰無仮説を棄却します。
例えば、上の例の男性の賃金において、帰無仮説を \sigma^2_0=6 とした時の \chi^2 統計量は2.49と計算されます。自由度9における確率は0.019のため、帰無仮説は棄却されます。

sigma2_0 = 6.0

chi2_value = (len(payment_male) - 1) * pd.Series(payment_male).var() / sigma2_0
print(chi2_value) # 2.485999999999999

stats.chi2.cdf(chi2_value, len(payment_male) - 1)
# 0.018742372504659505

この \chi^2 検定は、適合度の検定や独立性の検定などに発展されます。

適合度の検定

想定していた確率分布から得られる度数(期待度数)と、実際に得られた度数(観測度数)が、一致していると言えるのかどうかを検定するのが適合度の検定です。
例えば、サイコロを50回振って出た目をカウントし、そのサイコロが歪んでないか(均等に各目が1/6の確率で現れているのか)を確かめる、といった問いに答える方法です。

k個のカテゴリ(サイコロであればk=6)について、理論確率 p_i から得られる各カテゴリの期待度数を np_i 、観測度数を f_i とすると、以下の式で \chi^2 統計量が計算できます。この統計量は自由度k-1の \chi^2 分布に従います。

$$ \chi^2=\sum_i=1^k \frac{(f_i-np_i)^2}{np_i} $$

サイコロを50回振った時の観測度数が以下の表で得られたとします。公平なサイコロであれば、期待度数は全て50/6です。このときに理論確率(全ての目が1/6の確率で出る)に対する帰無仮説 [\tex: H_0: p_1=1/6, ... , p_6=1/6 ] を立て、有意水準95%で検定します。

サイコロの目 1 2 3 4 5 6
観測度数 10 7 8 11 6 8
期待度数 50/6 50/6 50/6 50/6 50/6 50/6

SciPyでは、chisquare()の引数f_obsに観測度数、f_expに期待度数を与えることで、 \chi^2 統計値とp値を計算できます。
p値は0.84( \chi^2 は2.08)のため、帰無仮説は棄却できません。

stats.chisquare([10, 7, 8, 11, 6, 8], 
                f_exp=[50/6, 50/6, 50/6, 50/6, 50/6, 50/6])
# Power_divergenceResult(statistic=2.08, pvalue=0.8379684711477735)

独立性の検定

代数学と解析学の成績の関連など、2つの属性間の関連の有無についても [\tex: \chi^2 ] 検定で確認できます。独立性の検定と呼ばれます。

k個のカテゴリからなる属性Aと、m個のカテゴリからなる属性Bについて、観測度数を列挙します。

A_1 ... A_k 計
B_1 f_{1 1} ... f_{1 k} f_{1 \cdot}
... ... ... ... ...
B_m f_{m 1} ... f_{m k} f_{m \cdot}
計 f_{\cdot 1} ... f_{\cdot k} n

帰無仮説 H_0: \mbox{全てのi, jについて} P(A_i \cap B_j)=P(A_i)P(B_j) を立て、以下で \chi^2 統計量を計算します。この値は自由度(k-1)(m-1)の \chi^2 分布に従います。

$$ \chi^2=\sum_i \sum_j \frac{(nf_{ij}-f_{i \cdot}f_{\cdot j})^2}{nf_{i \cdot}f_{\cdot j}} $$

ここで、喫煙習慣の有無と6年後の生存/死亡を調べた表を用いて、2つが独立であることを帰無仮説として、有意水準95%で検定します(赤本の練習問題12.6からの抜粋)。

非喫煙者(non-smoker) 喫煙者(smoker) 計
生存(alive) 950 348 1298
死亡(dead) 117 54 171
計 1067 402 1469

SciPyではchi2_contingency()に表で渡すことで計算できます。
返り値は4つの要素を持つタプルとなっており、( \chi^2 統計量, p値, 自由度, 完全に独立の場合の期待度数)が格納されています。p値が0.22のため、帰無仮説が棄却できず、2つの属性が独立していないと判定できません。

smoking_df = pd.DataFrame({
    'non-smoker': [950, 117],
    'smoker': [348, 54]
}, index=['alive', 'dead'])

stats.chi2_contingency(smoking_df)
# (1.496886561580958,
#  0.22115103741056918,
#  1,
#  array([[942.79509871, 355.20490129],
#         [124.20490129,  46.79509871]]))

まとめ

上で紹介した以外にも、母集団の比率の検定や二項検定など、様々な検定がありますので随時触れていきたいと思います。