【ゼロからできるMCMC】メトロポリス法の失敗例について【Chapter 4】

ゼロからできるMCMCを現在読んでいます。
修士課程で量子モンテカルロ法を使った数値計算を行っていたのでMCMCについてはある程度知っているのですが初心に戻って基礎から勉強し直しています。
実は現在のお仕事でもSAを使って計算することがあり、MCMCを使ってSAを動かしたりもしています。

Chpater4ではメトロポリス法について説明がされています。 p56の失敗例を自分でも試してみました。
テキストでは1次元の標準正規分布をメトロポリス法を使ってシミュレートしています。
注意すべき点はk番目の座標$x^{(k)}$からk+1番目の座標$x^{(k+1)}$に遷移するときの幅$\Delta x$の値域の設定を誤ると間違った確率分布になることです。
($x^{(k+1)} = x^{(k)} + \Delta x$)
$\Delta x \in [-0.5, 0.5]$では成功しており、$\Delta x \in [-0.5, 1]$では失敗しています。
前者では以下の詳細釣り合いの条件が成立しており、後者では破れています。

$P(x)T(x\to x^{\prime}) = P(x^{\prime})T(x^{\prime} \to x)$

改めて$x^{\prime} = x + \Delta x$とすると、$x\to x^{\prime}$に遷移する場合は$\Delta x$の幅だけ変化する必要があります。
また$x^{\prime} \to x$の場合は$-\Delta x$だけ変化する必要があります。
これら両方の遷移確率が釣り合っていることが必要になります。
$\Delta x \in [-0.5, 0.5]$であれば$0.5$の幅の範囲で上記の条件を満たしています。
一方で$\Delta x \in [-0.5, 1]$の場合、$\Delta x=1$の時は逆の$-\Delta x=1$も必要になりますが$\Delta x$の最小値が$-0.5$であるため遷移できません。
($P(x^{\prime} \to x)=0$になってしまいます。)
このため上の詳細釣り合いの条件が成立しなくなります。
$\Delta x$に設定する値域が非対称の場合は上記の理由により詳細釣り合いが満たされません。
とりあえずテキストにも記載されている様に$\Delta x\in [-c, c]$と値域を対称的に設定しておけば良いと思います。

テキストで紹介されている成功例・失敗例をpythonで実装してみました。
c言語で書いたコードは著者のgithubで公開されています。
MCMC-Sample-Codes/Gaussian_Metropolis.c at master · masanorihanada/MCMC-Sample-Codes · GitHub

以下のコードは私のgithubのリポジトリでも公開しています。

github.com

まずはモジュールのインポート。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

メトロポリス法を行うMCMCのコードです。
初期位置から最終位置までに遷移した状態がdataとして返ってきます。

niter: MCMCのiteration回数
x: 初期位置
step_left: $\Delta x$の下限値
step_right: $\Delta x$の上限値
seed: シード

def do_MCMC(niter=10**7, x=0.0, step_left=-0.5, step_right=0.5, seed=1111):
    naccept = 0
    data = []
    np.random.seed(seed)
    for i in range(niter):
        backup_x = x
        action_init = 0.5 * x**2
        # dx in [step_left, step_right]
        dx = (step_right - step_left) * np.random.rand() + step_left
        x += dx
        action_fin = 0.5 * x** 2
        if np.exp(action_init - action_fin) > np.random.rand():
            # accept
            naccept += 1
        else:
            x = backup_x
        data.append(x)
    return data

比較用の正規分布を生成するためのコード。

center: 正規分布の平均値
left: 用意するデータの下限値
right: 用意するデータの上限値

def get_normal_dist(center=0.0, left=-4, right=4):
    width = (right - left)
    num_grids = 10*5
    diff = width / num_grids
    pos = [left + diff * i for i in range(num_grids)]
    norm = np.sqrt(2 * np.pi)
    normal = [1.0/norm*np.exp(-0.5 * (left + diff * i - center)**2) for i in range(num_grids)]
    return pos, normal

まずは成功例である$\Delta x \in [-0.5, 0.5]$で計算を行うためのコードと結果になります。

data = do_MCMC()
pos, normal = get_normal_dist()
plt.plot(pos, normal)
plt.hist(data, bins=50, density=True)
plt.show()
plt.savefig('./fig/succeeded_iter10_7.jpg')

f:id:butch416:20200804020554j:plain — 成功例

青い線がget_normal_distで生成した真の標準正規分布です。
オレンジ色のヒストグラムがMCMCで生成した分布であり両者は一致していることが分かります。

次に失敗例である$\Delta x \in [-0.5, 1.0]$で計算を行います。
step_right=1.0に設定することで$\Delta x$の上限値を変更できます。

data = do_MCMC(step_right=1.0)
pos, normal = get_normal_dist(right=7)
plt.plot(pos, normal, label='Analytic(mean=0)')
plt.hist(data, bins=50, density=True, label='MCMC')
pos, normal = get_normal_dist(center=np.mean(data), right=7)
plt.plot(pos, normal, label='Analytic')
plt.legend()
plt.savefig('./fig/failed_iter10_7.jpg')
plt.show()

f:id:butch416:20200804022758j:plain — 失敗例

青い線とオレンジ色のヒストグラムは先ほどと同様です。
明らかに両者にはズレがあるので$\Delta x \in [-0.5, 1.0]$では標準正規分布を再現できないことが分かります。
とはいえオレンジ色のヒストグラムは正規分布っぽい形状をしている気がします。
なのでMCMCで得られたデータの平均値を平均値として推定した正規分布を緑色の線で重ねてみました。
微妙に推定した正規分布からもズレが生じている気がします。

以上より$\Delta x$の値域を適切に設定しないと目的の分布が得られないことが確かめられました。

butch’s blog

メモ置き場。

【ゼロからできるMCMC】メトロポリス法の失敗例について【Chapter 4】