ハミルトンモンテカルロ

2次元確率分布のハミルトンモンテカルロサンプリング

ハミルトニアン・モンテカルロ法(元々はハイブリッド・モンテカルロとして知られていた)は、マルコフ連鎖モンテカルロ法の一種であり、直接サンプリングすることが困難な目標確率分布に収束する分布を持つランダムサンプルの系列を求める。この系列は、期待値モーメントといった目標分布の積分値を推定するために用いられる。

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

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

アルゴリズム

サンプリングする対象分布が( )であり、サンプルのチェーンが必要であるとします。 f×{\displaystyle f(\mathbf {x} )}×Rd{\displaystyle \mathbf {x} \in \mathbb {R} ^{d}}d1{\displaystyle d\geq 1}X0X1X2{\displaystyle \mathbf {X} _{0},\mathbf {X} _{1},\mathbf {X} _{2},\ldots }

ハミルトン方程式

d×dtHpそしてdpdtH×{\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}}p{\displaystyle p_{i}}{\displaystyle i}H{\displaystyle H}M{\displaystyle M}

H×pあなた×+12pTM1p{\displaystyle H(\mathbf {x} ,\mathbf {p} )=U(\mathbf {x} )+{\dfrac {1}{2}}\mathbf {p} ^{\text{T}}M^{-1}\mathbf {p} }

ここでは位置エネルギーである。標的の位置エネルギーは次のように与えられる。 あなた×{\displaystyle U(\mathbf {x} )}

あなた×lnf×{\displaystyle U(\mathbf {x} )=-\ln f(\mathbf {x} )}

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

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

pnt+Δt2pntΔt2あなた×|××nt{\displaystyle \mathbf {p} _{n}\left(t+{\dfrac {\Delta t}{2}}\right)=\mathbf {p} _{n}(t)-{\dfrac {\Delta t}{2}}\nabla \left.U(\mathbf {x} )\right|_{\mathbf {x} =\mathbf {x} _{n}(t)}}
×nt+Δt×nt+ΔtM1pnt+Δt2{\displaystyle \mathbf {x} _{n}(t+\Delta t)=\mathbf {x} _{n}(t)+\Delta tM^{-1}\mathbf {p} _{n}\left(t+{\dfrac {\Delta t}{2}}\right)}
pn(t+Δt)=pn(t+Δt2)Δt2U(x)|x=xn(t+Δt){\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+\Delta t)}}

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

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

系を熱力学的平衡分布へと導く一つの方法は、メトロポリス・ヘイスティングス・アルゴリズムを用いて粒子の状態を変化させることです。まずリープフロッグ・ステップを適用し、次にメトロポリス・ヘイスティングス・ステップを適用します。

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

Xn+1|Xn=xn={xn(LΔt)with probability α(xn(0),xn(LΔt))xn(0)otherwise{\displaystyle \mathbf {X} _{n+1}|\mathbf {X} _{n}=\mathbf {x} _{n}={\begin{cases}\mathbf {x} _{n}(L\Delta t)&{\text{with probability }}\alpha \left(\mathbf {x} _{n}(0),\mathbf {x} _{n}(L\Delta t)\right)\\\mathbf {x} _{n}(0)&{\text{otherwise}}\end{cases}}}

どこ

α(xn(0),xn(LΔt))=min(1,exp[H(xn(LΔt),pn(LΔt))]exp[H(xn(0),pn(0))]).{\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).}

完全な更新は、まず運動量をランダムにサンプリングし(以前の反復とは独立して)、次に運動方程式を積分し(例えばリープフロッグ法を用いて)、最後にメトロポリス・ヘイスティングスの受理/拒否ステップから新しい構成を得るという手順で行われます。この更新メカニズムを繰り返して を取得します。 p{\displaystyle \mathbf {p} }Xn+1,Xn+2,Xn+3,{\displaystyle \mathbf {X} _{n+1},\mathbf {X} _{n+2},\mathbf {X} _{n+3},\ldots }

Uターンサンプラーなし

No U-Turn Sampler (NUTS) [ 5 ]は、ステップ数を自動的に制御する拡張です。チューニングは非常に重要です。例えば、1次元の場合、ポテンシャルは であり、これは単振動子のポテンシャルに相当します。値が大きすぎると粒子が振動し、計算時間を無駄にします。値が小さすぎると、粒子はランダムウォークのように振舞います。 L{\displaystyle L}L{\displaystyle L}N(0,1/k){\displaystyle {\text{N}}(0,1/{\sqrt {k}})}U(x)=kx2/2{\displaystyle U(x)=kx^{2}/2}L{\displaystyle L}L{\displaystyle L}

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

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

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

(xn+xn)pn<0or.(xn+xn)pn+<0{\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}

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

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

または

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

ここで、たとえば、。 δ=1000{\displaystyle \delta =1000}

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

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

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

参照

参考文献

  1. ^デュアン, サイモン; ケネディ, アンソニー D.; ペンドルトン, ブライアン J.; ロウエト, ダンカン (1987). 「ハイブリッドモンテカルロ」. Physics Letters B. 195 ( 2): 216– 222. Bibcode : 1987PhLB..195..216D . doi : 10.1016/0370-2693(87)91197-X .
  2. ^ Neal, Radford M. (1996). 「モンテカルロ法の実装」.ニューラルネットワークのためのベイズ学習. 統計学講義ノート. 第118巻. Springer. pp.  55– 98. doi : 10.1007/978-1-4612-0745-0_3 . ISBN 0-387-94724-8
  3. ^ Gelman, Andrew; Lee, Daniel; Guo, Jiqiang (2015). 「Stan: ベイズ推論と最適化のための確率的プログラミング言語」.教育行動統計ジャーナル. 40 (5): 530– 543. doi : 10.3102/1076998615606113 . S2CID 18351694 . 
  4. ^ Betancourt, Michael (2018-07-15). 「ハミルトニアン・モンテカルロの概念的入門」. arXiv : 1701.02434 [ stat.ME ].
  5. ^ Hoffman, Matthew D; Gelman, Andrew (2014). 「No-U-turnサンプラー:ハミルトンモンテカルロ法における適応的なパス長設定」 . Journal of Machine Learning Research . 15 (1): 1593– 1623. 2024年3月28日閲覧

さらに読む

  • ベタンコート、マイケル、ジロラミ、マーク (2015).「階層モデルのためのハミルトンモンテカルロ法」。サティアンシュ・クマール・ウパディヤイ他編『ベイズ法の最新動向とその応用』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
  • ニール、ラドフォード・M (2011). 「ハミルトン力学を用いたMCMC」(PDF) . スティーブ・ブルックス、アンドリュー・ゲルマン、ガリン・L・ジョーンズ、シャオ・リー・メン編.マルコフ連鎖モンテカルロハンドブック. チャップマン・アンド・ホール/CRC. ISBN 9781420079418