【ゼロからできるMCMC】練習問題4.2(途中でステップ幅が変化するMCMC)【Chapter 4】

ゼロからできるMCMCの練習問題4.2について考えます。
シミュレートしたい真の分布は

$P(x) = \frac{1}{2\sqrt{2\pi}}\biggl(e^{-\frac{x^2}{2}} + e^{-\frac{(x-100)^2}{2}} \biggr)$

です。分布の山が$x=0, x=100$に集中していることを考慮して、偶数番目のステップと奇数番目のステップで遷移幅$\Delta x \in [-c, c]$の値域cを変更しています。
具体的には偶数番目のステップでは$c=1$、奇数番目のステップでは$c=100$に設定しています。

テキストには上記の更新方法でも分布をシミュレートできる理由について解説されています。
実際の数値計算でも問題ないことを確かめました。
ですが私自身よく理解していないこともあるのでその点についてまとめます。

まずはステップの偶奇を表す変数$y=0, 1$を導入しています。
そして状態を$(x, y)$で指定することにしています。
$y$は$0, 1, 0, 1,\cdots$と周期的に変化していきます。
MCMCの必要条件として非周期性がありましたがこの場合は成立しなくなります。

$y$については上記の様に非周期性が成り立ちませんが現在は$(x, y)$で考えているので必ずしもそう言えるのでしょうか？
$(x, y)$で考えると、$(x^{(0)}, 0), (x^{(1)}, 1), (x^{(2)}, 0), (x^{(3)}, 1),$と$x^{(k)}$も同時に周期的になっていなければいいんじゃないかと思いました。

テキストでは$y$を改めて$y=0\to 1 \to 0$と再定義することでこの問題点を回避しています。
これは常に$y$が同じ状態しか遷移しないので2回の遷移で元に戻る、3回の遷移で元に戻る、・・・と周期が1になるため非周期性が満たされると理解しました。

これに伴い遷移確率にも修正が必要になります。
$x\to x^{\prime}$の遷移について考えます。
$x$を偶数番目とすると、$x^{\prime}$は奇数番目のステップになります。

詳細釣り合い条件を書き出すと

$P(x) T_{c=1}(x\to x^{\prime}) = P(x^{\prime}) T_{c=100}(x^{\prime} \to x)$

になります。
これが成立していれば嬉しいのですが成立しません。
$T_{c=100}(x^{\prime} \to x)$に注目すると$x=x^{\prime}+\Delta x$で$\Delta x \in [-100, 100]$になります。
$\Delta x = 100$になった時、逆の遷移$x\to x^{\prime}$を考えると偶数番目から奇数番目の遷移なので$1$の幅でしか遷移できないので遷移確率が0になります。
以上の理由より詳細釣り合い条件が敗れているので修正が必要になります。
新しい遷移確率

$\begin{align} T(x\to x^{\prime \prime}) &= \int dx^{\prime} T_{c=1}(x\to x^{\prime}) T_{c=100}(x^{\prime} \to x^{\prime\prime}) \\ \tilde{T}(x\to x^{\prime \prime}) &= \int dx^{\prime} T_{c=100}(x\to x^{\prime}) T_{c=1}(x^{\prime} \to x^{\prime\prime}) \end{align}$

を導入すると、以下の条件が成立することが分かります。

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

詳細釣り合い条件っぽいですがそうではないみたいです。
この式の導出を以下で行います。
まずは$c=1, 100$それぞれに固定された世界で成立する詳細釣り合い条件を書き出しておきます。

$\begin{align} P(x) T_{c=1}(x\to x^{\prime}) &= P(x^{\prime}) T_{c=1}(x^{\prime} \to x)\\ P(x) T_{c=100}(x\to x^{\prime}) &= P(x^{\prime}) T_{c=100}(x^{\prime} \to x) \end{align}$

これを利用すると、

$\begin{align} P(x) T(x\to x^{\prime}) &= \int dx^{\prime\prime} \biggl(P(x)T_{c=1}(x\to x^{\prime\prime})\biggr) T_{c=100}(x^{\prime\prime} \to x^{\prime}) \\ &= \int dx^{\prime\prime} \biggl(P(x^{\prime\prime})T_{c=1}(x^{\prime\prime}\to x)\biggr) T_{c=100}(x^{\prime\prime} \to x^{\prime}) \\ &= \int dx^{\prime\prime} T_{c=1}(x^{\prime\prime}\to x) \biggl( P(x^{\prime\prime}) T_{c=100}(x^{\prime\prime} \to x^{\prime}) \biggr) \\ &= \int dx^{\prime\prime} T_{c=1}(x^{\prime\prime}\to x) \biggl( P(x^{\prime}) T_{c=100}(x^{\prime} \to x^{\prime\prime}) \biggr) \\ &= P(x^{\prime}) \int dx^{\prime\prime} T_{c=100}(x^{\prime} \to x^{\prime\prime}) T_{c=1}(x^{\prime\prime}\to x) \\ &= P(x^{\prime}) \tilde{T}(x^{\prime} \to x) \end{align}$

となり修正した遷移確率で詳細釣り合い条件もどきが成立していることが分かります。
詳細釣り合い条件は成立していませんが釣り合い条件が成立していることがテキストで説明されています。

$\begin{align} \int dx^{\prime} P(x) T(x\to x^{\prime}) = \int dx^{\prime} P(x^{\prime}) T(x^{\prime}\to x) \end{align}$

つまり、

$\begin{align} P(x) = \int dx^{\prime} P(x^{\prime}) T(x^{\prime}\to x) \end{align}$

それぞれの$c$でも定常分布に収束することが分かっているので

$\begin{align} P(x) = \int dx^{\prime} P(x^{\prime}) T_{c}(x^{\prime}\to x) \end{align}$

となることを利用して修正した繊維確率$T$でも釣り合い条件が成立することを説明しています。

私が疑問に思ったのが実装するときには$T_c, T(\tilde{T})$どちらの遷移確率で計算するべきかという点です。
テキストの説明的には

と定義した繊維確率$T(\tilde{T})$を使うべきだとは思いますが積分計算が間に入っているので計算がめんどくさそうです。
単純に偶数番目で$c=1$、奇数番目で$c=100$とした$T_c$を用いて計算したコードと結果を記載します。
コード自体はこちらの記事で使用したものとほとんど同じです。
以下のコードは私のgithubのリポジトリでも公開しています。

github.com

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

def get_normal_dist(center1=0.0, center2=100.0, left=-4, right=120):
    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) * 2
    normal = [1.0/norm*(np.exp(-0.5 * (left + diff * i - center1)**2) + np.exp(-0.5 * (left + diff * i - center2)**2)) for i in range(num_grids)]
    return pos, normal

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

data = do_MCMC()
pos, normal = get_normal_dist()
plt.plot(pos, normal, label='Analytic')
plt.hist(data, bins=1000, density=True, label='MCMC')
plt.legend()
plt.savefig('./fig/practice4_2_iter10_7.jpg')
plt.show()