はじめに
こんにちわ、研究開発部の島田です。今回はMCMC法の一種である、ランダムウォーク メトロポリス・ヘイスティングス法(MH法)についてのお話です。
これまた小難しそうな名前ですが、実はMCMC法の中では最もシンプルなアルゴリズムとして有名です。今回はこちらを理解した上で実装を行い結果を見てみたいと思います。
おさらい
MH法を理解するにはまずは前回記事の内容が大事になってくるためおさらいをします。
前回記事では、マルコフ連鎖について保育士さんの行動モデルを例にして理解しました。 特にある状態が確率的に変化することを遷移核と呼び、マルコフ連鎖を扱うにあたっては非常に重要な概念であることがわかりました。そして、前回記事では遷移核があらかじめわかっている状態でマルコフ連鎖を使いました。しかし、現実問題は遷移核がわからないケースがほとんどであるため、遷移核をどのように設定すればよいのかという問題を解決する必要があります。
遷移核を設定する際に、考えるべきことの一つとしてどのような条件を満たせば定常分布へ収束するのかです。これを詳細つり合いと言い、詳細つり合いを満たすと確率分布が定常分布になります。
そしてこの詳細つり合いをを満たす遷移核は簡単にはわからないため、代わりに条件つき確率分布と乱数によるサンプリング手法を使って、詳細つり合いになるように補正をかけてあげる、これがMH法の考え方になります。乱数によるサンプリング手法については、下記記事を参考にするとイメージがわきやすいです。
ランダムウォーク メトロポリス・ヘイスティングス法(ランダムウォークMH法)
前項でも説明した通りMH法の考え方として、詳細つり合いをを満たす遷移核は簡単にはわからないため、サンプリングが簡単な遷移核（提案分布）で代用します。そして詳細つり合いを満たすように提案分布をあとで補正をかけます。
さて、ここで気になるのが、提案分布はどのように決めるのか？ということですね。
よく使われるのがランダムウォークという考え方で、数学や物理学の世界で頻繁に登場します。ランダムウォークという面白い名前なのですが、由来となっているのはなんと酔っ払いの千鳥歩きだそうです。これからの忘年会シーズンでよく目にする酔っ払いがフラフラ歩いている姿、そうあれです。そしてランダムウォークを使ったMH法のことを、ランダムウォークMH法と呼びます。
では、酔っ払いの千鳥歩き（ランダムウォーク）を数学的に表現してみます。まずは一番シンプルな一次元のランダムウォークは下図のように表せます。最初に豚の酔っ払いさんがいる地点から左に行くのか右に行くのかは確率的な問題です。
新しい地点は下記のように表せます。
ここで、はランダムウォークが動く幅を表すパラメータです。
は0〜1の正規分布を仮定しています。
さて、ランダムウォークについて理解したところでいよいよランダムMH法の中身について理解していきます。
まずはMH法の一般的な考え方についてですが、先に詳細つり合いを形式化します。
左辺については、例えばある確率が発生します（実際には事後分布）。そうすると、確率
は
に移ります（遷移核）。右辺は今度はある
が発生して、
が
に移ることを示しています。詳細つり合いを満たすということはすなわちこれらの確率の流れが等しくなることになります。
MH法では、サンプリングが簡単な提案分布で代用しますが、このとき詳細つり合いは厳密には満たされないために、
というような式で表せます。
そこで、詳細つり合いを満たすように補正係数を導入します。
これをγについて解くと下のようになります。
すなわち、このという補正係数を導入することで詳細つり合いを満たすような遷移が実現できることになります。
ここで、ランダムウォークの概念を取り入れると、と
は等しくなるはずです。なぜならば上の一次元ランダムウォーク図を見るとわかる通り、豚さんが右に移動する確率と移動された地点から元の場所に移動する確率も等しいはずですよね。そのため、補正係数
は下記のようにシンプルに表せます。
少し説明が長くなりましたが、実際にこのアルゴリズムを実装してみます。
ランダムウォークMH法の実行
今回もベータ分布で確認してみます。まずはベータ分布の確認です。
np.random.seed() a, b = 1.8, 2.9 x = np.linspace(beta.ppf(0.001, a, b), beta.ppf(0.999, a, b), 100) plt.plot(x, beta.pdf(x, a, b))
次にランダムウォークMH法のコアな部分ですが、流れとしては下のフローのようになります。
初期値θを決定する（適当で良い）
ランダムウォークを使って次の新しいθ’を探索
次の条件式の判定。
Step3の条件式が、
・真の場合、現在位置の方が次の位置よりも確率分布の高い位置であると言えるため、で補正する。区間[0, 1]の一様乱数を生成し、
と比べて小さい場合は受理、それ以外の場合は棄却する。
・偽の場合、補正せずに新しい乱数を受理する。
これを実装して、結果を確かめてみます。
def rw_MH(): NUM_MC_STEP = 30000#モンテカルロステップ数 theta = 0.6#ランダムウォークの初期値 E = 0.5#ランダムウォークの幅 theta_mc_step = [theta] for i in range(NUM_MC_STEP): re_theta = theta + E * np.random.randn() if beta.pdf(theta, a, b) > beta.pdf(re_theta, a, b): gamma = beta.pdf(re_theta, a, b) / beta.pdf(theta, a, b) if np.random.rand() < gamma: theta = re_theta else: theta = re_theta theta_mc_step.append(theta) return theta_mc_step
オレンジ色の曲線が元のベータ分布ですが、ランダムMH法によるサンプリングも同じような曲線が得られていることがわかります。
最後に
今回はMCMC法の最もシンプルなランダムウォークMH法を理解し、実装・結果を見てみました。遷移核があらかじめわかっていないケースにおいてもサンプリングが可能になりました。他にもギブスサンプラーやハメルトニアンモンテカルロ法といった様々なMCMC法があるので、是非今後も試してみたいと思います。
参考文献
・RとStanではじめるベイズ統計モデリングによるデータ分析入門（講談社）
