【ベイズの定理】事後分布から推定量を導く方法について|python
今回はベイズ統計学の基礎の基礎を解説いたします。
ベイズの定理の理解の後、ベイズ更新のやり方を学んでいく形にします。
ベイズ統計学
ここでは、ベイズ統計学の考え方をざっくりと学びましょう。
ベイズの定理に基づくベイズ統計は、経験を数値化して数学的に扱う手法だと言えます。
現在の事象をもとに、過去の考えを改めたり、過去の考えの正しさを強めたり、「経験的」な性質が強いです。この考え方を「ベイズ更新」を呼びます。
この考え方自体は、推測統計学よりも人間の考え方に近いと思います。
例えば、道の角から人が飛び出してくる事前確率は、0.5ですが、とある朝の通勤時間に道の角から人が飛び出してぶつかってしまいました。 そのため、朝は道の角から人が飛び出してくる確率がどうやら高いようなので、明日はもっと気をつけよう。
このような経験をもとに考えを更新していくことは、日常ではたくさんあると思います。
これを数学的に扱うのがベイズ統計学です。
ベイズの定理
では、ここでベイズの定理について解説します。
その前に、事前分布と事後分布について理解をする必要があります。
事前分布：データを観測する前のパラメータの分布
事後分布：データを観測した後のパラメータの分布
尤度に関しては、以下をご覧ください。
では、例を見てみましょう。
例：「トトロいたもん！」 メイちゃんは、どうやらトトロをみたらしいです。 しかし、森の中を探しても全然見つかりませんでした。 この時、メイちゃんが嘘つきである確率をベイズの定理を使って求めてください。 情報： メイちゃんが嘘つきの場合、「実際にトトロが居てもいなくても」、「トトロいたもん！」と出力します。 メイちゃんは普段嘘をつく方ではないので、メイちゃんが嘘つきである事前確率は0.1です。 メイちゃんが正直者の時に、トトロが見つかってしまう確率を0.8とします。 一方、メイちゃんが嘘つきの時にトトロが見つかってしまう確率は0.3とします。
今回求めるのは、以下の条件付き確率です。
$$P(嘘つき｜観測できない)$$
今回の確率は、ベイズの定理を使って以下のように求めることができます。
分母：\(P(観測できない) = P(観測できない | 嘘つき)*P(嘘つき) + P(観測できない | 正直者)*P(正直者) = 0.1 * 0.7 + 0.9 * 0.2 = 0.25\)
分子：\(P(嘘つき)*P(観測できない | 嘘つき) = 0.1 * 0.7 = 0.07 \)
よって、求める確率は\(0.28\)です。
ベイズの定理は一般的には、以下のような式です。
求めたい条件付き確率は、逆の条件付き確率を使えば求められる、ことになります。
$$P(A|B)={\frac{P(A)P(B|A)}{P(B)}}$$
「逆」というのは、AとBを原因と結果に分けた時にわかりやすくなると思います。
以下のように、「原因がわかれば結果がわかる」確率を、「結果がわかれば原因がわかる」確率に加工できます。
今回の例でいえば、「トトロが居なかった」という結果から、
1「メイちゃんは正直者だが、トトロが逃げた」
2「メイちゃんが嘘つきだった」
という原因の確率を求める、ということになります。
$$P(原因|結果)={\frac{尤度*P(結果|原因)}{P(結果)}}$$
ただ、P(結果)は往々にして計算が大変な場合があります。
今回のように単純な例を除き、事後分布の計算には「マルコフ連鎖モンテカルロ法(MCMC法)」という、「尤度と事前分布の積だけで事後分布を近似させる」という手法がよく使われます。
【完全ガイド】MCMC法についてわかりやすく解説|ベイズ推定
ベイズの定理における、尤度と事前分布と事後分布の関係は以下の通りです。
尤度という「メイちゃんが嘘つきの時に、観測できかった確率」を使って、事前分布(確率)という「嘘つきの確率」が
事後分布(確率)という「観測できなかったときに、嘘つきの確率」に変わりました。
$$事後分布= \frac{事前分布*尤度}{P(結果)}$$
さて、トトロを観測できなかった時にメイちゃんが嘘つきの確率が0.1→0.28とわかりました。
‥しっくりこないでしょうか。
それは、「トトロが逃げる確率」や「メイちゃんが嘘つきの事前確率」を決めうちで設定したから、それを元にした確率が微妙に感じられるのは最もでしょう。
しかし、ベイズ統計学が試行を繰り返していくことで真価を発揮します。
次のベイズ更新で見ていきましょう。
ベイズ更新
では、もう一度メイちゃんが「トトロいたもん！」と言ってきましたが、どれだけ探してもトトロを観測できませんでした。 この時の、メイちゃんが嘘つきの確率を求めてみましょう。 ただし、先ほどの事後確率の情報を使います。 メイちゃんが嘘つきの確率は0.1から0.28にベイズ更新されています。
この時求める確率は以下のようになります。
$$P(観測されない|嘘つき)= \frac{0.28*0.7}{0.28*0.7+0.72*0.2}$$
これを計算すると、0.576になります。
メイちゃんが嘘つきである確率が、0.1→0.28→0.58となりました。
何回も「トトロいたもん！」と言われても、実際にトトロを観測できないときは、メイちゃんが嘘つきである可能性がどんどん高まることになります。
このようにベイズ統計学は、試行結果を繰り返し使って確率を確からしいものに近づける考え方になります。
より一般化すると
例えば、観測した結果を \(D\) とします。
観測結果 \(D\) によって、確率変数 \(X\) に関する主観確率が変化します。
観測結果 \(D\) のもとで確率変数 \(X\) に関する主観確率を \(Pr(X|D)\) と書き、事後確率 (posterior) と呼んでいます。
この事後確率分布を求めるのがベイズ推量の目的です。
確率変数 X(以後パラメータと呼びます) の元では、観測結果 \(D\) が確率 \(Pr(D|X)\) に従って生成されるはずであることが分かっていたとします。
\(Pr(D|X)\) を \(D\) が観測された時の \(X\) の尤度 (likelihood) と呼びます。
さて、これまで紹介した事前確率 \(Pr(X)\)、事後確率 \(Pr(X|D)\)、尤度 \(Pr(D|X)\) の間には以下の関係が成り立ちます。
$$P(X|D)=\frac{Pr(D|X)Pr(X)}{\int Pr(D|X)Pr(X)}$$
CODE
ここでは簡単なベイズ更新のコードを紹介します。
まず事前分布として一様分布が仮定されています。
import pymc3 as pm
# 事前分布
with pm.Model() as model:
p = pm.Uniform('p', 0, 1)
# 観測データ
data = [1, 0, 0, 1, 1, 0, 1, 0, 1, 1]
# 事後分布
with model:
likelihood = pm.Bernoulli('likelihood', p, observed=data)
trace = pm.sample()dataというリスト型の変数に観測データを入れています。
観測データが0か1のベルヌーイ試行に従うものとしています。
また、確率pに対して一様分布を事前分布としています。事後分布は、pm.sample()を使ってサンプリングすることで求めることができます。
発展：共役事前分布について
先ほどのベイズの定理で事後分布を更新していくプロセスですが、事後分布が必ずしも解析的に扱いやすい確率分布ではないことがあります。
ただ、ある事前分布に対して、尤度関数の形を保ちつつ事後分布を計算しやすくする事前分布があり、それを自然共役事前分布と呼びます。
例えば、以下のようなベルヌーイ分布の尤度関数に対して
$$p(x|θ)=θ^x(1-\theta)^{1-x}$$
以下のようなβ分布が自然共役事前分布となります。
$$p(\theta)=θ^{\alpha-1}(1-\theta)^{\beta-1}$$
発展：MAP推定（最大事後確率推定）について
これは例えばデータが少ない時のように、最尤法だとパラメータが信頼できません。
→そのデータがまれにしか得られないデータかもしれないからですね。
そこで、与えられたデータに基づいて事後分布を最大化するパラメータの値を見つけるアプローチをMAP推定と呼びます。
$$\hat{\theta}_{MAP}=argmax_{\theta}p(\theta|D)=argmax_{\theta}p(D|\theta)p(\theta)$$
上の式で見るとわかりやすいですね。
最尤法だと、あるパラメータが与えられた時にデータが得られる確率（これを尤度と呼びます）\(p(D|\theta)\)が、最大化されます。
補足：頻度論統計学との違い
一般的な「統計学」と呼ばれる学問は、「頻度論的な考え」に基づいています。
ここでは、頻度論統計学とベイズ統計学の違いをざっくり解説します。
頻度論：関心のあるパラメータは未知の定数として扱う
ベイズ統計学：関心のあるパラメータは確率変数として扱う
ベイズ統計の大きな特徴は、「正真正銘の事前分布(事前確率)は、主観で決まる」ということです。
この点は、俗にいう頻度論者に嫌われる点です。試行を重ねるごとに、事前分布は客観的な情報を反映して信頼性を増してはいくものの、最初の事前確率は「決め打ち」なのです。
