言葉で理解するメトロポリスヘイスティング法
自己紹介
医学部に通う学生です、数学は頑張りましたが統計は"完全なる"素人から始めました。式をたくさん振り回すのではなく言葉を尽くしてまとめていきます、直感的に理解できるのが目標です。真面目な解説は他の人の記事の方が圧倒的に良いと思ます、非理系がどうにか概念を掴むための備忘録です。統計検定準一級取得、次はJDLA E資格に向けて勉強中です。
1.事後分布を理解する
まず、ベイズの基本の式$$f(\theta|x)=\frac{f(x|\theta)f(\theta)}{f(x)}$$
について、これはつまり元々$f(\theta)$という分布が自らの手元にあるときに、$x$というデータを手に入れたときに分布$f(\theta)$が「更新」されて$f(\theta|x)$となることを意味していました。
【具体例】
迷惑メールが手元に来る確率を$f(\theta)$としたときに、たくさんの観測結果から迷惑メールであった時に$x$という言葉が入っている確率$f(x|\theta)$を私たちは客観的に手に入れることが出来ます。そのときに、「$x$という言葉が入っているときにそれが迷惑メールである確率」(逆確率とも呼ばれます)は$f(\theta|x)$と表せて、これは$f(\theta|x)=\frac{f(x, \theta)}{f(x)}=\frac{f(x|\theta)f(\theta)}{f(x)}$と書けるはずです。
2.事後分布から得たい情報とは
さて、頑張って観測した情報$x$から新たに事後分布$f(\theta|x)$を手に入れました。しかし、その結果どのようなデータが手に入ったかを知るにはその特徴量(平均、分散)なども知りたいはずです。そのためには例えば平均値であれば$$\hat{\theta}_{eap}=\int\theta f(\theta|x)d\theta$$という積分を求めなければなりませんが、この積分が難しい場面が多々あるのです(なんで難しいのかは一旦おいておきましょう)。
まあその積分が難しかったとしてもですよ、**もしここで事後分布$f(\theta|x)$に従う平均値$\mu$,分散 $\sigma$(母数と呼びます)などをランダムに抜き出してたくさんもらえたら？**たとえば$\mu={1.51,1.48,1.39,1.52,1.53...,1.50}$と1000個くらい与えられたとしたら、そのもらった1,000個の平均値の平均をとってしまえばそれが事後分布の平均としても良さそうです。でも、そんなものを都合よくもらえるはずがない...どうしよう(泣)
3.ランダムに母数をゲットする仕組みを作る
さて、もう一度話をまとめましょう。何かしらの母数$\theta$(例えば平均値$\mu)$が従う分布$f(\theta)$が手元にあったところにデータ$x$を手に入れたことによって事後分布$f(\theta|x)$を手に入れました。しかし、その推定量$\hat{\theta}{eap}$を簡単に求めることが出来ませんでした。そのため、$f(\theta|x)$に従う乱数をいっぱい手に入れればその平均を事後推定量$\hat{\theta}{eap}=\int\theta f(\theta|x)d\theta$としてもいいよね！というのが今の段階の発想です。
だから、事後分布$f(\theta|x)$からいっっっぱい値をとってきたいのです。だから、この関数から自由にたくさん値をとりつつもその分布が良くとる値をたくさんサンプリングして、あんまり取らない値はあんまサンプリングしないようにもしたいですよね！
4.事後分布が良くとる値を優先的に取る
さて、ここで詳細つり合い条件$$f(\theta|\theta^{'})f(\theta^{'})=f(\theta^{'}|\theta)f(\theta)$$
について考えましょう。これは比の式に直すと$$f(\theta):f(\theta^{'})=f(\theta|\theta^{'}):f(\theta^{'}|\theta)$$
とかけます。これは、$\theta$をとる確率$f(\theta)$が$\theta^{'}$をとる確率$f(\theta^{'})$より高いのであれば$\theta^{'}$から$\theta$に推移する確率$f(\theta|\theta^{'})$も高いということを意味しています。つまり、$\theta$がメジャーな値であるときには$\theta^{'}から$$\theta$に移りやすいという式なわけです❕
これって3で考えた、「分布が良くとる値をたくさんサンプリングして、あんまり取らない値はあんまサンプリングしないようにもしたい」というワガママを叶えてくれるものな訳です！
5.さあ、メトロポリスヘイスティング法が始まります
ここで、乱数をとってきたい事後分布$f(\theta|x)$を$f(\theta)$と簡略的に表現させてください(混乱しやすいので！)。
先に述べた、詳細つり合い条件にあった$f(~|~)$という関数ですが、遷移核と呼ばれます(マルコフ連鎖の話は調べてみてください...)がこんな遷移核を見つけることもまた困難なわけです。したがって、その代わりに僕らが簡単にわかる関数(正規分布など)を提案分布$q(\theta)$と呼び、これを$f(~|~)$の代わりに用います。
でも、そんなことをすると$$f(\theta|\theta^{'})f(\theta^{'})=f(\theta^{'}|\theta)f(\theta)$$という等式が成り立たなくなります。なのでこの等号をを成立させるために次のようにします。
$$f(\theta|\theta^{'})=cq(\theta|\theta^{'})$$$$f(\theta^{'}|\theta)=c^{'}q(\theta^{'}|\theta)$$このように係数を導入すればきちんと詳細つり合い条件を満たせます。書き換えると$$cq(\theta|\theta^{'})f(\theta^{'})=c^{'}q(\theta^{'}|\theta)f(\theta)$$となります。でもここで、この係数$c,c^{'}$が２つもあるの嫌なので両辺$c^{'}$で割ります。そのとき$r=\frac{c}{c^{'}}, r^{'}=\frac{c^{'}}{c^{'}}=1$として$$rq(\theta|\theta^{'})f(\theta^{'})=r^{'}q(\theta^{'}|\theta)f(\theta)$$と書き換えます。
6.いよいよMH(メトロポリスヘイスティング)アルゴリズム
- まず、ある母数$\theta^{(t)}$が与えられたとき、次の乱数$\theta^{(t+1)}$の候補として提案分布$q(|\theta^{(t)})$から乱数$a$を発生させます。この$a$と$\theta^{(t)}$を比べるのがMHアルゴリズムです。そもそもこの目的は、事後分布$f(\theta)$がとりやすい値を優先的にサンプリングすることなので、最初に与えられた$\theta^{(t)}$が新たに発生した$a$よりも確率の高い値なら$\theta^{(t+1)}=\theta^{(t)}$としてもう一回提案分布から乱数をもらえばいいのです。その一方で、新たに発生した$a$の方が確率が高い値なら$\theta^{(t+1)}=a$とするべきですよね！
2.さて、上に述べたことより$\theta^{(t)}$から$a$に推移しやすいのかどうかを判断するには詳細つり合い条件を見返せばよいことを思い出してください。
$$rq(a|\theta^{(t)})f(\theta^{(t)})=r^{'}q(\theta^{(t)}|a)f(a)$$
さて、もし$\theta^{(t)}$から$a$に推移しやすい($a$の方が起こる確率が高い)のであれば$$q(a|\theta^{(t)})f(\theta^{(t)})>q(\theta^{(t)}|a)f(a)$$となり、詳細つり合い条件を満たすためには$r$が0~1の範囲にあることになります。この$r$による補正がかかって初めて詳細つり合い条件が満たされることになります。つまり、確率$r$で$q(a|\theta^{(t)})f(\theta^{(t)})>q(\theta^{(t)}|a)f(a)$が成立($\theta^{(t)}$から$a$に推移して$\theta^{(t+1)}=a$とする)し、確率$1-r$で$\theta^{(t+1)}=\theta^{(t)}$とします。逆に、$$q(a|\theta^{(t)})f(\theta^{(t)})<q(\theta^{(t)}|a)f(a)$$のときは、確率補正は$r^{'}=1$なので必ず$a$を受容して$\theta^{(t+1)}=a$とします。
【まとめ】
以上のアルゴリズムをまとめると、「提案された候補点$a$を確率$min(1,r)$で受容($\theta^{(t+1)}=a$)し、さもなくばその場に留まる($\theta^{(t+1)}=\theta^{t}$)ことを繰り返す」といえます。
7.さいごに
なるべく、高度な数式を振り回さずに進めてきましたが何となく概念を掴んでいただいて今後の書籍や他の記事を読むための土台としてもらえると幸いです。また、適宜他の記事も作っていこうと思いますのでよろしくお願いいたします。いつか、医学系の人や医療メーカーの人たちで統計や機械学習などを勉強するコミュニティが出来るといいですね...