MCMC名前はかっこいいですよね。使いこなせるとカッコいいんですが、なかなか初見で理解するのが難しいです。以前自分で勉強したときも段階的に理解できるような記事がなくて苦労した記憶があるので、今回はそのような記事を作るというのが目的です。
本日の内容です。
- MCMCとは？
- そもそもMCMCで何がしたいの？
- ベイズの定理と共役事前分布
- 事後分布が解析不能な形な時の方法について
- 現在の情報を使ったサンプリング
- MCMCのアルゴリズムを見てみる
- MCMCで事後分布の情報を取れる理由
- 実際の挙動を見てみる
- マルコフ連鎖モンテカルロ法って結局なんなの？
- さらなる学習のために
MCMCとは？
Markov Chain Monte Carlo法（マルコフ連鎖モンテカルロ法） のことです。一般的な記事ではマルコフ連鎖とモンテカルロ法という風に分解され〜から入ってマルコフ連鎖とモンテカルロ法のそれぞれの解説が始まる記事が多いのですが、あまり本質ではないので今回は一旦置いておきます。
そもそもMCMCで何がしたいの？
ベイズモデリングにおいて、事後分布の性質について知りたい時に使う方法です。ベイズモデリングのお気持ちをざっくり話すと、「適当な」事前分布を用意し、情報が与えられることで情報を更新し、事後分布を生成し、この事後分布を使って意思決定や分析、検定などを行う、ということを行います。つまり、事後分布についての情報ができるだけ欲しいということになります。
事後分布の生成にはベイズの定理を使います。
ベイズの定理と共役事前分布
一般的には事後分布を知りたいときにはベイズの定理を用います。
を事前分布とすると、データ
が与えられた時の事後分布
は
の形で与えられます。細かいところはいつか記事になります。今理解しておくべきことは、このベイズの定理によって事前分布から事後分布が出ること、この計算は分母の積分が難しく、適切にモデリングしなければ解くことができない、ということの2点です。
一般的に事前分布に共役事前分布と呼ばれる特殊な形を使うと解析的に解けますが一般的には難しいです。一般的に積分が難しかったり、一般的な漸化式が解けなかったりとかと同じようなものです（実際ほぼ同じ）
事後分布が解析不能な形な時の方法について
事後分布が解析不能な時は、サンプリングしてなんとか情報が取れないかを検討します。これがMCMCの大事な発想につながっていきます。
サンプリングって何？
標本抽出を行うことです。標本抽出とは、母集団がどのような分布をしているかわからない時に一部の集団をとってくることです。例えば選挙とかの例が有名ですね。選挙に行ったうち一部の人にアンケートを実施して選挙行動の動向を聞くことで人口全体の選挙行動を推測して事前の趨勢を把握したり、開票が始まったタイミングで一部しかまだ開票が進んでいない段階でも全体の結果を推測して当確を出したりということができます。
今回は説明のために母集団が分かっている分布を用いて色々試してみます。
サイコロの出目のサンプリング
・例えば、サイコロを初めてみた人がいて、サイコロが、1/6ずつ等しく同じ目が出ると分かっていなかった場合、10000回サイコロを振るという試行をしてみます。
実際に振って結果をメモするのはめんどくさいのでPythonを使って試行をシミュレーションしてみます。
結果は以下の通りです。
サイコロが等しく出目が出る一様分布に従うということを知らなかったとしてもこの結果を見るとなんとなく「出目が等しく出そう」という結論を導くことができるのではないでしょうか？
正規分布からのサンプリング
サイコロの出目のような1, 2, 3, 4, 5, 6 といった離散でなくとも、連続分布からもサンプリングができます。
ただし連続分布において1.49359203といった特定の値を取る確率は0なので、サンプリングの結果は1以上2未満に入ったなどと「幅」を持たせて集計します。今回のサンプリングは標準正規分布からのサンプリングを100万回行い-4から4までの区間を100分割して、それぞれの区間に入った個数を集計します。
これを元々の母集団である標準正規分布のグラフと比較すると次のようになります。
サンプリングによりかなりの確度で元の分布の性質を反映できていることがわかるのではないでしょうか。
現在の情報を使ったサンプリング
2状態、確率1/2の状態遷移
今度は、下図のような状態遷移図に従う、1秒ごとに確率1/2で状態を遷移し、確率1/2でその状態に留まるモデルを考えます。初期状態では状態1にいるものとします。これも理論上解くことができ、高校で習う確率漸化式の要領でn秒後の確率分布を求めることができますし、十分時間が経ったあとの状態1,2それぞれにいる確率は対称性から1/2と求めることができます。
さて、十分賢くない中学3年生に戻り、確率漸化式のことを学んでいない状態に戻ったとして、この十分経ったあと状態Aにいる確率がどれほどかを知りたくなったとします。
そうすると、コイントスをしながら実際にこの移動を10000回シミュレーションして、実際に何回Aにいたかを見てみることでなんとなく確率がわかることが想像できるでしょうか。コインを10000回トスしてメモしたくないので上記の設定でシミュレーションを行い、Pythonくんに結果だけ出してもらいます。
試しに最初の10回の遷移を書き出すと
['1', '2', '2', '1', '2', '1', '2', '1', '1', '2']
というような感じになります。
10000回のシミュレーションは以下のようになります。
ほぼ1/2ずつになっていますね。
このように時系列データに対しても長時間シミュレーションしてサンプリングすることにより収束した際の確率分布の情報を取ることができます。
MCMCのアルゴリズムを見てみる
ついに、本題に入ります。MCMCのアルゴリズムは次のようになっています。ここではMCMCのアルゴリズムのうち最も初歩的なMH（メトロポリス・ヘイスティングス法）について扱います。
①初期値を適当に設定する（定義域内ならなんでもいい）
②ループ回数を決める(多ければなんでもいい)
③遷移確率を適当に定義する（計算しやすければなんでもいい）（この状態移動に関する確率分布を提案分布と言います）
④以下を設定した回数ループする
❶現在位置を確認し、②で決めた遷移状態に従って次に移動する予定値を決める
❷実際にその予定値に移動するか、現在値に止まるかどうかを事前分布 と適当に決めた遷移ルール
で決まる値を用いて確率的に決める
❸移動したか移動していないかによらず、今いるところをメモし、サンプリングする
ご覧のようにMCMCというのは①〜③まではめちゃめちゃ適当です。MCMCの本質と言っても過言ではない一番重要なエッセンスは❷です。❷がなければ先ほどの2状態の遷移となんら変わりません。「状態遷移をして、その現在位置をサンプリングする」というのが先ほどの例ですが、MCMCは「確率的に状態遷移をしてその現在位置をサンプリングする」というのに過ぎません。
MCMCで事後分布の情報を取れる理由
MCMCでは❷によって、「事後分布において事後確率が高いところにいるとき、遷移を保留し、事後確率が低いところにいるとき、遷移を積極的に行うように移動する」というようになっています。
これにより事後確率が高いところで頻回のサンプリングができるように、事後確率が低いところではあまりサンプリングが行われないようになります。
これはMCMCでは詳細釣り合いという条件として付加されているのでその辺を参照ください。そしてこのアルゴリズムはループ回数を極限まで増やすことでサンプリングの結果と事後分布が一致することが知られています（ただし収束までのループ回数は保証されていませんが。。。）
実際の挙動を見てみる
事前分布が正規分布、正規分布から得られるデータが得られる時、これらは共役事前分布であり、事後分布は正規分布になることが解析的に求められます。
ここでMCMCを実際に実装してみて挙動を比較してみます。
折線が移動経路の挙動、右側のグラフが集計した事後分布のサンプリングになります。実際の挙動を見てみると、上に外れても、下に外れてもすぐに真ん中に帰ってきて、真ん中にいる間はしばらく留まっているということがわかります。
マルコフ連鎖モンテカルロ法って結局なんなの？
ここまでMCMCについて、一切の言葉の説明をせずに突っ走ってきました。マルコフ連鎖モンテカルロ法はマルコフ連鎖とモンテカルロ法に分割され、それを合成した言葉なわけで、実際ある意味MCMCを表しているということなのですが、それから解説を始めるとよくわからなくなってしまうということが往々にしてあるんだと思います。MCMCのわからなさって結局ここにあったんだと思います。本稿の最後に後方視的に意味を解説していきたいと思います。
マルコフ連鎖って？
マルコフ連鎖とは、現在の状況によってのみ次の移動条件が決まる確率過程のことです。現在の状態さえ分かれば、未来の状況が予測できる、過去の状態は未来を予測するのに不要ということをマルコフ性と言います。マルコフ性に従う確率変数の列がマルコフ連鎖です。
例えば、2項間漸化式、あれがマルコフ連鎖です。先ほどの1/2の確率で状態を移動する過程もマルコフ連鎖と言えます。MCMCにおいてはシミュレーションにおいて次の移動先を決めるのに現在の位置しか関係ないのでマルコフ性があるということになります。先ほどの「事後分布において事後確率が高いところにいるとき、遷移を保留し、事後確率が低いところにいるとき、遷移を積極的に行うように移動する」とかも現在の位置のみを参照していて、過去の位置については見ていないというのが大事なポイントでしたよね。遷移するかどうかもそうですが、具体的な遷移先の座標についても現在の位置のみを参照して決まっています。
モンテカルロ法って？
モンテカルロ法とは、確率的な試行（乱数を用いた試行）によってシミュレーションを行い推定を行う方法のことです。有名なものだとビュフォンの針などでしょうか。均一幅で引かれた直線群に対して針を投げると円周率が求まると言ったものです。針を投げるという確率的操作で円周率の推定ができることがモンテカルロ法に相当します。今回の説明で言えば❷における、移動するかしないかが確率的に決まることによって事後分布が推定される、というところにモンテカルロ法が使用されています。
マルコフ連鎖モンテカルロ法って？
結局、ベイズ推定において「現在の情報のみを参照しながら遷移先を決定し」、「確率的に遷移するかしないかを決める」ことで事後分布をシミュレーションする、という方法になります。それ以上の発展的アルゴリズムや、具体的な数式の内容については成書をご参照ください。
さらなる学習のために
いくつかベイズ統計のおすすめの本を置いた記事を書いておきます。ご自身のレベルに合わせてチョイスください。
syleir.hatenablog.com
