新時代の MCMC を迎えるために

A Blog Entry on Bayesian Computation by an Applied Mathematician

$$

$$

「新時代の MCMC を迎えるために」統数研オープンハウス（タップで PDF 閲覧）

以下は，5月24日 10:30~12:30 （コアタイム：10:30~11:10）に行われた 2024年度統数研オープンハウス ポスターセッション（掲載 No. E1）に於て発表されたポスターに関する解説記事です．

1 導入

1.1 MCMC 小史

現状，HMC (Hamiltonian Monte Carlo) という約 40 年前に提案された MCMC 手法が，Stan などの確率的プログラミング言語のデフォルト MCMC 手法として採用されています．¹

この手法はもともと (Duane ほか, 1987) が場の量子論に特化した Monte Carlo 法として提案したものであったところを，(Neal, 1994) が一般の統計モデルに適用可能な形式に翻訳する形で提案されたものでした．

ということで，HMC は，オリジナルの MCMC が物理学者 (Metropolis ほか, 1953) に由るように，物理学において着想された MCMC 手法であったのです．

そのHMC が，提案から 40 年目を迎える前に，更なる効率的な手法によって代替されようとしています．

そのきっかけ (Peters と de With, 2012) も，やはり，物理学（正確には物質科学）からの着想でした．

1.2 MCMC とは何か？

MCMC とは，確率変数をシミュレーションする際に用いられる汎用的アルゴリズムです．

一様分布や正規分布などの名前がついた分布ではない場合，どのようにすればその分布に従う確率変数をシミュレーションできるのか？は，古くからの問題でした．

実際，「MCMC では空間を探索するマルコフ連鎖を構成し，その足跡を辿るとちょうど確率変数のシミュレーションになっている」と種明かしを聞いても，「なぜそのような回りくどい方法を使うのか？」「もっと良い方法はないのか？」と思っても当然でしょう．

ですが，MCMC を，発明された経緯を辿り，物理学の問題意識から見てみると，実は極めて自然な発想に思えてくるかもしれません．

以降，MCMC の起源である物理系のシミュレーション（第 2 節）を例に取り，分子動力学法（第 2.1 節），Metropolis 法（第 2.2 節）を復習します．

分子動力学法の出力（第 2.1 節）²

Metropolis 法の出力（第 2.2 節）

これを基礎として，近年提案された非対称な MCMC 手法（第 3 章），そして最新の連続時間 MCMC 手法（第 4 章）を紹介します．

非対称 MCMC 法（Lifted Metropolis 法）の出力（第 3 章）

連続時間 MCMC 法（Zig-Zag サンプラー）の出力（第 4 章）

1.3 自己相関・軌跡の一覧

Metropolis 法の自己相関関数（第 2.2 節）

非対称 MCMC 法（Lifted Metropolis 法）の自己相関関数（第 3 章）

連続時間 MCMC 法（Zig-Zag サンプラー）の自己相関関数（第 4 章）

Metropolis 法の軌跡（第 2.2 節）

非対称 MCMC 法（Lifted Metropolis 法）の軌跡（第 3 章）

連続時間 MCMC 法（Zig-Zag サンプラー）の軌跡（第 4 章）

2 MCMC の起源

よりみち：どうして MCMC が必要だったのか？

(Metropolis ほか, 1953) では，温度 \(T\) 一定の条件下で \(N\) 粒子系をシミュレートし，任意の物理量 \(F\) に対してその相空間上の平均 \[ \langle F\rangle=\frac{\int Fe^{-\frac{E}{kT}}dp}{\int e^{-\frac{E}{kT}}dp} \] を効率的に計算する汎用アルゴリズムが提案された．これが現在では Metropolis 法と呼ばれている．³

(Metropolis ほか, 1953) では \(N\) が数百になる場合を考えており（時代を感じるスケール感），当然愚直な数値積分は現代の計算機でも実行可能ではない．そこで Monte Carlo 法を考えることになるが，当時 Monte Carlo 法といえば，一様乱数を用いた計算法の全般を指し，具体的には \(\langle F\rangle\) を重点サンプリング推定量 \[ \widehat{F}=\frac{\sum_{n=1}^NF(\omega)e^{-\frac{E(\omega)}{kT}}}{\sum_{n=1}^Ne^{-\frac{E(\omega)}{kT}}} \] で推定することを指した．⁴

しかしこれでは，高エネルギーな状態・低エネルギーな状態を全く区別せず，状態 \(\omega\in\Omega\) を完全に一様に生成するため，その分だけ非効率である．

これを低減することが出来れば Monte Carlo 法の更なる効率改善に繋がる．こうして，Gibbs 分布 \(\frac{1}{Z}e^{-\frac{E}{kT}}\) から直接的サンプリングする方法が模索されたのである．

(Metropolis ほか, 1953) では，エネルギー \(E\) を持つ系の Boltzmann-Gibbs 分布 \(\frac{1}{Z}e^{-\frac{E}{kT}}\) から直接サンプリングする方法が探求されました．

ここでは簡単のため，１粒子が次のようなポテンシャル \[ U(x)=\frac{x^2}{2}+\frac{x^4}{4} \] に従って運動する場合を考えましょう：

図 1: ポテンシャル \(U\) のプロット

このポテンシャルに関する Boltzmann-Gibbs 分布 \(\pi\,\propto\,e^{-\beta U}\) は次のような形になります：⁵

図 2: ポテンシャル \(U\) が定める Botlzmann-Gibbs 分布のプロット

２つのプロットを見比べると，低エネルギー状態ほど出現確率が高く，エネルギーが上がるにつれて急激に出現確率が下がることがわかります．以降，\(\beta=1\) としましょう．⁶

2.1 分子動力学法

統計力学によれば，\(\beta=1\) で定まる温度とポテンシャル \(U\) を持つ Boltzmann-Gibbs 分布 \(e^{-U}\) は，温度 \(T=\frac{1}{k_B\beta}\) を持つ熱浴に接している力学系を，長時間シミュレーションして時間平均を取ることでサンプリングできるはずです．

このように，力学に基づいて物理過程を数値シミュレーションをすることを通じてサンプリングを達成する方法を 分子動力学法 といいます．

これを実際にやってみます．図 1 で定めたポテンシャルを持つ粒子を考えます．⁷

続いてこれを温度 \(T=\frac{1}{k_B\beta}\) を持つ熱浴と相互作用させます．例えば，ポテンシャル 1 の \(x=0\) の位置に半透性の壁を置き，確率 \(1/2\) でこの温度 \(T\) の壁の粒子と弾性衝突するとします．（残りの確率 \(1/2\) では衝突せずに通過する）．

壁の粒子の速度は Maxwell の境界条件から与えられるものとすれば，次のようにして粒子の位置 \(x\) がシミュレートできます：⁸

import numpy as np

np.random.seed(2024)

def U(x):
    return x**2/2 + x**4/4

def pi(x):
    return np.exp(-U(x))

def md(num_samples, initial_state, initial_velocity, timestep=0.1):
    samples = [initial_state]
    current_state = initial_state
    current_velocity = initial_velocity
    current_time = 0

    for _ in range(num_samples - 1):
        proposed_state = current_state + current_velocity * timestep
        current_time += timestep
        if current_state * proposed_state < 0 and np.random.rand() < 1/2:
            current_velocity = (-1) * np.sign(current_velocity) * np.sqrt(-2 * np.log(np.random.rand() + 1e-7))
        else:
            current_velocity = current_velocity - ( current_state + current_state ** 3 ) * timestep
            current_state = proposed_state
        samples.append(current_state)

    return np.array(samples)

# サンプル数と初期条件を固定
num_samples = 10000
initial_state = 0.0
initial_velocity = 1.0

samples_MD = md(num_samples * 10, initial_state, initial_velocity, timestep=0.01)

図 3: 分子動力学法からのサンプル

この方法は極めて収束が遅く，イテレーション数を \(10^6\) 以上に取らないと目標分布 \(e^{-U}\) の良い近似とならないことを思い知りました（上図も \(10^6\) サンプルで生成しています）．なお，以降の MCMC 法ではいずれもイテレーション数は一桁少ない \(10^5\) としています．

たしかに，目標分布 \(e^{-U}\) に収束しそうですね．

2.2 Metropolis 法

もちろん，分布 \(e^{-U}\) をサンプリングするために，必ずしも背景にある物理過程まで戻ってシミュレーションをする必要はありません．

そこで，シミュレーションは簡単なランダムウォークで行い，その結果を適切に修正することで目標分布に収束させる方法が (Metropolis ほか, 1953) で考えられました．

(Metropolis ほか, 1953) の手法は，現在では random walk Metropolis-Hastings 法と呼ばれます．

この背後の物理現象から離陸する一歩が，分子動力学法と MCMC 法とを分けるものでした．

def metropolis(num_samples, initial_state, verbose=False):
    samples = [initial_state]
    current_state = initial_state

    accept = []

    for _ in range(num_samples - 1):
        proposed_state = current_state + np.random.uniform(-2,2)
        acceptance_ratio = pi(proposed_state) / pi(current_state)
        if np.random.rand() < acceptance_ratio:
            current_state = proposed_state
            accept.append(True)
        samples.append(current_state)

    if verbose:
        rate = len(accept) / num_samples
        print(f'acceptance rate : {rate}')

    return np.array(samples)

図 4: Metropolis 法からのサンプル

サンプル数は分子動力学法の \(1/10\) であるにも拘らず，目標分布 \(e^{-U}\) の良い近似を得ています．

一般に，MCMC からのサンプルの質の良さは，自己相関関数 を見ることで評価できます．⁹

Metropolis 法の自己相関関数を計算してみると，横軸の Lag が大きくなればなるほど Autocorrelation の値は小さくなっています．

図 5: Metropolis 法の自己相関関数

図 6: Metropolis 法の軌跡

上図は Metropolis 法で構成される Markov 連鎖の軌跡を表しています．行ったり来たりしているのがわかります．棄却率は５割弱です．¹⁰

2.3 統計学への応用

こうして MCMC が発明されれば，すぐにイノベーションとして理解されたかというとそうではありませんでした．

この Metropolis の手法が極めて賢いシミュレーション手法であることは一目瞭然でも，一般の確率分布からのサンプリングに使える汎用アルゴリズムになっているという抽象的な観点が得られるまでには時間を要しました．これを成し遂げたのが (Hastings, 1970) でした．¹¹

さらに，Hastings のこの結果も見過ごされたと言って良いでしょう．真にMCMC を統計学界隈に広め，現代におけるベイズ統計学の興隆の契機となったのは階層モデリングにおける Gibbs サンプリングの有用性を強調した (Gelfand と Smith, 1990) だと言われます．¹²

当時，代替手法としては複雑な数値アルゴリズムしかなかったベイズ統計学において，MCMC は汎用的で実装も容易であることが周知され，ベイズ統計学が普及するきっかけとなりました．

3 非対称化への試み

3.1 対称性という制約

ここでもう一度 Metropolis 法の軌跡 図 6 を見てみましょう．

図 6 Metropolis 法の軌跡

最初の 50 サンプルしか表示していませんから，運が悪いとうまく見つからないかもしれませんが，「一度歩んだルートを，その後すぐに逆に戻ってしまう」という事象が発生しやすいことが観察できますでしょうか？

これを 対称性 (reversibility) または 可逆性 と言います．

Metropolis 法は構成上，この対称性を持つことが必要ですが，対称であるが故に一箇所に長時間とどまってしまうことが多くなります．

その結果，対象分布が複雑で多峰性を持つ場合は，もっといろんなモード（峰）からもサンプリングをしてほしいのに，長時間１つの峰から離れられずにいることがあります．

コーヒーに砂糖を溶かすことを考えてみましょう．砂糖の粒が拡散するのに任せておくと，最終的には均一に溶けるでしょうが，莫大な時間がかかります．スプーンで混ぜるなどして，砂糖が元の場所にとどまらずに移動し続けるようにすれば，はるかに速く平衡状態に到達できるでしょう．

これが 非対称化 のアイデアです．数ある Metropolis 法の改良の方向の中でも，この対称性を破るという試みは特に注目されてきました．

3.2 リフティング

Metropolis 法を非対称化するアプローチに，リフティング (Chen ほか, 1999) と呼ばれる方法があります．

これは，元々の状態空間を２つの「モード」\(+1\) と \(-1\) に分裂させ，\(+1\) のモードではひたすら右側に，\(-1\) のモードではひたすら左側に移動するようする方法です．

２つのモード \(+1,-1\) の間を遷移する確率を調整することで，最終的な不変分布は変わらないようにすることができます．

こうすることで，対称性を破り，一度「この方向に行く！」と決めたら行き続けるようにしながら，収束先は変わらないように変更することが出来るのです．

実際に Metropolis 法に適用した Lifted Metropolis-Hastings 法 (Turitsyn ほか, 2011) を実装してみましょう：

def lifted_metropolis(num_samples, initial_state, verbose=False):
    samples = [initial_state]
    current_state = initial_state
    lifting_variable = 1
    accept = []

    for _ in range(num_samples - 1):
        delta = np.random.uniform(0,2)
        proposed_state = current_state + lifting_variable * delta
        acceptance_ratio = pi(proposed_state) / pi(current_state)

        if np.random.rand() < acceptance_ratio:
            current_state = proposed_state
            accept.append(True)
        else:
            lifting_variable = (-1) * lifting_variable

        samples.append(current_state)
    
    if verbose:
        rate = len(accept) / num_samples
        print(f'acceptance rate : {rate}')

    return np.array(samples)

新しく追加されたリフティング変数 \(\sigma\in\{\pm1\}\) に依存して，\(\sigma=+1\) の場合には右方向に，\(\sigma=-1\) の場合は左方向にしか提案を出さない MH 法と見れます．

図 7: 非対称 Metropolis 法からのサンプル

図 8: 非対称 Metropolis 法の自己相関関数

自己相関関数を見ると，Metropolis 法よりも急速に減衰していることがわかります．むしろ，過減衰のように自己相関関数が負になっていることもあります．

これは，一度「この方向に行く！」と決めたら行き続けるように設計したために，正の値が出たしばらくあとは負の値が，負の値が出たしばらくあとは正の値が出やすいようになってしまっているためです．

したがってこれは１次元の分布を考えていることに起因するため，殊更問題とすべきではないでしょう．

図 9: 非対称 Metropolis 法の軌跡

3.3 リフティングの有用性

今回のような単純なポテンシャル \(U\) （図 1） だけでなく，統計力学における磁性体のモデルである Curie-Weiss モデルのハミルトニアン

\[ H_n(x)=-\frac{d\beta}{2n}\sum_{i,j=1}^nx_ix_j-h\sum_{i=1}^nx_i, \] \[ h\in\mathbb{R},x\in\{\pm1\}^n,\quad n,d=1,2,\cdots \]

が定める Boltzmann-Gibbs 分布 \(e^{-H}\) に対する Lifted Metropolis-Hastings も，単純な Metropolis-Hastings 法よりも効率的であることが知られています．¹³

具体的には，モデルのパラメータ数 \(n\) に対して，緩和時間を \(\sqrt{n}\) のオーダーだけ改善することが，(Turitsyn ほか, 2011) では数値実験で，(Bierkens と Roberts, 2017) では理論的に検証されているのです．

4 新たな MCMC

4.1 背後の物理現象からの更なる離陸

第 2 章で，分子動力学法（第 2.1 節）から，提案分布を背後の物理現象とは全く関係ないランダムウォークとすることで，Metropolis 法（第 2.2 節）は一気に効率的なサンプリング法となったことを見ました．

しかし，Metropolis 法はまだ思考が物理に引っ張られているのかも知れません．平衡統計力学において，ミクロの状態は等価で，ミクロなダイナミクスは可逆と考えられます．その前提が，知らず知らずのうちにまだ埋め込まれたままだと言えるでしょう．

そこで，スプーンでかき混ぜるように，遷移を非対称にすることで，より効率的なサンプリング法となることを前章 3 で見ました．

ここでは，さらに暗黙の思い込みから解き放たれようとします．それは，シミュレーションするにあたって，必ずしも離散時間ステップに囚われる必要はない ということです．

もう一度，Lifted Metropolis-Hastings 法の軌跡を見てみましょう：

図 9 非対称 Metropolis 法の軌跡

この軌跡から得られる情報のほとんどは，「どこで折り返したか？」です．

ですから，この軌跡をシミュレーションするにあたって，一歩一歩採択-棄却のステップを繰り返す必要はなく，「どこで折り返すか？」を先に計算できてしまえば，あとは好きなステップサイズで切り出してサンプルとすれば良いのです．

実は，「折り返す地点だけを効率的に計算する」ことが可能であり，それが 連続時間 MCMC のアイデアです．

4.2 連続時間 MCMC

Lifted Metropolis-Hastings の適切な連続時間極限 \(\Delta t\to0\) を考えることで，「折り返す」という事象（が起こった回数）は Poisson 過程に従うことが導けます．

すると，「折り返す」事象が起こるまでの待ち時間 (interarrival time) は指数分布に従うことがわかります．これに基づいて，「折り返す」事象が起こる時刻を計算し，そこまでの軌跡を直線で補間すれば，Lifted Metropolis-Hastings 法（の連続時間極限）の軌跡が模倣できることになります．

最終的に得られる過程は，ランダムな時刻に「折り返す」事象が起こり，その間は確定的な動き（等速直線運動）をするというもので，このような過程を 区分確定的 Markov 過程 (PDMP: Piecewise Deterministic Markov Process) と呼びます．

このような PDMP は，Lifted Metropolis-Hastings 以外にも種々の MCMC 法の極限から見つかっており，その中でも特に有名なのが次の Zig-Zag sampler です：

import math

def zigzag(num_samples, initial_state, step=1):
    samples = [initial_state]
    trajectory = [initial_state]
    current_state = initial_state
    lifting_variable = 1
    t = 0

    while t < num_samples * step:
        state_event = lifting_variable * np.sqrt(-1 + np.sqrt( 1 - 4 * np.log(np.random.rand()) ))
        t_event = t + np.abs(state_event - current_state)
        for _ in np.arange(np.ceil(t/step)*step, np.ceil(t_event/step)*step, step):
              samples.append(current_state + lifting_variable * (_ - t))
        current_state = state_event
        trajectory.append(current_state)
        lifting_variable = (-1) * lifting_variable
        t = t_event

    return np.array(samples), np.array(trajectory)

samples_zigzag, trajectory_zigzag = zigzag(num_samples, initial_state, step=2)
plt.figure(figsize=(3.5, 3))
plt.hist(samples_zigzag, bins=50, density=True, alpha=0.7, color=minty)
plt.show()

図 10: Zig-Zag サンプラーからのサンプル

図 11: Zig-Zag サンプラーの自己相関関数

自己相関関数は，Lifted Metropolis-Hastings 法と同様に急激に下がって負の値に突き抜けたあとは，少し振動が残っているのがわかります．

全３サンプラーの比較は第 1.3 節をご覧ください．

次の軌跡を見て分かる通り，モードである \(x=0\) を中心に激しく往復するので，直後のサンプルとは負の相関が出やすいようです．

図 12: Zig-Zag サンプラーの軌跡

連続時間極限 \(\Delta t\to0\) をとっているということは，「極めて小さいステップサイズでの random walk Metropolis 法（第 2.2 節）」に相当します．従って，一度折り返したら，原点 \(x=0\) を超えるまでは絶対に棄却されません．

そのため，このように往復するような軌跡が得られます．

4.3 連続時間 MCMC の美点

前節では，必ずしも PDMP 法である Zig-Zag sampler が，Lifted Metropolis-Hastings 法より，自己相関関数の観点で良いとは言い切れないことを見ました．

しかし，今回の設定は１次元という特殊な条件下であることを考慮に入れる必要があります．

１次元なので Zig-Zag サンプラーは行き来することしか出来ていませんが，¹⁴ ２次元以上，特に高次元の場合は，Zig-Zag サンプラーは極めて効率的に状態空間を探索できることが期待されます．

例えば，標準正規分布に対する２次元での軌跡は次の通りです：

図 13: Zig-Zag Sampler in \(\mathbb{R}^2\)

そして何より，軌道が効率的な空間の探索に有利であるだけでなく，正確なサブサンプリングを取り入れることが可能です．すなわち，ほとんどの他手法と違って，バイアスを導入することなく，データの一部のみを用いてアルゴリズムを走らせることができます．

したがって，従来の MCMC 法が採択-棄却のステップにおいて尤度を評価する必要があり，データサイズ \(n\) に対して \(O(n)\) の計算量を要するのに対して，\(O(n)\) に比例する焼き入れのステップを除けば，\(O(1)\) の複雑でほとんど i.i.d. なサンプルを得ることができます (Bierkens ほか, 2019)．

5 終わりに

こうみると，MCMC は物理学の問題意識から生まれた手法でありながら，背後の物理現象を模倣することから離れていくことで，計算手法としての効率を高めていく一途を辿っていることがわかります．

そう見ると，新時代の大規模データと大規模モデルが課す MCMC の次なる脱皮は，連続時間 MCMC で間違いないような気がしてくるのですが……．まだ筆者にははっきりとは見えてきません．

また本稿では１つの流れしか取り上げておらず，例えば HMC とその非対称化がどのような位置づけにあるかもまだ考慮中です．

統計力学，統計学，機械学習が交差するなんとも面白いテーマです．

5.1 執筆のきっかけ

本ポスター，そして本解説記事の執筆のきっかけは，MLSS でのポスター発表 で連続時間 MCMC のことが機械学習の界隈では全く知られていないことを知ったことと，情報統計力学の研究集会に出席し，連続時間 MCMC が物理学の方からも研究されていることを知ったことでした．

統計界隈では（本稿で解説しました通り） PDMP や連続時間 MCMC と呼ばれる手法は，物理学界隈では event-based simulation，rejection-free と呼ばれる手法群の１つとして活発に研究されていました．

全く同じ問題を解こうとしているのに，用語法が全く異なることに驚かされました．

２つの分野の相互理解と知見の交換が進むことを目指し，これからも研究していきたいと考えます．

5.2 アルゴリズムのステップサイズについて

コードからわかります通り，提案分布のスケールは次のようになっています：

３手法のスケーリング

手法	提案分布	平均移動距離15
MH	\(\mathrm{U}([-2,2])\)	\(1\)
LMH	\(\mathrm{U}([0,2])\)	\(1\)
Zig-Zag	なし	\(2\)

ですので，平均した隣接（提案）サンプル間距離について，Zig-Zag サンプラーはズルをしているともいえます．

しかし，必ずしもアンフェアな比較をしていたわけではありません．

もし，３手法で計算量を揃えるならば，Zig-Zag サンプラーにとっての１回のループは方向転換をするまでであり，２つの方向転換の間の平均距離はだいたい \(1.93\) になります．

diffs = np.diff(trajectory_zigzag)
abs_diffs = np.abs(diffs)
mean_abs_diff = np.mean(abs_diffs)
print(mean_abs_diff)

1.9232887144203465

MH, LMH については，ほとんど最適なスケーリングになるように調節してありますが，Zig-Zag サンプラーにおいては最適なスケーリングという概念は存在しません．あらかじめ 図 12 の軌跡があり，どれくらい距離を空けてサンプルとするか，という問題しかありません．

ここでは，アルゴリズムのループ１回が１回のサンプリングになるように，３つの手法を揃えて比較の根拠としました．

なお，このように，Zig-Zag サンプラーなどの連続時間 MCMC と，従来の離散時間 MCMC の直接の性能比較は，微妙な問題が多く，筆者もまだ十分に説明できる準備はありません．

しかし，次年度以降のオープンハウスで，より詳細な解説を行う予定です．

Zig-Zag サンプラーのサンプリングステップを \(1\) や \(1.5\) にすると，たしかに自己相関は悪化します．

6 参考文献

本稿では，用いたコードの導出を一切触れませんでした．これについては，文献 (Tartero と Krauth, 2024) をご参照ください．非調和振動子を系にとり，正準集団とみなすことで，分子動力学法，メトロポリス法からそのリフティングまで，種々のサンプラーを同じ題材で比較するアイデアをもらいました．こんなにわかりやすい解析ができるのかと心底驚きました．

続いて，Metropolis-Hastings 法 → 非対称 MCMC → 連続時間 MCMC という発展の過程を，背後の物理過程の模倣からの離陸という視点で統一的に捉えることが出来るということは，(Turitsyn ほか, 2011) の魅力的なイントロダクションで気付かされました．本文献はリフティングによる MCMC の非可逆化を抽象的に定式化して数値実験で検証したものであり，「ねじれ詳細釣り合い条件」を導入した点で，アイデアの宝庫といえます．

リフティングによる Metropolis 法の非対称化について，(酒井佑士, 2017) は貴重な日本語文献です．当該文献の第３章（の第２節）にここで紹介した内容が詳しくまとめられています．

第 4.3 節で紹介しましたように，Zig-Zag サンプラーを用いたサンプリングではそのスケーラビリティが魅力です．この点については，Zig-Zag サンプラーが提案された論文 (Bierkens ほか, 2019) でも，前面に押し出して解説されています．

参考文献

Bierkens, J., Fearnhead, P., と Roberts, G. (2019). The Zig-Zag Process and Super-Efficient Sampling for Bayesian Analysis of Big Data. The Annals of Statistics, 47(3), 1288–1320.

Bierkens, J., と Roberts, G. (2017). A Piecewise Deterministic Scaling Limit of Lifted Metropolis-Hastings in the Curie-Weiss Model. The Annals of Applied Probability, 27(2), 846–882.

Chen, F., Lovász, L., と Pak, I. (1999). Lifting Markov chains to speed up mixing. Proceedings of the Thirty-First Annual ACM Symposium on Theory of Computing, ページ 275–281. New York, NY, USA: Association for Computing Machinery.

Duane, S., Kennedy, A. D., Pendleton, B. J., と Roweth, D. (1987). Hybrid Monte Carlo. Physics Letters B, 195(2), 216–222.

Gelfand, A. E., と Smith, A. F. M. (1990). Sampling-Based Approaches to Calculating Marginal Densities. Journal of the American Statistical Association, 85(410), 398–409.

Gelman, A., Roberts, G. O., と Gilks, W. R. (1996). Efficient Metropolis Jumping Rules. Bayesian Statistics 5: Proceedings of the Fifth Valencia International Meeting. Oxford University Press.

Hastings, W. K. (1970). Monte Carlo Sampling Methods Using Markov Chains and Their Applications. Biometrika, 57(1), 97–109.

Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., と Teller, E. (1953). Equation of State Calculations by Fast Computing Machines. The Journal of Chemical Physics, 21(6), 1087–1092.

Neal, R. M. (1994). An Improved Acceptance Procedure for the Hybrid Monte Carlo Algorithm. Journal of Computational Physics, 111(1), 194–203.

Peters, E. A. J. F., と de With, G. (2012). Rejection-free Monte Carlo Sampling for General Potentials. Physical Review E, 85(2).

Robert, C., と Casella, G. (2011). A Short History of Markov Chain Monte Carlo: Subjective Recollections from Incomplete Data. Statistical Science, 26(1), 102–115.

Roberts, G. O., Gelman, A., と Gilks, W. R. (1997). Weak Convergence and Optimal Scaling of Random Walk Metropolis Algorithms. The Annals of Applied Probability, 7(1), 110–120.

Tartero, G., と Krauth, W. (2024). Concepts in Monte Carlo Sampling. American Journal of Physics, 92(1), 65–77.

Turitsyn, K. S., Chertkov, M., と Vucelja, M. (2011). Irreversible Monte Carlo algorithms for Efficient Sampling. Physica D-Nonlinear Phenomena, 240(5-Apr), 410–414.

Vasdekis, G., と Roberts, G. O. (2022). A note on the polynomial ergodicity of the one-dimensional Zig-Zag process. Journal of Applied Probability, 59(3), 895–903.

酒井佑士. (2017). マルコフ連鎖モンテカルロ法における詳細つり合い条件の破れの効果と応用 (phdthesis). 東京大学. 読み込み から https://repository.dl.itc.u-tokyo.ac.jp/records/50422

脚注

1. Hamiltonian Monte Carlo の名称は (Neal2011-HMC?) からで，元々は Hybrid Monte Carlo と呼ばれていました．分子動力学法 (Molecular Dynamics) と MCMC のハイブリッド，という意味でした．Stan で実装されている MCMC アルゴリズムについては こちら を参照．

2. 特に収束が遅く，他の手法と比べてイテレーション数を 10　倍にしています．

3. 統計学界隈では (Hastings, 1970) を入れて，Metropolis-Hastings 法とも呼ばれる．

4. ただし，配置 \(\omega\in\Omega\) は空間内にランダム（一様）に粒子 \(N\) 個を配置することで生成することとする．

5. １粒子系なので相互作用はなく，\(E=U\)．

6. \(\beta=1\) と約束することは，系の温度を \(T=k_B^{-1}\) に固定することにあたります．

7. 図 1 で定めたポテンシャルを持つ力学系には，代表的なものは（非調和）振動子や，あるいは \(U\) の形をした谷を行ったり来たりするボールを考えても構いません．

8. 詳しい議論は (Tartero と Krauth, 2024) をご参照ください．大変教育的な入門です．

9. 自己相関関数が大きいほど，その Markov 連鎖を用いて構成した Monte Carlo 推定量の漸近分散が大きくなります．加えて，自己相関関数の裾が重すぎると，例えエルゴード性を持っており大数の法則が成り立とうとも，中心極限定理が成り立たなくなります．換言すれば，\(n^{-1/2}\) よりも遅い収束レートになってしまいます．

10. 一般には，ランダムウォーク MH 法において，採択率を \(0.2\le\alpha\le0.4\) 前後に抑えるのが良いとされています (Roberts ほか, 1997)．これは状態空間の次元が無限に漸近する設定下での，漸近論的な結果ですが，低次元の場合でも極めて良い指標になることが (Gelman ほか, 1996) で実証されています．また今回も，１次元であるにも拘らず，たしかに棄却率が半分を越さないほうが，自己相関が小さくなる傾向が確認されました．しかし今回は対象分布の裾が極めて軽いので，あまり大きなムーブは要らず，ステップサイズの最大値を \(2\)，採択率は \(0.5\) 近くにしました．他の手法，LMH と Zig-Zag もステップサイズの最大値が \(2\) になるように統一しました．

11. “While (Metropolis ほか, 1953) proposed the use of MCMC sampling to compute particular integrals in statistical mechanics, it was the Hastings paper that elevated the concept to a general one, and introduced it to the broader statistics community.” (Martin+2023-history?) 3.5節．

12. (Martin+2023-history?) 4節，(Robert と Casella, 2011, p. 102) など．

13. ただし，\(e^{-H}\) が多峰性を示す低温領域では，LMH の方が効率的であるというはっきりとした理論保証はまだありません．

14. 今回は対象分布の減衰が極めて激しかったために差が現れにくかったのだと考えられます．Zig-Zag サンプラーは１次元でも，（広い設定の下で）（そして特に目標分布の裾が重いときに）ランダムウォーク・メトロポリス法や Metripolis-adjusted Langevin algorithm よりも速い収束レートを持ちます (Vasdekis と Roberts, 2022)，

15. 提案分布の．実際の軌跡の１ステップでの移動距離とは違う．

著作権

Copyright Hirofumi Shiba 2024. All Rights Reserved