論文「Markov Chain Monte Carlo Method without Detailed Balance」の実装

かなり革新的なマルコフ連鎖モンテカルロ法(MCMC)の論文があるという話を友人に聞いて、実際に読んで見たらこれはすごいぞと。・・・ということでさっそく簡易実装。
書いたプログラムは書き捨てのつもりで書いたんでクラスとか使ってない。

その論文というのはこれ。

Physical Review Letterに掲載されているということなので、もの凄くインパクトのある結果だということですね、はい。
内容としては詳細釣り合いの条件を課さないでも棄却率が最小になるようにうまくマルコフ連鎖を設計しているというもの。
彼らはポッツモデルを使ってこの手法の優位性を合わせて検証している。

このBLOGでの実装は適当なオレオレモデル*1
マルコフ連鎖を複数本用意してシミュレーションを実行、各マルコフ連鎖の同ステップでの”ある状態の確率分布”をプロット。
結果はこちら。

どちらも不変分布*2にはいってるっぽいけど、提案手法のほうが早くいくとかそういう形ではないっぽいな・・・論文ではある物理量の期待値の収束が速いと言っているのか・・・?
この辺は実装を間違っているのかなんなのかもうちょい調べた方が良さそう。
できれば金融工学でいう所のSVモデルあたりのキャリブレーションとかに応用できないか試したい。

# -*- coding: utf-8 -*-
import numpy
import matplotlib.pyplot
#不変分布(比率のみ, 最初の要素に一番ウェイト高い奴を持ってくる)
INVARIANT_WEIGHT = [5.0] + [float(i+1) for i in range(4)]
#↑から計算される適当なグローバル変数(累積和・状態数)
INVARIANT_WEIGHT_CUM = numpy.cumsum(INVARIANT_WEIGHT)
NUM_OF_STATE = len(INVARIANT_WEIGHT)
#乱数生成機(seedは適当に設定)
GENERATOR_RAND = numpy.random.mtrand.RandomState(100)

#メトロポリス法で次の状態計算する関数
def Metropolis(current):
    prob        = GENERATOR_RAND.random_sample()    
    candidate   = GENERATOR_RAND.randint(0, NUM_OF_STATE)
    weight_candidate = INVARIANT_WEIGHT[candidate]
    weight_current   = INVARIANT_WEIGHT[current] 
    if prob < min(1.0,  weight_candidate / weight_current):
        return candidate
    else:
        return current
#Suwa-Todoらによる手法で次の状態を計算する関数,逆変換法を使ってる
def SuwaTodo(current):
    prob = GENERATOR_RAND.random_sample()
    prob_sum = 0.0    
    for candidate in range(NUM_OF_STATE):
        prob_sum += (v_SuwaTodo(current, candidate) / INVARIANT_WEIGHT[current])
        if prob < prob_sum:
            return candidate 
#Suwa-Todoらの手法で使用されている遷移確率×不変分布の値。論文ではv_ijとなってる
def v_SuwaTodo(i, j) :
    weight_i = INVARIANT_WEIGHT[i]
    weight_j = INVARIANT_WEIGHT[j]
    if j ==0 :
        invariant_cum_j = INVARIANT_WEIGHT_CUM[-1]
    else:
        invariant_cum_j = INVARIANT_WEIGHT_CUM[j - 1]
    delta = INVARIANT_WEIGHT_CUM[i] - invariant_cum_j + INVARIANT_WEIGHT[0]
    return max(0, min(delta, weight_i + weight_j - delta, weight_i, weight_j))
    

if __name__ == '__main__':
    #モンテカルロステップ数。列方向の次元
    NUM_OF_STEP = 40
    #MCMCの鎖の数。行方向の次元
    NUM_OF_CHAIN = 1000
    #初期状態
    initial_state = 3
    #MCMCの各ステップ・各鎖ごとのデータをためこむための行列
    metropolis = numpy.zeros((NUM_OF_CHAIN, NUM_OF_STEP))
    suwatodo   = numpy.zeros((NUM_OF_CHAIN, NUM_OF_STEP))
    #MCMCやる
    for row in range(NUM_OF_CHAIN):
        state_metropolis = initial_state
        state_suwatodo   = initial_state
        for col in range(NUM_OF_STEP):
            state_metropolis = Metropolis(state_metropolis)
            metropolis[row, col] = state_metropolis
            state_suwatodo = SuwaTodo(state_suwatodo)
            suwatodo[row, col]   = state_suwatodo

    #各ステップでの状態ごとの確率分布
    probability_series_metropolis = numpy.zeros((NUM_OF_STATE, NUM_OF_STEP))
    probability_series_suwatodo   = numpy.zeros((NUM_OF_STATE, NUM_OF_STEP))
    for col in range(NUM_OF_STEP):
        for state in range(NUM_OF_STATE):
            probability_series_metropolis[state, col] = len([i for i in metropolis[:, col] if i == state]) / float(NUM_OF_CHAIN)
            probability_series_suwatodo[state, col]   = len([i for i in suwatodo[:, col]   if i == state]) / float(NUM_OF_CHAIN)
    
    #ある特定の状態が出現する確率のステップ(時系列)推移
    matplotlib.pyplot.plot(probability_series_suwatodo[3].T, 'r',  label='Suwa-Todo',linewidth=3)
    matplotlib.pyplot.plot(probability_series_metropolis[3].T, 'b',  label='Metropolis',linewidth=3)
    matplotlib.pyplot.legend()    
    matplotlib.pyplot.show()    

※このコードを使った一切の結果は未保証です

*1:不変分布を適当にリストとして与えてモデルとした

*2:ここではコード内でINVARIANT_WEIGHTとして設定している p = 3 / (5 + 1 + 2 + 3 + 4) = 1/5 = 0.2が不変分布になってる。分子の3はプロット関数の所で適当に設定