マルコフ連鎖モンテカルロ法 - メトロポリス・ヘイスティングス法（１）

概要

マルコフ連鎖モンテカルロ法とは、サンプリングしたい分布（目標分布）が不変分布となるようなマルコフ連鎖を設定することでからサンプリングする手法です。マルコフ過程で、ある適当な推移核を定めたとき、各状態に遷移する確率が不変分布になる例を見ました。今回は逆で、不変分布が既知、つまり目標分布であるです。そしてが不変分布となるような推移核を求めることが目標です。

マルコフ連鎖モンテカルロ法を実現するアルゴリズムの１つが、メトロポリス・ヘイスティングス法です。なおマルコフ連鎖モンテカルロ法は略してMCMC法と表記されることが多いようです。

詳細釣り合い条件

マルコフ過程では、状態が離散的なものを考え推移核として行列を用いましたが、連続的な分布では、状態からへ遷移する推移核をと表します。 マルコフ連鎖が不変分布に収束する十分条件が次の詳細釣り合い条件です。

例えばだったとしたとき、となりますから、生起しやすい状態（）から生起しにくい状態()へ遷移する確率は低く、その逆の確率は高いことを意味しています。下図のようなイメージです。

ここで式(1)の両辺をで積分すれば、左辺はになりますから、

となります。式(2)は、あらゆる点から推移核をによってへ遷移する確率の和がそのものであることを示しています。つまり詳細釣り合い条件が成り立つ推移核によってサンプリングを続ければ、それは分布に収束するはずです。*1

メトロポリス・ヘイスティングス法のアルゴリズム

概要で書いた通り、いまサンプリングしたい目標分布はですから、これは既知です。推移核がわかっていれば良いのですが不明な状態です。そこでサンプリング可能な既知の提案分布を用います。この提案分布は欲しい推移核ではありませんから、式(1)の詳細釣り合い条件を満たしておらず、例えば

などのように等号が成立していません。ここで、式(3)の等号が成立すればからサンプリングできるのですから、次式のようにある確率で補正し、からへ遷移する確率を下げてやれば良いはずです。

つまりからへ遷移したとき、それを確率で採択すれば詳細釣り合い条件が満たされる、つまりからのサンプリングが実現できるということになります。 式(4)より、は

です。なお式(3)の等号は逆かもしれず、その場合は式(5)のが1以上となります。確率が大きく補正されるわけですから、無条件で採択すれば良いはずです。

以上をまとめます。

メトロポリス・ヘイスティングス法のアルゴリズム

1. 目標分布を、提案分布をとする。
2. 提案分布からをサンプリングする。
3. の値を計算する。
4. が1以上ならを採択する。が1未満なら、その確率でを採択する。採択されなければ同じ値（）をサンプリング点とする。
5. 手順２に戻り、再び次の候補点をサンプリングする。

なお、のような提案分布を用いれば、

です。例えばランダムウォークと呼ばれるものはが成り立ちます。

式(5)で考えれば、イメージはしやすいアルゴリズムですよね。既知の提案分布から次の候補点をサンプリングし、それが目標分布において確率が高い方向なら採択、低い方向なら、低い度合いに合った確率で採択する、という感じでしょうか。ただし採択されなかったら捨てるのではなく、同じ値がサンプリングされたとするところが棄却サンプリングの採択部分とは違います。また、初期値は任意に決めますが、発生確率が低いところから出発してしまうと、本来サンプルされにくい部分からいきなり始まってしまいます。初期値への依存を薄めるため、最初のいくつかのサンプルは捨ててしまうのが良いようで、これをバーンインと呼ぶようです。

サンプリング実験

コードを書いて実験してみました。とりあえず簡単のために１次元の分布です。提案分布は、正規分布を用いたランダムウォークです。目標分布はです。

ちゃんとサンプルできました。

次に以下のような目標分布で試してみました。

MCMC法の性質上、サンプル点は前回の値の影響を受けます。したがって、このような確率の高い山が複数ある分布からサンプリングをすると、１つの山から集中的にサンプルし、ある確率で別の山に移るという挙動になるはずです。それを確認したのが次のアニメーションです。50サンプルごとにヒストグラムを描画しています。やはり各山から交互にサンプルしているような動きになっていました。

しかしサンプル数を十分に大きくすれば、目的の分布に近いものになりました。もっと意地悪な形状の分布にするとサンプリングがうまくいかない場合などもあるかもしれません。

次回は2次元の目標分布からサンプリングの実験をしてみたいと思います。（次回：マルコフ連鎖モンテカルロ法 - メトロポリス・ヘイスティングス法（２））　

今回のコードです。提案分布の分散を変えるとサンプリングが失敗することもありました。これも1次元程度ならいろいろ調整しながらやれますが多次元だとうまくいったかどうか、どう判断すればいいんでしょうかね。

# マルコフ連鎖モンテカルロ法（メトロポリスヘイスティングズ法）のテスト
import numpy as np
import matplotlib.pyplot as plt

def randomwalk(x):
    z = np.random.normal(0,0.5)
    return x + z

# サンプルしたい関数
def p(x):
    if x < -1 or 1 < x:
        return 0
    return np.cos(0.5*np.pi*x)**2

# MCMC
def MCMC(x):
        # サンプル候補
        z = randomwalk(x)
        if p(z)/p(x) > np.random.uniform(0, 1):
            return z
        else:
            return x


# ランダムシードを固定
np.random.seed(0)

# サンプル数
N = 100000

# バーンイン
B = 1000

# 初期値
x = 0.1

samples = np.empty(N)

# MCMC
for i in range(N):
    x = MCMC(x)
    samples[i] = x

plt.hist(samples[B:], bins=50, normed=True, alpha=0.5)
x = np.linspace(-1, 1, 500)
y = np.vectorize(p)(x)
plt.plot(x, y)
plt.show()

*1:記事の中でからへ遷移したり、からに遷移したりいろいろな書き方をしてしまってややこしいですが、どっちでも同じことです