ベイズ推定のキホン

やること

この投稿では「ベイズの定理」導出から、それを「ベイズ推定」として応用するまでの考え方を追っていきます。理解の助けに、簡単な例題と Python の基本ライブラリのみを使ったコードスニペットを示します。ベイズ統計モデリングのベースとなる理論的な展開ポイントを説明しますが、MCMC を使ったモデリングの実践的な話には至りませんのでご了承ください。

「ベイズの定理」の導出

「ベイズの定理」は「条件付確率」から次のように導出します。

条件付確率(ある事象 \( A \) が起こったという条件の下で別の事象 \( B \) が起こる確率):

\[ P(B|A) = \frac{P(A \cap B)}{P(A)} \tag{1} \]

両辺に \( P(A) \) を掛けて左辺と右辺を入れ替えると、

乗法定理:
\[ P(A \cap B) = P(B|A)P(A) \tag{2} \]

が得られます。ここで右辺について、\( A \) の役割を \( B \) に担わせれば(単に記号だからどっちでもいいので)
\( P(A \cap B) = P(A|B)P(B) \) に変換でき、上記 乗法定理 の右辺とも等値なので、
\( P(A|B)P(B) = P(B|A)P(A) \) が得られ、これを \( P(A|B) \) について解くと、

ベイズの定理(基本版):
\[ P(A|B) = \frac{P(B|A)P(A)}{P(B)} \tag{3} \]

が導出されます。尚、ベイズ統計では、ある事象 \( A \) が起こる確率 \( P(A)\) のことを「事前確率」、
その条件の下で別の事象 \( B \) が起こる確率 \( P(A|B) \) のことを「事後確率」といいます。

※式を変換しただけで、「ベイズの定理 は 条件付確率 そのものである」と理解するのがポイント。

「ベイズの定理」を実用的な式に変換する

上記で得られたベイズの定理(基本版)を発展させて、より実用的な形に変形させていきます。
複数事象の組み合わせの中で「ベイズの定理」を扱う準備として、これまで単体で扱っていた事象 \( A \) を、
事象 \( A_1 \)、\( A_2 \)、\( A_3 \) という3つの排反事象として分割させてみます。

このとき、ある事象が\( a \)個に分割されているとき、それぞれの確率を足し上げると1(= 100%)になる、つまり、

\[ \sum_{i=1}^{a}p(A_i) = 1 \tag{4} \]

となる点に注意してください。ベイズの定理(基本版)の \( A \) を \( A_1 \) に書き換えると単に

\[ P(A_1|B) = \frac{P(B|A_1)P(A_1)}{P(B)} \tag{5} \]

と書けます。ここで、事象 \( B \) は事象 \( A_1 \)、\( A_2 \)、\( A_3 \) いずれかに属するものとすると、\( P(B) \) は次のように3つの和で表現できます。

\[ P(B) = P(B \cap A_1) + P(B \cap A_2) + P(B \cap A_3) \tag{6} \]

上式の右辺は 乗法定理 により次のように変換できます。

\[ P(B) = P(B|A_1)P(A_1) + P(B|A_2)P(A_2) + P(B|A_3)P(A_3) \tag{7} \]

これを最初の \( (5) \) 式に代入すると、

\[ P(A_1|B) = \frac{P(B|A_1)P(A_1)}{P(B|A_1)P(A_1) + P(B|A_2)P(A_2) + P(B|A_3)P(A_3)} \tag{8} \]

と書け、さらに一般化すると、

ベイズの定理(応用版):
\[ P(A_i|B) = \frac{P(B|A_i)P(A_i)}{\sum_{i=1}^{a}P(B|A_i)P(A_i)} \tag{9} \]

を得ます。上記では\( A_1 \)、\( A_2 \)、\( A_3 \) と事象を3つの排反事象に分割した例を示しましたが、数学の試験問題とかでは文章をきちんと読んで、排反事象をいくつに分割するか考えることが解答のポイントになります。あとは問題の文章からそれぞれの確率を拾って公式的に代入して解いてしまえばOK。

例題
赤玉と白玉が入った袋が2つ存在する。(それぞれA,Bとする)
袋Aには赤玉が2つ、白玉が1つはいっている。袋Bには赤玉が1つ、白玉が4つはいっている。
ある人がどちらかの袋から玉をひとつ取り出し、あなたにその玉の色だけを伝える。
そのとき選ばれた袋がAである確率を推定せよ。

解答
– 袋Aから赤玉を取り出す確率: \( P(r|A) = 2/3 \) ← 問題文より
– 袋Aから白玉を取り出す確率: \( P(w|A) = 1/3 \) ← 問題文より (※ \( 1 – P(r|A) \) )
– 袋Bから赤玉を取り出す確率: \( P(r|B) = 1/5 \) ← 問題文より
– 袋Bから白玉を取り出す確率: \( P(w|B) = 4/5 \) ← 問題文より (※ \( 1 – P(r|B) \) )
– 袋Aから赤玉を取り出す確率: \( P(A) = 1/2 \) ← 理由不十分の原則より
– 袋Bから赤玉を取り出す確率: \( P(B) = 1/2 \) ← 理由不十分の原則より (※ \( 1 – P(A) \) )

取り出した玉が赤玉だった場合、選ばれた袋がAである確率は
\[ \begin{align}
P(A|r) &= \cfrac{P(r|A)P(A)}{P(r|A)P(A) + P(r|B)P(B)} \\
&= \cfrac{2/3 \times 1/2}{(2/3 \times 1/2) + (1/5 \times 1/2)} \\
&= \cfrac{10}{13} \approx 76.92 \%
\end{align} \]

取り出した玉が白玉だった場合、選ばれた袋がAである確率は
\[ \begin{align}
P(A|w) &= \cfrac{P(w|A)P(A)}{P(w|A)P(A) + P(w|B)P(B)} \\
&= \cfrac{1/3 \times 1/2}{(1/3 \times 1/2) + (4/5 \times 1/2)} \\
&= \cfrac{5}{17} \approx 29.41 \%
\end{align} \]

統計モデルに「ベイズの定理」を適用する

「確率」に関する「ベイズの定理」を「分布」に関する「ベイズの定理」に拡張します。ここからは、5つのステップで拡張に至るロジックを確認していきます。

Step 1. ベイズ統計における「原因」と「結果」(逆確率を理解する)

通常の条件付確率で論じられる場合、因果関係はその時系列をたどって \( P(結果|原因) \) と、「原因から結果」の確率を求めようとします。一方で、ベイズの定理では、時間の流れに逆らって「結果から原因」の確率 \( P(原因|結果) \) を求めようとします。「ベイズの定理」で結果から原因を推量した確率であるとき「事後確率」のことを「逆確率」ということがあります。

逆確率: 「結果 → 原因」を推し量る。 \( P(原因|結果) \)

Step 2. 原因(の仮定 \( H \) )と結果(のデータ \( D \) )を考える

事象\( A \) を「原因(仮定: \( H \) ypothesis)」、事象 \( B \) を「その結果の事象(データ: \( D \) ata)」と読み替えると、ベイズの定理は次のように書き換えて解釈できます。

\[ P(H|D) = \cfrac{P(D|H)P(H)}{P(D)} \tag{10} \]

\[ P(H_i|D) = \frac{P(D|H_i)P(H_i)}{\sum_{i=1}^{a}P(D|H_i)P(H_i)} \tag{11} \]

このように変換することで、「データ \( D \) という結果が与えられたとき、その原因が、特定の仮定 \( H_1 \) である確率」を求める式として解釈しやすくなります。

Step 3. 原因(仮定)\( H \) を確率分布のパラメーター \( \theta \) に置き換える

統計学では、確率分布の 母数(パラメーター) を推定することが重要なテーマとなります。確率分布の性質を決定づけることが出来れば、そこから生成される値によって現象の説明や予測が出来るからです。「ベイズの定理」を統計学で役立てるには、仮定 \( H \) を「母数が \( \theta \) の値をとるときの確率分布」と読み替える必要があります。そこで、次のように書き換えます。

\[ P(\theta|D) = \cfrac{P(D|\theta)P(\theta)}{P(D)} \tag{12} \]

\[ P(\theta_i|D) = \frac{P(D|\theta_i)P(\theta_i)}{\sum_{i=1}^{a}P(D|\theta_i)P(\theta_i)} \tag{13} \]

文章としての仮定 \( H \) を、数値としての確率分布のパラメーター \( \theta \) と解釈し直して「ベイズの定理」を利用するのがポイント。

Step 4. パラメーター \( \theta \) を連続変数として解釈する

確率分布においてパラメーターは連続変数であることが一般的です。しかし、\( \theta \) が連続変数のときには、\( P(\theta|D) \) 、\( P(D|\theta) \) 、\( P(\theta) \) は 0 となってしまうため、これを「確率」と解釈できなくなります。
※ イメージ:\( \theta = 1.111111… \) である場合の \( P(\theta) \) は限りなく 0。無限に足し算をしなければいけない。

以後、基本的には連続的な確率変数を扱うので、これまで「確率」としていたものを「確率密度関数」または「確率分布」と読み替えます。離散変数と連続変数を区別するため、次のように表記を分けます。

– 事前確率 \( P(\theta) \) → 事前分布 \( \pi(\theta) \)
– 尤度 \( P(D|\theta) \) → 尤度 \( f(D|\theta) \)
– 事後確率 \( P(\theta|D) \) → 事後分布 \( \pi(\theta|D) \)

したがって、確率変数が連続の場合の「ベイズの定理」は次のように表記になります。

\[ \pi(\theta|D) = \cfrac{f(D|\theta)\pi(\theta)}{P(D)} \tag{14} \]

Step 5. 計算に必要な要素をまとめる(「カーネル」について)

\( (14) \) 式にある分母 \( P(D) \) は、「データ \( D \) の得られる確率」を意味します。ただ、推定したい \( \theta \) とは関係なく、実際に分析する場面ではデータ \( D \) が所与であるので、よりコンパクトに理解するために定数 \( k \) をつかって

\[ \pi(\theta|D) = kf(D|\theta)\pi(\theta) \tag{15} \]

とおきます。定数 \( k \) は確率の総和が1、すなわち \( \theta \) のすべてについて和(積分)が1になる、という性質に調整して決定されるものです。この定数 \( k \) のことを「正規化定数」と呼びます。この式全体の積分の値が 1 となるのは偶然ではなく、恣意的に正規化定数 \( k \) の値を選んで決定している点がポイントです。ちょっと乱暴な操作に感じるので、「そんなことしちゃっていいのかな」と思いがちですが、分布の性質(形状)を維持して、それが「確率密度関数(面積全部足したら100%)」として働いてもらうために必要な処理となります。
尚、確率分布において、母数 \( \theta \) を含んでいるの部分( \( (15) \) 式では右辺の \( f(D|\theta)\pi(\theta) \) )のことを「カーネル」といい、その確率分布の本質的な性質を決定付けています。調整役の正規化定数 \( k \) をいちいち式に明示するのではなく、比例定数として次のように記載することでカーネルを強調するような書き方をとることがよくあります。

事後分布 \( \propto \) 尤度 \( \times \) 事前分布
\[ \pi(\theta|D) \propto f(D|\theta)\pi(\theta) \tag{16} \]

※ \( \propto \) :「左辺は右辺に比例する」という意味の記号(左辺 is proportional to 右辺)「 \( \propto \) 」を見たら、『総和が1となるような「正規化定数 \( k \) 」が隠れているだけ』と理解するのがポイント。

以上で「分布」に関する「ベイズの定理」の拡張が完了しました。いろいろな操作がありましたが、各ステップを辿ることで、「ベイズ推定」のスタートラインに立つことが出来ます。

では、実際の「ベイズ推定」の利用イメージを膨らませるために、シンプルな例題を Python で実装しながら解いていきます。

例題
表が出る確率が\( \theta \) である1枚のコインがある。
このコインを5回投げた時、{表、表、裏、表、裏}の順に出た。
このとき、表の出る確率 \( \theta \) の事後分布 \( \pi(\theta|D) \) を求めよ。

解答
まずは \( \pi(\theta|D) \propto f(D|\theta)\pi(\theta) \) を思い出します。

– コインの表が出る確率: \( f(表|\theta) = \theta \)
– コインの裏が出る確率: \( f(裏|\theta) = 1 – \theta \)

初めに「理由不十分の原則」より、事前分布 \( \pi(\theta) \) を「一様分布」とします。
事前分布: \( \pi(\theta) = 1 \space (0 \leq \theta \leq 1) \)

実際にベイズ更新が起こるプロセスを matplotlib で示すと次のようになります。

    import numpy as np
    import matplotlib.pyplot as plt
    plt.style.use('ggplot')
    %config InlineBackend.figure_format = 'retina'

    # 事前分布として一様分布を設定
    def pi_0(theta):
        return 1.0

    def veiw_pi(pi):
        '''事後分布の形を描画'''
        fig = plt.figure(figsize = (10, 5))
        ax = plt.axes()

        theta = np.arange(0,1.1,.1)
        ax.plot(theta, [pi(t) for t in theta])
        ax.set_ylabel(' $\pi(\\theta)$ ')
        ax.set_xlabel(' $\\theta$ ')
        ax.grid(True);

    veiw_pi(pi_0)
[Out:]
        

この事前分布に対して、「1回目に表が出た」という観測データを取り込みます。
「1回目に表が出る」という事象を\( D_1 \) とすると、尤度は\( f(表|\theta) = \theta \) で示されるので、

1回目の事後分布(): \( \pi(\theta|D_1) \propto \)\( \theta \)\( \times 1 = \)\( \theta \)

※反映箇所を赤く示しています。

\( (0 \leq \theta \leq 1) \) の定義域で確率の総和が1という条件から、正規化定数は \( k = 2 \) と求められます。
また、1回目のデータを取り込んだという意味で、\( \pi_1(\theta) \) と表記します。すなわち、

\[ \pi_1(\theta) = 2\theta \]

となり、グラフは次のようになります。

    # ベイズ更新 1
    def pi_1(theta):
        return 2 * theta

    veiw_pi(pi_1)
[Out:]
        

さらに、「2回目に表が出る」という事象を \( D_2 \) とすると、1回目と同様に、尤度は \( f(表|\theta) = \theta \) で示されるので、2回目の事後分布は1回目の事前確率に掛け合わせることで求まります。

2回目の事後分布(): \( \pi(\theta|D_2) \propto \)\( \theta \)\( \times 2\theta = 2\)\(\theta^2 \)

\( (0 \leq \theta \leq 1) \) の定義域で確率の総和が1という条件から、正規化定数は \( k = 3 \) と求められ、

\[ \pi_2(\theta) = 3\theta^2 \]

となり、グラフは次のようになります。

    # ベイズ更新 2
    def pi_2(theta):
        return 3 * theta ** 2

    veiw_pi(pi_2)
        [Out:]
        

このように、事後分布となった結果を事前分布において、新たなデータを取り込んで次の事後分布を計算することを「ベイズ更新」と言います。さらに「3回目:裏」「4回目:表」「5回目:裏」というデータを取り込んでベイズ更新を進めます。それぞれの事象を \( D_3 \)、\( D_4 \)、\( D_5 \) とすると、これまでと同様に、

3回目の事後分布(): \( \pi(\theta|D_3) \propto \)\((1 – \theta) \)\( \times 3\theta^2 = 3 \)\((1 – \theta) \)\(\theta^2 \)
規格化の条件より、正規化定数は \( k = 12 \) なので、

\[ \pi_3(\theta|D_3) = 12(1-\theta)\theta^2 \]

4回目の事後分布(): \( \pi(\theta|D_4) \propto \)\(\theta \)\( \times 12(1-\theta)\theta^2 = 12(1-\theta)\)\(\theta^3 \)
規格化の条件より、正規化定数は \( k = 20 \) なので、

\[ \pi_4(\theta|D_4) = 20(1-\theta)\theta^3 \]

5回目の事後分布(): \( \pi(\theta|D_5) \propto \)\((1 – \theta) \)\( \times 20(1-\theta)\theta^3 = 20\)\((1-\theta)^2 \)\(\theta^3 \)
規格化の条件より、正規化定数は \( k = 60 \) なので、

\[ \underline{\pi_5(\theta|D_5) = 60(1-\theta)^2\theta^3} \]

上の式がこの例題の答えになります。ちなみに、正規化定数 \( k \) は代数計算ライブラリの SymPy を使い、積分計算で次のように求めることが出来ます。

    import sympy as sym
    theta = sym.Symbol('theta')
    sym.integrate((1 - theta) * theta ** 2, (theta, 0, 1))
    [Out:] 1/12
    sym.integrate((1 - theta) * theta ** 3, (theta, 0, 1))
    [Out:] 1/20
    sym.integrate(((1 - theta) ** 2) * (theta ** 3), (theta, 0, 1))
    [Out:] 1/60

また、それぞれのグラフは次のようになります。

    # ベイズ更新 3
    def pi_3(theta):
        return 12 * (1 - theta) * theta ** 2

    # ベイズ更新 4
    def pi_4(theta):
        return 20 * (1 - theta) * theta ** 3

    # ベイズ更新 5
    def pi_5(theta):
        return 60 * ((1 - theta) ** 2) * (theta ** 3)
    veiw_pi(pi_3)
        [Out:]
        
    veiw_pi(pi_4)
        [Out:]
        
    veiw_pi(pi_5)
        [Out:]
        

「ベイズ推定」(3つの点推定)

上記の例題で得られた事後分布について、グラフでの出力結果をみれば、推定値 \( \hat{\theta} \) は大体 \( 0.6 \) くらいだろうというのが見て取れます。解析的に点推定値を求める際、ベイズ推定では次の3つの方法が知られています。

1. EAP推定値(事後期待値、Expected A Posteriori)

事後分布 \( \pi(\theta | D ) \) による期待値 \( \theta_{eap} \) を求めるやり方。

\[ \begin{align}
\hat{\theta}_{eap} &= E[\theta|D] \\
&= \int \theta \pi(\theta|D)d\theta \\
&= \int \theta \frac{f(D|\theta)\pi(\theta)}{f(D)}d\theta \tag{17}
\end{align} \]

# 事後分布を定義
f = 60 * ((1 - theta) ** 2) * (theta ** 3)

# 期待値を計算する
EAP = sym.integrate(theta * f, (theta, 0, 1))
float(EAP)
    [Out:] 0.5714285714285714

2. MAP推定値(事後最大値、Maximum A Posteriori)

事後分布 \( \pi(\theta | D ) \) による確率密度が最大となる \( \hat\theta_{map} \) を求めるやり方。

\[ \hat\theta_{map} = arg\,max\,\pi(\theta|D) \tag{18} \]

一般的には、尤度関数について最尤推定法がとる解法と同様に

  1. 事後分布の対数を計算
  2. 母数で微分して 0 と置く
  3. 母数に関して方程式を解く

の順に求めます。

# 事後分布の対数を計算
log_f = sym.log(f)

# 母数で微分して 0 と置く
eq = sym.Eq(sym.diff(log_f), 0)
        
# 母数に関して方程式を解く
MAP = sym.solveset(eq).args[0]
float(MAP)
    [Out:] 0.6

※ 事前分布として一様分布を利用すると、MAPは、最尤推定の結果と同じになります。主観によらず、所与のデータから客観的に導き出される最尤推定量と一致するという意味で、一様分布は主観に頼らない事前分布として 無情報的な事前分布 であると言われます。

3. MED推定値(事後中央値、Posterior Median)

事後分布 \( \pi(\theta | D ) \) の積分がちょうど半分の 0.5(メディアン)になる \( \theta_{med} \) を求めるやり方。

\[ \int_{-\infty}^{\hat\theta_{med}}\pi(\theta|D)d\theta = 0.5 \tag{19} \]

# 分布関数が 0.5 になるとき
eq_2 = sym.Eq(sym.integrate(f), 0.5)

# 定義域を0から1に指定して解く
MED = sym.solveset(eq_2, theta, sym.Interval(0, 1)).args[0]
MED
    [Out:] 0.578592809309287

\( \hat\theta_{EAP} \)、\( \hat\theta_{MAP} \)、\( \hat\theta_{MED} \) のどれが優れているというわけではなく、それぞれ事後分布の解釈に使います。各値の違を視覚的に示すと次のようになります。

fig = plt.figure(figsize = (10, 5))
ax = plt.axes()

x = np.arange(0, 1, .01)
pi_t = [pi_5(t) for t in x]

ax.plot(x, pi_t, label = ' $\pi(\\theta)$ ')
ax.vlines(x = EAP, ymin = 0, ymax = 2.3, colors = 'green', linestyles = 'dotted', label = ' $\hat\\theta_{EAP}$ ')
ax.vlines(x = MAP, ymin = 0, ymax = 2.3, colors = 'blue', linestyles = 'dotted', label = ' $\hat\\theta_{MAP}$ ')
ax.vlines(x = MED, ymin = 0, ymax = 2.3, colors = 'orange', linestyles = 'dotted', label = ' $\hat\\theta_{MED}$ ')
ax.set_ylabel(' $\pi(\\theta)$ ')
ax.set_xlabel(' $\\theta$ ')
ax.legend(loc = 2, fontsize = 15)
ax.grid(True);
    [Out:]
    

参考図書

この投稿の内容について次の本を参考にしました。

道具としてのベイズ統計
涌井 良幸
日本実業出版社
売り上げランキング: 147,352

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です