ハミルトンモンテカルロ

ハミルトンモンテカルロアルゴリズム（元々はハイブリッドモンテカルロと呼ばれていました）は、直接サンプリングすることが難しい目標確率分布に分布が収束するランダムサンプルのシーケンスを取得するためのマルコフ連鎖モンテカルロ法です。このシーケンスは、期待値やモーメントなどの目標分布の積分を推定するために使用できます。

ハミルトンモンテカルロは、メトロポリス-ヘイスティングスアルゴリズムのインスタンスに対応し、時間可逆で体積を保存する数値積分器（通常はリープフロッグ積分器）を使用してシミュレートされたハミルトン動力学の進化により、状態空間内の新しい点への移動が提案されます。メトロポリス-ヘイスティングスアルゴリズムでガウスランダムウォーク提案分布を使用する場合と比較して、ハミルトンモンテカルロは、シンプレクティック積分器を使用する場合のシミュレートされたハミルトン動力学の近似エネルギー保存特性により、高い受け入れ確率を維持する遠い状態への移動を提案することにより、連続するサンプル状態間の相関を減らします。相関が減るということは、特定のモンテカルロ誤差に対するターゲット確率分布に関して積分を近似するために必要なマルコフ連鎖サンプルが少なくなることを意味します。

このアルゴリズムはもともと、Simon Duane、Anthony Kennedy、Brian Pendleton、Duncan Rowethによって1987年に格子量子色力学の計算用に提案されました。^[1] 1996年に、Radford M. Nealは、この手法をより広範な統計問題、特に人工ニューラルネットワークに使用できることを示しました。^{[2]しかし、}ベイジアンネットワークの勾配を提供しなければならないという負担により、統計学やその他の定量的分野でのアルゴリズムの幅広い採用は遅れ、2010年代半ばにStanの開発者がHMCを自動微分と組み合わせて実装しました。^[3]

サンプルするターゲット分布が( )であり、サンプルのチェーンが必要であるとします。 ${\displaystyle f(\mathbf {x} )}$${\displaystyle \mathbf {x} \in \mathbb {R} ^{d}}$${\displaystyle d\geq 1}$${\displaystyle \mathbf {X} _{0},\mathbf {X} _{1},\mathbf {X} _{2},\ldots }$

ハミルトン方程式は

${\displaystyle {\frac {{\text{d}}x_{i}}{{\text{d}}t}}={\frac {\partial H}{\partial p_{i}}}\quad {\text{and}}\quad {\dfrac {{\text{d}}p_{i}}{{\text{d}}t}}=-{\dfrac {\partial H}{\partial x_{i}}}}$

ここで、およびはそれぞれ位置ベクトルと運動量ベクトルの 番目の成分であり、はハミルトニアンである。対称かつ正定値の 質量行列を とすると、ハミルトニアンは${\displaystyle x_{i}}$${\displaystyle p_{i}}$${\displaystyle i}$${\displaystyle H}$${\displaystyle M}$

${\displaystyle H(\mathbf {x} ,\mathbf {p} )=U(\mathbf {x} )+{\dfrac {1}{2}}\mathbf {p} ^{\text{T}} M^{-1}\mathbf {p} }$

ここで は位置エネルギーです。ターゲットの位置エネルギーは次のように与えられます。 ${\displaystyle U(\mathbf {x} )}$

${\displaystyle U(\mathbf {x} )=-\ln f(\mathbf {x} )}$

これはボルツマン因子から来ています。指数確率の重みは明確に定義する必要があるため、この定式化ではハミルトニアンは無次元であることに注意してください。たとえば、有限温度でのシミュレーションでは、因子(ボルツマン定数) は直接およびに吸収されます。 ${\displaystyle H}$${\displaystyle \exp \left(-H\right)}$ ${\displaystyle T}$${\displaystyle k_{\text{B}}T}$ ${\displaystyle k_{\text{B}}}$${\displaystyle U}$${\displaystyle M}$

このアルゴリズムでは、リープフロッグステップの数に正の整数、ステップサイズに正の数が必要です。チェーンが にあると仮定します。 とします。まず、ランダムなガウス運動量が から引き出されます。次に、粒子はハミルトン力学の下で時間 の間動きます。これは、リープフロッグアルゴリズムを使用してハミルトン方程式を数値的に解くことによって行われます。リープフロッグアルゴリズムを使用した時間後の位置と運動量ベクトルは次のとおりです。^[4]${\displaystyle L}$${\displaystyle \Delta t}$${\displaystyle \mathbf {X} _{n}=\mathbf {x} _{n}}$${\displaystyle \mathbf {x} _{n}(0)=\mathbf {x} _{n}}$ ${\displaystyle \mathbf {p} _{n}(0)}$${\displaystyle {\text{N}}\left(\mathbf {0},M\right)}$${\displaystyle L\Delta t}$${\displaystyle \Delta t}$

${\displaystyle \mathbf {p} _{n}\left(t+{\dfrac {\Delta t}{2}}\right)=\mathbf {p} _{n}(t)-{\dfrac {\デルタ t}{2}}\nabla \left.U(\mathbf {x} )\right|_{\mathbf {x} =\mathbf {x} _{n}(t)}}$

${\displaystyle \mathbf {x} _{n}(t+\Delta t)=\mathbf {x} _{n}(t)+\Delta tM^{-1}\mathbf {p} _{n}\left(t+{\dfrac {\Delta t}{2}}\right)}$

${\displaystyle \mathbf {p} _{n}(t+\Delta t)=\mathbf {p} _{n}\left(t+{\dfrac {\Delta t}{2}}\right)-{\ dfrac {\Delta t}{2}}\nabla \left.U(\mathbf {x} )\right|_{\mathbf {x} =\mathbf {x} _{n}(t+\デルタ t)}}$

これらの方程式をおよび回適用しておよび を取得します。 ${\displaystyle \mathbf {x} _{n}(0)}$${\displaystyle \mathbf {p} _{n}(0)}$ ${\displaystyle L}$${\displaystyle \mathbf {x} _{n}(L\Delta t)}$${\displaystyle \mathbf {p} _{n}(L\Delta t)}$

リープフロッグ アルゴリズムは、相互作用しない古典粒子の運動に対する近似解です。正確であれば、古典的な位置エネルギー場が存在する場合、各粒子のエネルギーは保存されるため、解はランダムに生成された初期のエネルギー分布を変更することはありません。熱力学的平衡分布に到達するには、粒子が、たとえば周囲の熱浴と何らかの相互作用をする必要があります。これにより、システム全体がボルツマン分布に従った確率で異なるエネルギーを引き受けることができます。

システムを熱力学的平衡分布に近づける方法の 1 つは、メトロポリス-ヘイスティングス アルゴリズムを使用して粒子の状態を変更することです。まず、リープフロッグ ステップを適用し、次にメトロポリス-ヘイスティングス ステップを適用します。

からへの移行は ${\displaystyle \mathbf {X} _{n}=\mathbf {x} _{n}}$${\displaystyle \mathbf {X} _{n+1}}$

${\displaystyle \mathbf {X} _{n+1}|\mathbf {X} _{n}=\mathbf {x} _{n}={\begin{cases}\mathbf {x} _{n} (L\Delta t)&{\text{確率付き }}\alpha \left(\mathbf {x} _{n}(0),\mathbf {x} _{n}(L\Delta t)\right)\\\mathbf {x} _{n}(0)&{\text{otherwise}}\end{cases}}}$

どこ

${\displaystyle \alpha \left(\mathbf {x} _{n}(0),\mathbf {x} _{n}(L\Delta t)\right)={\text{min}}\left(1,{\dfrac {\exp \left[-H(\mathbf {x} _{n}(L\Delta t),\mathbf {p} _{n}(L\Delta t))\right]}{\exp \left[-H(\mathbf {x} _{n}(0),\mathbf {p} _{n}(0))\right]}}\right).}$

完全な更新は、最初に運動量をランダムにサンプリングし（以前の反復とは無関係に）、次に運動方程式を積分し（たとえば、リープフロッグを使用して）、最後にメトロポリス-ヘイスティングスの受け入れ/拒否ステップから新しい構成を取得することから構成されます。この更新メカニズムが繰り返されて、 が得られます。 ${\displaystyle \mathbf {p} }$${\displaystyle \mathbf {X} _{n+1},\mathbf {X} _{n+2},\mathbf {X} _{n+3},\ldots }$

No U-Turn Sampler (NUTS) ^[5]は、自動制御による拡張です。チューニングは重要です。たとえば、1次元の場合、ポテンシャルは、単純な調和振動子のポテンシャルに対応します。大きすぎると、粒子が振動し、計算時間が無駄になります。小さすぎると、粒子はランダムウォークのように動作します。 ${\displaystyle L}$${\displaystyle L}$${\displaystyle {\text{N}}(0,1/{\sqrt {k}})}$${\displaystyle U(x)=kx^{2}/2}$${\displaystyle L}$${\displaystyle L}$

大まかに言えば、NUTS は U ターン条件が満たされるまで、ハミルトン力学をランダムに時間的に前方と後方の両方に実行します。U ターン条件が満たされると、パスからランダムなポイントが MCMC サンプル用に選択され、その新しいポイントからプロセスが繰り返されます。

詳細には、二分木を構築して、蛙跳びのステップの経路をトレースします。MCMC サンプルを生成するために、反復手順が実行されます。スライス変数がサンプリングされます。前方粒子の位置と運動量をそれぞれとします。同様に、後方粒子についても、とととします。各反復では、二分木は、前方粒子を時間的に前方に移動するか、後方粒子を時間的に後方に移動するかをランダムに一様に選択します。また、各反復では、蛙跳びのステップの数は 2 倍に増加します。たとえば、最初の反復では、前方粒子は 1 つの蛙跳びのステップを使用して時間を前方に移動します。次の反復では、後方粒子は 2 つの蛙跳びのステップを使用して時間を後方に移動します。 ${\displaystyle U_{n}\sim {\text{Uniform}}(0,\exp(-H[\mathbf {x} _{n}(0),\mathbf {p} _{n}(0)]))}$${\displaystyle \mathbf {x} _{n}^{+}}$${\displaystyle \mathbf {p} _{n}^{+}}$${\displaystyle \mathbf {x} _{n}^{-}}$${\displaystyle \mathbf {p} _{n}^{-}}$

この反復手順はUターン条件が満たされるまで継続されます。つまり、

${\displaystyle (\mathbf {x} _{n}^{+}-\mathbf {x} _{n}^{-})\cdot \mathbf {p} _{n}^{-}<0\quad {\text{or}}\quad .(\mathbf {x} _{n}^{+}-\mathbf {x} _{n}^{-})\cdot \mathbf {p} _{n}^{+}<0}$

あるいはハミルトニアンが不正確になったとき

${\displaystyle \exp \left[-H(\mathbf {x} _{n}^{+},\mathbf {p} _{n}^{+})+\delta \right]<U_{n}}$

または

${\displaystyle \exp \left[-H(\mathbf {x} _{n}^{-},\mathbf {p} _{n}^{-})+\delta \right]<U_{n}}$

ここで、例えば、. ${\displaystyle \delta =1000}$

Uターン条件が満たされると、次のMCMCサンプルは、二分木によって描かれたリープフロッグパスを均一にサンプリングすることによって得られ、これは次式を満たす。 ${\displaystyle \mathbf {x} _{n+1}}$${\displaystyle \{\mathbf {x} _{n}^{-},\ldots ,\mathbf {x} _{n}(-\Delta t),\mathbf {x} _{n}(0),\mathbf {x} _{n}(\Delta t),\ldots ,\mathbf {x} _{n}^{+}\}}$

${\displaystyle U_{n}<\exp \left[-H(\mathbf {x_{n+1}} ,\mathbf {p_{n+1})} \right]}$

残りの HMC パラメータが適切であれば、これは通常満たされます。

• 動的モンテカルロ法
• モンテカルロ分子モデリング用ソフトウェア
• HMC を実装した確率的プログラミング言語Stan 。
• PyMC は、HMC を実装した確率的プログラミング言語です。
• メトロポリス調整ランジュバンアルゴリズム

• Betancourt, Michael; Girolami, Mark (2015)。「階層モデルのためのハミルトンモンテカルロ」。Upadhyay, Satyanshu Kumar; et al. (eds.)。ベイジアン手法の最新動向とその応用。CRC Press。pp. 79–101。ISBN 978-1-4822-3511-1。
• ベタンコート、マイケル (2018)。「ハミルトンモンテカルロの概念的入門」。arXiv : 1701.02434 [ stat.ME]。
• バルブ、エイドリアン; チュー、ソンチュン (2020)。「ハミルトニアンとランジュバンモンテカルロ」。モンテカルロ法。シンガポール：シュプリンガー。pp. 281–326。ISBN 978-981-13-2970-8。
• Neal, Radford M (2011)。「ハミルトン力学を用いた MCMC」(PDF)。Steve Brooks、Andrew Gelman、Galin L. Jones、Xiao-Li Meng (編)。マルコフ連鎖モンテカルロハンドブック。Chapman and Hall/ CRC。ISBN 9781420079418。

• ベタンコート、マイケル。「ハミルトンモンテカルロによる効率的なベイズ推論」。MLSSアイスランド 2014 – YouTube経由。
• McElreath、リチャード。「マルコフ連鎖モンテカルロ」。統計的再考 2022 – YouTube経由。
• ハミルトンモンテカルロ法をゼロから学ぶ
• 最適化とモンテカルロ法