【完全ガイド】MCMC法についてわかりやすく解説|ベイズ推定
こんにちは、青の統計学です。
今回は、ベイズ理論を使ったパラメータ推定手法であるMCMC法(Markov Chain Monte Carlo法 マルコフ連鎖モンテカルロ法)について解説いたします。
ベイズの定理から丁寧に解説していきますので、ベイズの考え方に馴染みがなくてもOKです。
ブックマーク推奨です！
MCMC法|Markov Chain Monte Carlo法
まずは、ベイズ統計学の考え方について学びましょう。
ベイズの定理
ベイズの定理について、具体例を元に解説した、より平易なコンテンツは以下になります。
さて、ベイズ推量の世界では、全ての確率は主観的な確率 (subjective probability) だとされています。
普通、確率というのは対象の不確かさを示す量とされますが、ベイズ推量では不確かさはその対象を観測する人の知識にあると考えており、それを主観確率と呼んでいます。
では、事前確率と事後確率について説明します。
ベイズの定理で使う事前確率 Pr(X) とは、観測者が観測以前に持っている主観確率を確率分布関数で表現したものです。
ちなみに、事前分布に無情報一様分布を設定すると、事後分布と尤度関数が一致するので、MAP推定量は最尤推定量と一致します。
この辺りは頻度論と整合性があって素敵ですね。
一方、事後確率はなんでしょうか？
例えば、観測した結果を \(D\) とします。
観測結果 \(D\) によって、確率変数 \(X\) に関する主観確率が変化します。
観測結果 \(D\) のもとで確率変数 \(X\) に関する主観確率を \(Pr(X|D)\) と書き、事後確率 (posterior) と呼んでいます。
この事後確率分布を求めるのがベイズ推量の目的です。
尤度|lilelihoodとは
ベイズ統計では、尤度という概念が必要です。
詳しく知りたい方は以下のコンテンツをご覧ください。
【AICで使う】KL divergence(カルバック-ライブラー情報量)をわかりやすく解説|python
確率変数 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)}$$
これをベイズの定理と呼びます。
また分母の値を特に エビデンス(evidence) と呼びます。
ベイズ統計では、データを元にしてどんどん事後分布を更新していくアプローチを取ります。
次元の呪い|なぜモンテカルロ法だと推定が難しいか
パラメータの数が増えると、共同事後分布の複雑さが増し、サンプリングが難しくなります。
なので、まずモンテカルロ法というランダムなサンプリング法を考えてみます。
統計検定準1級でも、正方形の中に90度の扇形を入れて、ランダムにサンプリングして、円の面積を求めていくという問題がありました。
このように、乱数を用いた試行を繰り返すことで近似解を求める手法で、確率論的な事象についての推定値を得る場合を特に「モンテカルロシミュレーション」と呼びます
では、単純モンテカルロ積分を見てみましょう。
提案分布を一様分布\(U(0,1)\)として、提案分布から乱数\(X〜N(0,1)\)を発生させていきます。
\(g(x)\)はよくわからないので、それを覆うような提案分布からサンプリングしていくということですね。
$$θ=\int_{0}^{1}g(x)dx$$
モンテカルロ積分は事後確率分布 \(Pr(X|D)\) からサンプル \(X\) を取り出すことで期待値\(E[g(X)]\)を近似値的に評価していきます。
大数の弱法則により、サンプリング回数が無限大になると、\(θ\)に確率収束します。
$$E[g(X)]=\hat{θ}=\frac{1}{n} \int_{0}^{1}g(X_i)dx$$
ランダムサンプルを生成して、それらの平均を取ることで積分値を近似します。
ただし、問題点があります。
1回ごとのサンプリングが独立なので、高次元だと殆ど棄却されてしまうのです。
高次元になると、提案分布に対して採択される的が小さいイメージを持っていただけると良いと思います。
高次元空間では、分布が疎になることが多く、効率的なサンプリングが困難になる可能性があるということですね。
これを次元の呪いと呼びます。
多項式曲線のフィッティングの問題ですね、係数の数が冪乗に増加していき、高次元になると幾何的な直感が当たらないです。
では効率よくサンプリングするには、一期前のサンプリングを使えば良いのではという発想が自然と生まれるはずです。
これが、マルコフ連鎖モンテカルロ法の着想です。
マルコフ連鎖
では、事後確率分布をどう推定すれば良いでしょうか？
次にマルコフ連鎖という考え方を使ってサンプリングをします。
時刻 \(t,t ≥ 0\) において、系列の現在の状態 \(X_t\) だけに依存する分布 \(Pr(X_{t+1}|X_t)\) から次の状態 \(X_{t+1}\)をサンプルして、ランダムなパラメータの列 \({X_0, X_1, X_2, . . . }\) を生成することを考えます。
つまり、\(X_t\) を決めると次の状態 \(X_{t+1}\) はそれより過去の系列 \({X_0, X_1, . . . ,X_{t−1}}\) には依存しないとします。
$$Pr(X_{n+1}=j|X_{n},X_{n-1},…,X_0)=Pr(X_{n+1},X_n=i)$$
このような確率過程をマルコフ連鎖と呼びます。
ここでは系列は時間的に一様で\(t\) に依らないとすると、上のマルコフ連鎖は単一の定常分布 (stationary distribution) に収束することが知られています。
この辺は、統計検定準一級でも頻出ですね！
【最短合格】統計検定準一級のチートシート|難易度や出題範囲について
さて、この定常分布を \(π(X)\) とすると、定常分布\(π(X)\) が正確に事後確率分布 \(Pr(X|D)\) に近づくような操作を行うのが MCMC 法の目的です。
\(t\) が増加するにつれてサンプル点は定常分布からサンプルされるようになり、初期条件の情報は失われます。
つまり定常分布は初期状態には依らないことになります。
→主観によって設定される事前分布の影響は受けにくい！！！
繰り返しになりますがマルコフ連鎖は、各ステップが前のステップのみに依存する性質を持ち、時間が経つにつれて分布が収束します。
これにより、複雑な確率分布からのサンプリングが可能になるというわけです。
具体的には、MCMCは以下のような流れになります。
①初期値\(X_0\)を決める
②乱数で今の\(X\)から新しい\(X_{new}\)を探してくる
③次の条件式を判定する。
尤度の大小関係を調べます。
$$f(X_{new}|D)>f(X|D)$$
④真であれば、状態を更新し、偽であればある程度の確率で受容する。
⑤②-④を繰り返す。
上図のように、尤度が高い部分をフラフラとサンプリングしていくようになります。
詳細釣り合いと遷移核
先ほど、「マルコフ連鎖が定常になる〜」と簡単に説明しましたが、どのような条件を満たせばマルコフ連鎖は定常になるのでしょうか？
結論から申し上げると、下のように一期前と現在の遷移核と事後分布の積が等しくなる場合です。
$$f(X’|X)f(X’)=f(X|X’)f(X)$$
\(f(X|X’)\)：遷移核(移りやすさのようなイメージ)
→パラメータ\(X’\)から\(X\)にどれだけ移りやすいか
\(f(X)\)：事後分布(起こりやすさ)
とはいえ、この条件は十分条件なので、マルコフ過程が定常分布に収束するための必須の条件ではないです。
メトロポリス・ヘイスティングアルゴリズム|Metropolis-Hastings algorism
上で説明したような、詳細釣り合いを満たすための遷移核のバリエーションはいくつかあります。
・メトロポリス・ヘイスティングス法
・ギプスサンプリング
・ハミルトニアンモンテカルロ法
ここでは、MHアルゴリズムをご紹介します。
上の遷移核の議論で思った方も多いと思いますが、「事後分布の形が複雑で直接サンプリングできないという話だったのに、遷移核がわかっているのはなぜ？」となるはずです。
その通りで、遷移核\(f(X’|X)\)を見つけるのは簡単ではなく、遷移核に従うような乱数を発生させることも同様に難しいです。
なので、サンプリングが簡単な遷移核\(q(X’|X)\)で代用することがMHのスタートです。
→これを提案分布と呼びます。
ただ、提案分布は真の遷移核ではないので、詳細釣り合いは満たされないので下のような補正係数\(r\)を導入します。
この\(r\)について解けば良いことになります。
$$rq(X’|X)f(X’)=q(X|X’)f(X)$$
そして、\(q(X’|X)f(X’)>q(X|X’)f(X)\)ならば、\(r\)を使って確率的に補正をしていきます。
一方で、上の条件式が偽であるならば、新しいパラメータを受け入れた方が詳細釣り合いを満たすので受容します。
この試行を繰り返していきます。
ちなみに提案分布が\(q(X’|X)=q(X|X’)\)のようなランダムウォーク(つまり、パラメータの移りやすさが1期前と同じ)であるならば、補正係数\(r\)は一期前と現在の事後分布の比になることがわかります。
何が嬉しいかというと、計算量が少なくなるのに加えて、探索空間の広範囲カバーが挙げられます。
というのも、ランダムウォークは事後分布の異なる領域を探索するのに役立ちます。
局所的な最大値に囚われるリスクを軽減し、より全体的な分布の特徴を捉えることができます。
提案分布からのサンプリングを例に、ヤコビアンを説明したコンテンツはこちら
【数理統計学】ヤコビアンをわかりやすく解説|MCMCでの使用例
