ベイズ統計の仮説検定
この記事ではベイズ統計で仮説検定をする手法を紹介し、最後に実際にpythonを用いて検定をしてみます。筆者はベイズ統計もpythonも初心者なので、内容に不備があるときはコメントをお願いします。
1.仮説検定とは
$ $大雑把に言うと、仮説検定とは未知のパラメータ$\theta$がある区間$I_i$に属するという仮説$H_i$を立てて$H_i$という仮説を採択してよいか、つまり未知のパラメータ$\theta$が区間$I_i$に属するといってよいのかどうか分析する手法です。
ベイズ統計では主に以下の3つの手法が代表的なものとなっております。
- 事後確率による検定
- 事後オッズ比による検定
- ベイズ因子による検定
$ $以下では確率変数$X$の確立密度関数を$p(X)$, 事象$X$の起こる確率を$P(X)$で表します。また、$\Omega$を標本空間とし、集合$A, B$に対して$A \cap B = \phi$が成り立つことを$A \bot B$で表します。
2.事後確率による検定
事後確率の評価による検定の考え方は非常に単純です。例えば以下のような仮説を考えてみましょう。
\left\{
\begin{split}
H_0: &q \in I_0 \\
H_1: &q \in I_0 \\
&\vdots \\
H_k: &q \in I_k
\end{split}
\right.
$ $( ただし、$I_i \subset \Omega$であり、$I_i$同士は領域を共有してもかまわないことに注意してください。)
$ $このとき、事象$X$が実現したとき仮説$H_i ( i = 0, 1, \dots k )$が成り立つ事後確率はそれぞれ
P(H_i|X) = P(\theta \in I_i|X) = \int_{I_i} p(\theta|X) d\theta
となります。
$P(H_i|X)$が十分に大きい(1に近い)ならば仮説$H_i$を採択し、小さい(0に近い)ならば仮説$H_i$を棄却するといった要領で検定できます。
3.事後オッズ比による検定
次に事後オッズ比を用いた検定です。事後確率による検定では好きなだけ仮説を立てることができましたが、事後オッズ比を用いた検定で用いる仮説は 2 つだけです。
以下のような仮説を考えます。
\left\{
\begin{split}
H_0: q \in I_0 \\
H_1: q \in I_1
\end{split}
\right.
$ $ここで、$I_0$ と$I_1$は$I_0 \bot I_1 , I_0 \cup I_1 = \Omega$が成り立つように設定します。このようにすることで、$H_0$ と$H_1$のうち一方を帰無仮説、他方を対立仮説とすることができます。
事後オッズ比は、
\begin{split}
事後オッズ比 &= \frac{P(H_0|X)}{P(H_1|X)} \\
&= \frac{P(\theta \in I_0|X)}{P(\theta \in I_1|X)} \\
&= \frac{\int_{I_0}p(\theta|X)d\theta}{\int_{I_1}p(\theta|X)d\theta}
\end{split}
$ $で定義されます。事後オッズ比が1より十分大きければ$H_0$を、十分小さければ$H_1$を採択します。
$ $以上が事後オッズ比を用いた検定の手法ですが、事後分布が事前分布の影響を受けることを考慮すると少し問題があります。事前分布がどちらか一方の仮説にとって極端に有利であった場合、事後分布が事前分布の影響を反映して正しく検定できないかもしれません。例えば、本当は$H_1$が正しいのに$H_0$に圧倒的に有利な事前分布を仮定したとすると、$P(H_0|X)$が実際より大きく$P(H_1|X)$が実際より小さく見積もられ、事後オッズ比がかなり大きく見積もられる可能性があります。こういったことを避けるためにベイズファクターを用いて検定することがあります。
4.ベイズファクターによる検定
$ $早速ですがベイズファクター$B_{01}$は以下の式で定義されます。
\begin{split} B_{01} &= 事後オッズ比 \div 事前オッズ比 \\
&= \frac{P(H_0|X)}{P(H_1|X)} \div \frac{P(H_0)}{P(H_1)} \\
&=\frac{\int_{I_0}p(\theta|X)d\theta}{\int_{I_1}p(\theta|X)d\theta} \div \frac{\int_{I_0}p(\theta)d\theta}{\int_{I_1}p(\theta)d\theta}
\end{split}
ではなぜこのベイズファクターによる評価が事後オッズ比による検定での問題を回避することになるのか考えてみましょう。
$ $事前分布の設定が圧倒的に$H_0$に有利で、事前オッズ比が100、事後オッズ比が20となったとします。このとき、事後オッズ比による検定では$H_0$が採択されることになります。また、念のため事後確立による検定でどうなるかも考えてみると$P(H_0|X) \simeq 0.95$となりこちらも$H_0$が採択されるといってもよいでしょう。ですがベイズファクター$B_{01}$を計算すると、$B_{01} = 0.2$となります。ベイズファクターの値は以下に示したJeffreysによる評価基準を用いて評価します。$1/10 \le 0.2 \le 1/3$ですから、この場合対立仮説である$H_1$を十分に支持することになります。
[出典:https://www.slideshare.net/masarutokuoka/mcmc-35640699]
以上の例は極端ですが、実際にベイズファクターを用いた検定をすることで事後オッズ比を用いた場合とは異なる結論がえられることがあり、ベイズファクターを用いた検定のほうが事前分布の影響を小さくしていてより誤った結論を導きにくいと考えられます。
5.実際にpythonを使ってサンプルデータを検定
実際に事後確率、事後オッズ比、ベイズファクターを計算して仮説検定をやってみました。今回は簡単のため事前分布はベータ分布を仮定し、ベルヌーイ分布に従うようなモデルに対して仮説検定を行ってみました。
まずサンプルデータ (data) 、事前分布のハイパーパラメータ (a0, b0) 、区間 (I) を引数に取り、事後確率 (post_prob)、事後オッズ比 (post_odds_ratio) 、ベイズファクター (BF) を返す関数を定義しました。
import scipy.stats as st
def ppb(data, a0, b0, I):
n = data.size
y = data.sum()
pre_p0 = st.beta.cdf(I[1], a0, b0) - st.beta.cdf(I[0], a0, b0)
pre_p1 = 1 - pre_p0
post_p0 = st.beta.cdf(I[1], y + a0, n - y + b0) - st.beta.cdf(I[0], y + a0, n - y + b0)
post_p1 = 1 - post_p0
post_prob = post_p0
post_odds_ratio = post_p0 / post_p1
BF = post_odds_ratio * pre_p1 / pre_p0
return post_prob, post_odds_ratio, BF
続いてパラメータ p = 0.7 のベルヌーイ分布からランダムにデータを生成し、サンプルサイズ100の標本を生成しました。
import numpy as np
p = 0.7
n = 100
np.random.seed(0)
data = st.bernoulli.rvs(q, size=n)
このデータに対して、ハイパーパラメータが
- ア) a0 = b0 = 1
- イ) a0 = 2 , b0 = 5
- ウ) a0 = 3 , b0 = 2
の場合について分析を行いました。帰無仮説として用いる区間(I)は I = [0.6, 0.8]としました。
a0 = [1,2,3]
b0 = [1,5,2]
I = [0.6,0.8]
s = ['ア','イ','ウ']
for i in range(len(a0)):
a, b, c = ppb(data, a0[i], b0[i], I)
print('{}) \n 事後確率: {:.5f} \n 事後オッズ比: {:.5f} \n ベイズファクター: {:.5f}'.format(s[i],a,b,c))
実行結果
ア)
事後確率: 0.79632
事後オッズ比: 3.90973
ベイズファクター: 15.63891
イ)
事後確率: 0.93228
事後オッズ比: 13.76698
ベイズファクター: 336.00387
ウ)
事後確率: 0.81888
事後オッズ比: 4.52127
ベイズファクター: 8.62195
ここで、先ほど示したJeffreysによるベイズファクターの評価基準を参考に検定するとア)帰無仮説を強く支持、イ)帰無仮説をとても強く支持、ウ)帰無仮説を十分に支持という結論が得られました。pの真の値を0.7に設定しているので、いずれも妥当な結論を導いたといえます。
また、以下のベータ分布のグラフを見ると、イ)が最も事前分布の影響を受けそうです。
実際に事後分布のグラフを描画したものが以下の画像です。確かにイ)はア)、ウ)に比べて事前分布の影響を大きく受けているように見えます。
さらに、帰無仮説の区間だけを[0.65, 0.75]に変えてみると、以下のような結果になりました。
ア)
事後確率: 0.34440
事後オッズ比: 0.52532
ベイズファクター: 4.72784
イ)
事後確率: 0.57232
事後オッズ比: 1.33818
ベイズファクター: 74.33720
ウ)
事後確率: 0.36755
事後オッズ比: 0.58115
ベイズファクター: 2.73401
条件が厳しくなり事後確率や事後オッズ比だけでは判断しにくいですが、ベイズファクターを用いればいずれの場合も帰無仮説を支持できなくもなさそうです。
最後に帰無仮説の区間を[0.4, 0.6]として真の値が区間に含まれない場合どのような結論が得られるのか実験してみました。
ア)
事後確率: 0.00019
事後オッズ比: 0.00019
ベイズファクター: 0.00078
イ)
事後確率: 0.00123
事後オッズ比: 0.00123
ベイズファクター: 0.00515
ウ)
事後確率: 0.00020
事後オッズ比: 0.00020
ベイズファクター: 0.00049
77
ベイズファクターを用いるまでもないですが、いずれの場合もきちんと帰無仮説は棄却されます。
6.最後に
事前分布をベータ分布で仮定し、ベルヌーイ分布に従うようなモデルしか扱いませんでしたが、今後もっと様々なモデルも扱っていきたいと思っています。
何か不備、誤り等ございましたらご教授いただければと思います。
最後まで読んでいただきありがとうございました。