8 マルコフ連鎖モンテカルロ（MCMC）法とベイズ統計モデル

8.1 例題：種子の生存確率（個体差なし）

例題：植物個体20個それぞれの種子8個を調べた時の生存種子数のデータをもとに，ある植物の種子の生存確率を推定する．

生存種子数$y_i$が二項分布にしたがうと仮定すると，ある個体$i$の生存種子数が$y_i$である確率は

p(y_i | q) 
= \begin{pmatrix}
8 \\ y_i
\end{pmatrix}
q^{y_i} (1-q)^{8-y_i}

尤度$L(q)$は，20個体分の「データが得られる確率」の積

L(q) = p(Y|q) = \prod_i p(y_i|q)

対数尤度は

\log L(q) = \sum_i \left[ y_i \log q + (8-y_i) \log (1-q) \right] + ({\rm 定数})

これを最大にする$q$が最尤推定値$\hat q$．今の場合，$d \log L(q)/dq = 0$から解析的に計算できる．

8.2 ふらふら試行錯誤による最尤推定

最尤推定量が解析的に求められない場合．
ここでは「効率が悪く精度も良くない試行錯誤による方法」を紹介．
生存確率を0.01きざみでとる．適当な初期値を考える．$q$=0.3．そこから$\pm$0.01動かして，対数尤度が高くなる方へ移動させる，という方法．
図8.5が結果．試行回数を重ねるにつれて$q$=0.46へ収束．$q$=0.6から始めた場合も．
これだと尤度関数がシンプルでない場合に局所的な解に落ち込んでしまう恐れありか．

8.3 MCMCアルゴリズムのひとつ：メトロポリス法

メトロポリス法．
尤度が大きければ更新，尤度が小さい場合は尤度の比を更新確率として更新するかどうかを決める．
このようなルールにもとづいて値を変化させる方法：MCMC
前の状態$q$にもとづいて新しい状態$q^{\rm new}$を作っている -> マルコフ連鎖
乱数を利用した計算アルゴリズム -> モンテカルロ法

8.3.1 メトロポリス法でサンプリングしてみる

図8.8に一例を示す．
MCMCの目的：特定の値の探索ではなく，ステップ数とともに変化するパラメータの値の生成（サンプリング）．
ステップ数が増えると，MCMCで得られたパラメータの頻度分布（マルコフ連鎖の「定常分布」）が確率密度関数に似たものになる．

8.3.2 マルコフ連鎖の定常分布

定常分布：ある変数$q$のマルコフ連鎖が一定の条件を満たしているときに，$q$がしたがう確率分布．
この章の例題のメトロポリス法によって作り出されたマルコフ連鎖はその条件を満たす（天下り的ですが）．
定常分布$p(q|Y)$．
ステップ数が小さい時は定常分布から大きく離れている．ステップ数が大きくなると定常分布に近づく．
図8.9．どの初期値から始めてもやがて定常分布に近づく．
効率の良いMCMCサンプリングの方法はいろいろあるらしい．
また，ステップ数の小さいときの結果を使わない，複数のMCMCサンプリングを比較する，といった方法も定常分布の推定の改善につながりそう（詳細は第9章以降）．

8.3.3 この定常分布は何をあらわす分布なのか？

この例題では，定常分布$p(q|Y)$は尤度$L(q)$に比例する確率分布．

p(q|Y) 
= \frac{L(q)}{\displaystyle \sum_q L(q)}

分母は規格化定数なので$p(q|Y) \propto L(q)$．
メトロポリス法によって得られる，十分に長いMCMCサンプリングは定常分布$p(q|Y)$からのランダムサンプリング．あるデータ$Y$に統計モデルをあてはめたときに$q$がとる値の確率分布と解釈することもできる．

8.4 MCMCサンプリングとベイズ統計モデル

種子の生存確率$q$を，ベイズの公式の形式で書くと

p(q|Y) 
= \frac{p(Y|q) p(q)}{\displaystyle \sum_q p(Y|q) p(q)}

左辺$p(q|Y)$は，データ$Y$が得られたときに$q$がしたがう確率分布．ベイズ統計学では事後分布と呼ぶ．
右辺$p(Y|q)$は，$q$の値が決まっているときにデータ$Y$が観測される確率．この例題の場合は二項分布の積である尤度$L(q)$で置き換えることができる．
$p(q)$は，データ$Y$がないときの$q$の確率分布．ベイズ統計モデルでは事前分布と呼ぶ．
右辺の分母は，$\sum_q p(Y|q) p(q) = p(Y)$で，$Y$というデータが得られる確率．
つまり，ベイズ統計モデルとは

{\rm 事後分布} 
= \frac{{\rm 尤度} \times {\rm 事前分布}}{{\rm データが得られる確率}}
\propto {\rm 尤度} \times {\rm 事前分布}

といった構造をもつ統計モデル．
8.3.3節の式と比較すると，事前分布$p(q)$が$q$の値によらずに定数であると無矛盾．図で表すと図8.13のようになる．

8.5 補足説明

8.5.1 メトロポリス法と定常分布の関係

メトロポリス法にしたがって，パラメータ$q$のサンプリングをすると，なぜ定常分布$p(q|Y)$からのランダムサンプルになるのか．

そのための条件

$q$が任意の初期値から定常分布$p(q|Y)$に収束する．
ある$q$が$p(q|Y)$にしたがっていて，メトロポリス法で$q^{\rm new}$を得た時に，この$q^{\rm new}$も$p(q|Y)$にしたがっている．

ここでは簡単のため2.の条件が成り立っていることだけを確認．

ちょっとした計算から詳細つり合いの条件が導かれる．

p(q^{\rm new} | Y) p(q^{\rm new} \to q) 
= p(q | Y) p(q \to q^{\rm new})

両辺$q$について和を取ると，左辺の最初の因子は$q$に無関係で，$\sum_q p(q^{\rm new} \to q) = 1$なので，

p(q^{\rm new} | Y) 
= \displaystyle \sum_q p(q | Y) p(q \to q^{\rm new})

これより，$q$が定常分布$p(q|Y)$にしたがっているなら，メトロポリス法によって選ばれる$q^{\rm new}$もまた定常分布$p(q^{\rm new}|Y)$にしたがっている¹．これはメトロポリス法に限らず，すべてのMCMCアルゴリズムに共通する性質．

8.5.1 ベイズの定理

8.4節の式はベイズの定理と呼ばれる．

p(Y,q) = p(q|Y) p(Y) \\
p(Y,q) = p(Y|q) p(Y)

からこれらの式の右辺が等しいとして

p(q|Y) = \frac{P(Y|q) p(q)}{p(Y)}

ここで$p(Y) = \sum_q p(Y|q) p(q)$と書けることに注意すると8.4節の式になる．

8.6 この章のまとめと参考文献

これまでの章では非ベイズの統計モデル．データのあてはまりの良さを尤度で評価し，その尤度を最大化するようにパラメータを選んだ．
ベイズ統計モデルでは，尤度だけでなくパラメータの事前分布も組み込まれ，事後分布があてはめの結果として得られる．

と書いてあるけどこの記述だけではいまいちよくわからない． ↩