はじめてのMCMC (ハイブリッド・モンテカルロ)
こんにちはtatsyです。
はじめてのMCMCも早くも3記事目ですが、自分的にはもう少し書きたいことがある感じです。というわけで今回はハイブリッド・モンテカルロです。
ちなみにこれまでの記事はこちらです。
ハイブリッド・モンテカルロとは？
ハイブリッド・モンテカルロは今までのM-H法やギブス・サンプリングと違い、サンプリングを行いたい確率変数$x$の他に拡張変数$u$を導入して、これらの変数のペア$(x, u)$サンプリングします。
ハイブリッド・モンテカルロは何か特定のやり方を指すというよりは、拡張変数を導入するMCMC全般を指すようです。今回の記事では、ハイブリッド・モンテカルロの中でも、経路積分を用いてサンプル変数を動かすハミルトニアン・モンテカルロ(HMC)というアルゴリズムを紹介します。
ハミルトニアン・モンテカルロ (HMC)
HMCは、これまでに見てきたM-H法やギブス・サンプリングのように提案分布によって提案された確率変数をそのまま採択・棄却するのではなく、「ハミルトン力学 (Hamiltonian dynamics)」の考え方に従って、より確率密度の高い点へと提案サンプルを移動します。
ハミルトン力学の言葉で置き換えると$x$は位置を表す変数 (position variable)であり、$u$は運動量を表す変数 (momentum variable)と呼ばれます。これらの変数に対し、ハミルトン力学ではポテンシャルエネルギー$U(x)$と力学的エネルギー$K(u)$の和としてハミルトニアン: $$ H(x_i, u_i) = U(x_i) + K(u_i) $$ を考えます。
この際、$x$が持つポテンシャルエネルギーはサンプルしたい確率密度関数(=不変分布)を用いて $$ U(x_i) = -\log p(x) $$ と定義し、力学的エネルギーについては、質量テンソル$M$の逆数を考えて、 $$ K(u) = \frac{1}{2} u^{T} M^{-1} u $$ のように定義します。
一方、ハミルトン力学においては、時間$t$と、上記の位置・運動量変数が以下のハミルトン方程式を満たします (ニュートン力学における運動方程式のようなもの)。 $$ \begin{aligned} \frac{\partial x}{\partial t} &= \frac{\partial H}{\partial u} \\ \frac{\partial u}{\partial t} &= -\frac{\partial H}{\partial x} \end{aligned} $$
このハミルトン方程式にポテンシャル・運動的エネルギーの定義を代入すると、 $$ \begin{aligned} \frac{\partial x}{\partial t} &= M^{-1} u \\ \frac{\partial u}{\partial t} &= \frac{\partial \log p(x)}{\partial x} \end{aligned} $$ という$(x, u)$の更新式が得られます。
この更新式に従って、初期の提案サンプル$(x_0, u_0)$を一定回数(ここでは$L$回とする)したものを$(x_L, u_L)$として、この値を元にした以下の採択確率を考えます。 $$ a(x_L, u_L | x_0, u_0) = \min \left[ 1, \frac{\exp H(x_L, u_L)}{\exp H(x_0, u_0)} \right] $$ なお、ハミルトン方程式に基づく$(x_l, u_l)$の更新には、オイラー法やLeapfrog積分などの数値解法を用います。
簡単のため、以下に示す実装では、質量テンソルが単位行列であると仮定します。今回は少し複雑なので、コードに行く前に簡単に流れをまとめておきます。
- 現在の点から運動量として$u_0$を平均0、分散1の正規分布からサンプリングする。
- Leapfrog積分: ハミルトン方程式に従い$x_0 = x^t$ならびに$u_0$を適当なステップ幅$\Delta$で$L$回更新する。
- $x_{l+\frac{1}{2}} = x_l + \frac{\Delta}{2} u_l$
- $u_{l+1} = u_l + \Delta \frac{\partial \log p}{\partial x}(x_{l+\frac{1}{2}})$
- $x_{l+1} = x_{l + \frac{1}{2}} + \frac{\Delta}{2} u_{l+1}$
- 採択確率$a(x_L, u_L | x_0, u_0)$に基づいて$(x_L, u_L)$を採択するかを決定する。
上記のコードでは、M-H法のサンプルと揃える目的で対数確率の微分を数値的に求めていますが、これが解析的に求まり、なおかつ簡単に評価できる場合には、より効率的に計算が可能になるはずです。
まとめ
ハイブリッド・モンテカルロはナイーブなM-H法とくらべ、一回の更新にかかる計算量はLeapfrog積分の分だけ増加します。その一方で、より高い確率で未知サンプルを得ることができるため、より複雑でサンプルが困難な確率密度分布に対しては有効な手法です (上記のサンプルだと効果が見えづらいですが…)。
次の記事ではスライス・サンプリングを紹介いたしますので、よろしければそちらもご覧ください。
最後までお読みいただきありがとうございました。