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

マルコフ連鎖モンテカルロ法とは、サンプリングしたい分布が不変分布となるようなマルコフ連鎖を設定することでからサンプリングする手法です。そしてその手法を実現するアルゴリズムの１つが、メトロポリス・ヘイスティングス法です。マルコフ連鎖モンテカルロ法は略してMCMC法と表記されることが多いようです。

先回のマルコフ過程では、状態が離散的なものを考えましたが、サンプリングしたい分布は離散的でなく連続的な分布です。離散的な場合は推移行列をかけることによって次の時刻の状態を表すことができました。連続の場合は、

と表されます。ここで は推移核と呼ばれ、からに推移する確率を表します。時刻に、にいる確率というのは、その1つ前の状態からに推移する全ての場合を積分すればよいということです。

推移核が以下の詳細釣り合い条件と呼ばれる式を満たすとき、は不変分布となるようです。

これはからに遷移する確率と、からに遷移する確率が等しいことを意味しています。この条件は、不変分布となる十分条件であり、必要条件ではありません。マルコフ過程で実験した不変分布もこれを満たしていません。

この辺りの理論はなかなか難解でしっかりと理解できていないのですが、メトロポリス・ヘイスティングス法のアルゴリズム自体はとてもシンプルです。サンプリングしたい分布（目標分布）をとしたとき、この分布の推移核がわかっていれば話は早いのですが、当然不明なので提案分布を用います。提案分布はサンプリング可能な既知のものとします。そしてこの提案分布からサンプリングし、それを目標分布から得られる採択確率に従って採択します。まとめると以下のような手順です。

1. マルコフ連鎖の初期値を適当に決める
2. 提案分布に従ってをサンプリングする
3. 以下の採択確率に従って採択する（棄却されれば再度提案分布からサンプリングする）

ここで、提案分布が詳細釣り合い条件を満たしていれば、式(3)の採択確率は

となります。

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

さて、早速コードを書いて実験してみました。とりあえず簡単のために１次元の分布です。提案分布は、正規分布を用いたものにします。左右対称な正規分布ですから詳細釣り合い条件を満たすはずです。目標分布はです。

ちゃんとサンプルできていますね。次に以下のようなものを試してみました。目標分布はです。棄却サンプリングで実験した分布です。 端っこのほうが目標分布どおりにサンプリングできていません。本来高い確率でサンプリングすべき部分ですが、範囲]を超えてしまう確率も高いので、目標分布に沿った採択率にならなかったということかな。

次回は2次元の目標分布からサンプリングの実験をしてみたいと思います。（次回：マルコフ連鎖モンテカルロ法 - メトロポリス・ヘイスティングス法（２））　本当はもっと数学的な背景をしっかり押さえていきたいんですが、ここに書いた上っ面の理解でちょっと心配。

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

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


# 提案分布
def q(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):
    while True:
        # サンプル候補
        z = q(x)
        # 受理するまでサンプル候補をとる
        if p(z)/p(x) > np.random.uniform(0, 1):
            return z


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

# サンプル数
N = 100000

# バーンイン
B = 1000


# 初期値
x = 0.5

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()