力学を使わないハミルトニアン・モンテカルロ法の説明

この記事は「KCSアドベントカレンダー Advent Calendar 2018」の12月9日の記事です。

はじめに

ある確率分布からサンプリングすることは、確率分布の関数が分かっていたとしても一般には簡単なことではありません。幅広い種類の確率分布から効率よくサンプリングを行う手法としてマルコフ連鎖モンテカルロ（MCMC）が知られています。今回のテーマであるハミルトニアン・モンテカルロ法（ハイブリッド・モンテカルロ法）(Duane, Kennedy, Pendleton, & Roweth, 1987) はMCMCの一種です。

機械学習を学んでいる人はベイズ学習で事後分布からサンプリングをするためにMCMCを学んでいる人が多いと思います。ハミルトニアン・モンテカルロ法はMCMCを学習していて躓きやすい分野の一つではないでしょうか。ハミルトニアンという力学の概念を応用したサンプリング方法であるため、解析力学の説明から始まったり、力学を用いて直感的な説明を行なっている文献が多くあります。歴史的に考えても力学を用いた説明が自然だと思うのですが、個人的には力学の概念であると言うことを意識して説明されるとハミルトニアン・モンテカルロの本質的な部分が分かりにくくなるのでは無いかと感じました。

この記事は力学のことを忘れてハミルトニアンの数学的性質のみに注目してハミルトニアン・モンテカルロ法を説明するという試みです。

一般的な力学から始まる説明はNeal (2011) やBishop (2006) など多くの文献で紹介されています。

前提

この記事は全く知識がなくても読めますがメトロポリス・ヘイスティング法を知っていることを前提としています。この辺りは公開されている文献も多いですし理解もしやすいと思います。たとえばcourseraの「Bayesian Methods for Machine Learning」はコンパクトにまとまっています。

エルゴート性と詳細つりあい条件

ハミルトニアン・モンテカルロ法より前の部分は前提とすると言いましたが、重要な部分を確認しておきましょう。この節の内容は証明など細かい部分には触れません。分からない点があった場合には、マルコフ連鎖について復讐をしてから以降の説明を読んだ方が理解しやすいと思います。

MCMCの目標はサンプリングしたい確率分布に対応したマルコフ連鎖を探し、そのマルコフ連鎖からサンプリングすることで目的の確率分布からの標本を得ることです。マルコフ連鎖からサンプリングを行うときには、初期状態を$z^{(0)}$として最初に$p(z^{(1)} \vert z^{(0)})$からサンプリングし、次に$p(z^{(2)} \vert z^{(1)})$からサンプリングをする、ということを繰り返して標本を得ていくことになります。このようにしてサンプリングを行うと一般には前の状態に依存した標本が得られてしまいますが、一旦そこには目を瞑りましょう。サンプリングを行いたい確率分布$p^{*}(z)$に収束する、つまり

$$\forall z^{(0)} \hspace{1em} \lim_{n\rightarrow \infty}p(z^{(n)}=z\vert z^{(0)}) = p^{*}(z)$$

が成り立つようなマルコフ連鎖を見つけることがMCMCの目標になります。

$p^*(z)$の分布を持つ状態から遷移確率$p(z \vert z')$で遷移をすると以下のように新しい分布が得られます。

$$p^{**}(z) = \int p(z \vert z')p^{*}(z') dz'$$

ここで、以下が成り立つような分布$p^{*}(z)$をマルコフ連鎖の定常分布と言います。すなわち、遷移をしても分布が変化しないものを定常分布と言います。

$$p^{*}(z) = \int p(z \vert z')p^{*}(z') dz'$$

一般にマルコフ連鎖の定常分布は一意ではありませんが、マルコフ連鎖がエルゴート的であれば定常分布を一意に持つということが知られています。マルコフ連鎖がエルゴート的であるとは、任意の始点からマルコフ連鎖を開始した時に（始点と終点に依存しない）ある有限回の推移で任意の点（終点）に推移する確率が0ではないということです。

そのうえ、エルゴート的なマルコフ連鎖は定常分布に収束するということが知られています。よって、定常分布として目標の確率分布$p^{*}(z)$を持つようなエルゴート的なマルコフ連鎖からサンプリングを行えば良いということが分かります。

ここで、あるエルゴート的なマルコフ連鎖が目標となる確率分布を定常分布として持っているかを判定する方法が必要になります。ある確率分布$p^{*}(z)$とマルコフ連鎖について以下の関係が成り立っていれば$p^{*}(z)$がこのマルコフ連鎖の定常分布であるということが知られています。この十分条件のことを詳細つりあい条件と言います。

$$p^{*}(z)p(z' \vert z) = p^{*}(z')p(z \vert z')$$

メトロポリス・ヘイスティング法の欠点

当たり前ですが、メトロポリス・ヘイスティング法などのMCMCの手法はエルゴート的で目標の分布を定常分布として持つマルコフ連鎖からサンプリングを行う手法なので、新しい手法は必要ないように思えます。なぜ、わざわざハミルトニアン・モンテカルロ法を考える必要があるのでしょうか。

メトロポリス・ヘイスティング法は正しいとは限らない遷移確率（提案分布という）$p(z\vert z')$からサンプリングをして、適切な割合で棄却する（得られた標本を捨てる）ことで詳細つりあい条件を満たすようにしてサンプリングを行うという手法でした。そのため、メトロポリス・ヘイスティング法において重要となるのが提案分布の選び方です。提案分布の分散が大きすぎる場合は適切ではないサンプリングが行われる確率が高くなるので棄却されるサンプルが多くなってしまいます。提案分布の分散が小さすぎる場合は狭い範囲からばかりサンプリングされてしまい、確率分布の全体からサンプリングするまでに何度もサンプリングを行う必要があります。いずれの場合もサンプリングに時間がかかってしまうことになりますし、一般に適切な提案分布を選ぶことは簡単ではありません。

ハミルトニアン・モンテカルロ法は棄却率が低くしたまま広い範囲からサンプリングを行えるようにした、メトロポリス・ヘイスティング法の改良版です。

メトロポリス・ヘイスティング法の確定的な遷移への一般化

メトロポリス・ヘイスティング法は提案分布$p(z' \vert z)$を用いていましたが、ここで確定的な関数による遷移を考えてみます。たとえば、$n$個の$z$の関数$g_1, \cdots, g_n$の中から$g_i$を確率$p_i$で選び、状態$z$に$g_i$を適用することで遷移するようなマルコフ連鎖を考えることできます。なんのためにこのようなことを考えるのか分からないかもしれないですが、ハミルトニアン・モンテカルロ法を考える上では重要なので先に学んでおきましょう。なお、それぞれの$g_i$について（ほとんど至る所で）逆関数が存在してヤコビアンが0ではないということを仮定します。

このような遷移の場合はエルゴート的であるかが問題となりますが、エルゴート的な遷移と組み合わせれば簡単にエルゴート的な遷移を構成することはできます。たとえば、確定的な遷移のあとにエルゴート的なメトロポリス・ヘイスティング法による遷移を行えば全体としてエルゴート性を保つことができます。

エルゴート性を満たすのであれば、問題となるのは詳細つりあい条件です。まず、確定的な遷移について$p(z' \vert z)$は以下のようになると考えることができます。$\delta$はデルタ関数、$J_{g_i}$は関数$g_i$のヤコビアンです。

$$p(z' \vert z) = \sum_{i=1}^n \frac{p_i \delta(g_i^{-1}(z') - z)}{\vert J_{g_i}\vert}$$

$z' = g(z)$以外の場合については$p(z' \vert z)=0$であり、写像による確率密度の変化がヤコビアンによって表されていると考えると直感的ですが、詳しい説明はAppendixを参照してください。

ここで、$z$から$z'$に遷移する際に$A(z', z)$の確率で標本が採用される場合には$p(z' \vert z)$は以下のように表されます。

$$p(z' \vert z) = A(z', z) \sum_{i=1}^n \frac{p_i \delta(g_i^{-1}(z') - z)}{\vert J_{g_i}\vert}$$

これを詳細つりあい条件に代入すると以下のように表されます。

$$\begin{equation} p^{*}(z) A(z', z) \sum_i \frac{p_i \delta(g_i^{-1}(z') - z)}{\vert J_{g_i}\vert} = p^{*}(z') A(z, z') \sum_i \frac{p_i \delta(g_i^{-1}(z) - z')}{\vert J_{g_i}\vert} \label{functiondetailedbalance} \end{equation}$$

詳細つりあい条件は少し見慣れない形ですが、メトロポリス・ヘイスティング法のアルゴリズム自体はほとんど同じです。選択された関数による$z'$への遷移の後に以下の確率$A(z', z)$で標本を採用します。

$$A(z', z) = \min \left(1, \frac{p^{*}(z') \sum_i \frac{p_i \delta(g_i^{-1}(z) - z')}{\vert J_{g_i}\vert}} {p^{*}(z) \sum_i \frac{p_i \delta(g_i^{-1}(z') - z)}{\vert J_{g_i}\vert}} \right)$$

この採用確率で行なった遷移が詳細つりあい条件を満たすことは以下のように示すことができます。

$$\begin{align*} p^{*}(z) & A(z', z) \sum_i \frac{p_i \delta(g_i^{-1}(z') - z)}{\vert J_{g_i}\vert} \\ &= p^{*}(z) \min \left(1, \frac{p^{*}(z') \sum_i \frac{p_i \delta(g_i^{-1}(z) - z')}{\vert J_{g_i}\vert}} {p^{*}(z) \sum_i \frac{p_i \delta(g_i^{-1}(z') - z)}{\vert J_{g_i}\vert}} \right) \sum_i \frac{p_i \delta(g_i^{-1}(z') - z)}{\vert J_{g_i}\vert} \\ &= \min \left(p^{*}(z) \sum_i \frac{p_i \delta(g_i^{-1}(z') - z)}{\vert J_{g_i}\vert}, p^{*}(z') \sum_i \frac{p_i \delta(g_i^{-1}(z) - z')}{\vert J_{g_i}\vert} \right) \\ &= p^{*}(z') A(z, z') \sum_i \frac{p_i \delta(g_i^{-1}(z) - z')}{\vert J_{g_i}\vert} \end{align*}$$

ハミルトニアン

いよいよハミルトニアン・モンテカルロ法の説明に移ります。目標となる確率分布を$p^{*}(z)$とします。

まず、確率分布を以下のように書き換えます。$Z_z$は正規化定数です。

$$p^{*}(z) = \frac{1}{Z_z} \exp (-E(z))$$

ハミルトニアン・モンテカルロ法では、新しい変数$r$との同時確率$p^{*}(z, r)$からサンプリングを行います。一般に$z$は多次元なので、$r$は$z$と同じ次元を持つ変数とします。$p^{*}(z, r)$からサンプリングをすることができれば、$r$について周辺化してしまえば$p^{*}(z)$からサンプリングができたことになりますね。

そもそも$p^{*}(z)$からのサンプリングが難しいという前提があるのに$p^{*}(z, r)$からのサンプリングなんて出来ないような気がしますが、良い性質の$p^{*}(z, r)$を構成して$p^{*}(z, r)$からサンプリングをするというのがハミルトニアン・モンテカルロ法のポイントとなります。

$r$は勝手に導入した変数なので、台が実数全体でサンプリングが簡単にできればどのような確率分布を持っていても構いません。例えば、標準正規分布などがよく使われます。$r$の確率分布も$p(z)$と同じように書き換えておきましょう。

$$p_r(r) = \frac{1}{Z_r} \exp (-K(r))$$

そして、$z$と$r$は独立であると仮定します。このとき同時分布$p^{*}(z, r)$は以下のように与えられます。

$$p^{*}(z, r) = p^{*}(z)p_r(r) = \frac{1}{Z_H} \exp(-(E(z) + K(r))) = \frac{1}{Z_H} \exp(-H(z, r))$$

ここで出てきた$H(z, r) = E(z) + K(r)$がハミルトニアンと呼ばれる形式です。ようやくハミルトニアンが出てきましたが、このままでは何が嬉しいのか全く分かりません。次の節で導かれるハミルトニアンの性質を見てもらえば、わざわざ変数$r$を導入してハミルトニアンで記述される同時分布を考えている理由が分かってきます。

等ハミルトニアン面上の移動

理想的な遷移

$p^{*}(z, r)$から$(z, r)$をサンプリングすることが目標となるので、まずはメトロポリス・ヘイスティング法でサンプリングをすることを考えましょう。恣意的ですが、確定的な遷移で提案分布が対称の場合を考えます。エルゴート性は前述の通り後からなんとかなるので一旦置いておきます。最も簡単な例として各$(z, r)$ごとに二つの関数$g_1$、$g_2$が等確率で選ばれ、$g^{-1}_1 = g_2$を満たすようなマルコフ連鎖が考えられます。

このような場合には式(\ref{functiondetailedbalance})の詳細つりあい条件を以下のように整理することができます。仮定より、$\delta(g_2^{-1}(z', r') - (z, r)) = 1$、$\delta(g^{-1}(z, r) - (z', r')) = 1$、および$\vert J_{g_1} \vert = 1 / \vert J_{g_2} \vert$であることに注意してください。

$$\begin{align*} p^{*}(z, r) A(z', z) \sum_{i=1}^2 \frac{\frac{1}{2} \delta(g^{-1}_i(z', r') - (z, r))}{\vert J_{g_i} \vert} &= p^{*}(z', r') A(z, z') \sum_{i=1}^2 \frac{\frac{1}{2} \delta(g^{-1}_i(z, r) - (z', r'))}{\vert J_{g_i} \vert} \\ \Leftrightarrow \hspace{1em} p^{*}(z, r) A(z', z) \frac{1}{\vert J_{g_2}\vert} &= p^{*}(z', r') A(z, z') \frac{1}{\vert J_{g_1} \vert} \\ \Leftrightarrow \hspace{1em} p^{*}(z, r) A(z', z) &= p^{*}(z', r') A(z, z') \frac{1}{\vert J_{g_1} \vert^2} \end{align*}$$

このとき、詳細つりあい条件が成り立つような採用確率$A(z', z)$は以下のように表されます。

$$\begin{equation} A(z', z) = \min \left(1, \frac{p^{*}(z', r')}{p^{*}(z, r)\vert J_{g_1} \vert^2} \right) \label{acceptancerate} \end{equation}$$

大きく移動する遷移を行ったときに棄却率が高くなってしまうのがメトロポリス・ ヘイスティング法の欠点でした。これを解決する方法としては、この例の場合は明らかな十分条件として$p^{*}(z', r')=p^{*}(z, r)$かつ$\vert J_{g_1} \vert=1$が成り立てば棄却率が0になるということが分かります。エルゴート性については後で考えないといけませんが、この条件を満たす遷移であれば詳細つりあい条件を満たしている上に、$(z, r)$から大きく離れた$(z', r')$に移動したとしても全く棄却されることはありません。

このような都合の良い遷移が簡単に見つかるようには思えませんが、ハミルトニアンについてはこれを満たす遷移が知られています。

等ハミルトニアン面上の移動

再び天下り式の説明になってしまいますが、新たな変数$t$（時間と考えると直感的に理解しやすいですね）を導入し、以下の性質を満たすような遷移を考えてみます。

$$\begin{equation} \begin{split} \frac{d z_i}{dt} &= \frac{\partial H}{\partial r_i} \\ \frac{d r_i}{dt} &= - \frac{\partial H}{\partial z_i} \end{split} \label{dynamics} \end{equation}$$

適当に$T$を決めて、$1/2$の確率で$T$もしくは$-T$の大きさだけ変数$t$について上の関係式の微分に沿って移動するという遷移を、先ほど考えていた決定的な遷移$g_1$、$g_2$として定義することができます。この遷移が$p^{*}(z', r')=p^{*}(z, r)$および$\vert J_{g_1} \vert=1$を満たせば、この遷移が理想的な遷移であるということになります。

$p^{*}(z', r')=p^{*}(z, r)$は簡単に示すことができます。$p(z, r) = \frac{1}{Z_H} \exp(-H(z, r))$だったので、変数$t$の変化に対してハミルトニアンが保存されることを示せば十分です。これは以下のように示されます。

$$\begin{align*} \frac{dH}{dt} &= \sum_i \left\{ \frac{\partial H}{\partial z_i} \frac{d z_i}{dt} + \frac{\partial H}{\partial r_i} \frac{d r_i}{dt} \right\} \\ &= \sum_i \left\{ \frac{\partial H}{\partial z_i} \frac{\partial H}{\partial r_i} - \frac{\partial H}{\partial r_i} \frac{\partial H}{\partial z_i} \right\} \\ &= 0 \end{align*}$$

ヤコビアンが1であることの証明は難しいのでここでは省略します。この事実はLiouvilleの定理として知られています。ReferencesにLiouvilleの定理を説明した記事があるので、きちんと理解したいという人は読んでみてください。とはいえ、下記のように実際には厳密に式(\ref{dynamics})を満たす遷移を計算することはできないので、式(\ref{dynamics})による遷移のヤコビアンが1であるということは実用上はそれほど重要ではありません。むしろヤコビアンが1になるように遷移を近似するということが重要です。

考えなくてはならないこと

以上より式(\ref{dynamics})を満たすような遷移は、全く棄却することなく詳細つりあい条件を満たすという理想的な遷移になっています。これがハミルトニアン・モンテカルロ法における本質的な部分なので、ここまで理解できたらハミルトニアン・モンテカルロ法は概ね理解したといっても過言ではありません。

しかし、この遷移には大きく2つの問題点があります。1つ目は式(\ref{dynamics})を満たすような遷移を実際に行うことは難しいという点です。一般にコンピュータで厳密な積分を行うことはできません。2つ目はエルゴート性を満たしていないということです。式(\ref{dynamics})を満たす遷移はハミルトニアンを保存する、つまり確率密度が同じ場所のみを遷移するのでした。もちろん、一般には確率密度が異なる部分が存在するので、この遷移だけでは全ての状態を訪れることはできません。

以降の節で、この二つの問題点を解決する方法を紹介します。

積分の近似

ヤコビアンが1の近似

式(\ref{dynamics})を満たすような遷移を厳密に行うことはできないので、近似を考えます。最も簡単な方法として、微小な$\epsilon$を用いて以下のように遷移することが考えられます（オイラー法）。

$$\begin{align*} z_i(t + \epsilon) &= z_i(t) + \epsilon \frac{\partial H}{\partial r_i} (z(t), r(t)) = z_i(t) + \epsilon \frac{\partial K}{\partial r_i} (r(t)) \\ r_i(t + \epsilon) &= r_i(t) - \epsilon \frac{\partial H}{\partial z_i} (z(t), r(t)) = r_i(t) - \epsilon \frac{\partial E}{\partial z_i} (z(t)) \end{align*}$$

しかしオイラー法による遷移は誤差が大きいうえにヤコビアンの値が1になるとは限りません。$p^{*}(z', r')$の値について多少の誤差は仕方ないですが、ヤコビアンの値が1でなくなると詳細つりあい条件が複雑になってしまいます。そこで、ハミルトニアン・モンテカルロ法では誤差が小さいうえにヤコビアンの値が1になるような近似が用いられます（リープ・フロッグ法）。

$$\begin{align*} r_i(t + \epsilon/2) &= r_i(t) - \frac{\epsilon}{2} \frac{\partial E}{\partial z_i} (z(t)) \\ z_i(t + \epsilon) &= z_i(t) + \epsilon \frac{\partial K}{\partial r_i} r_i(t + \epsilon/2) \\ r_i(t + \epsilon) &= r_i(t + \epsilon/2) - \frac{\epsilon}{2} \frac{\partial E}{\partial z_i} (z(t+\epsilon)) \end{align*}$$

オイラー法と大きく異なるようには見えませんが、それぞれの遷移は一部の変数（一番上では$r_i$）のみが変化し、変化量は固定された変数（一番上では$z$）のみに依存するということが重要な性質です。このような写像はせん断写像と呼ばれ、ヤコビアンが1であることが知られています。

これは簡単に示すことができるので、一番上の遷移について示しましょう。この遷移を関数$(g_1, \ldots, g_n, l_1, \ldots, l_n)$と考えます。つまり$g_k$と$l_k$ $(k = 1, \ldots, n)$は以下のような関数になっています。

$$\begin{align*} z_k(t + \epsilon/2) &= g_k(z(t), r(t)) \\ r_k(t + \epsilon/2) &= l_k(z(t), r(t)) \end{align*}$$

このとき、$l_i$のみが$l_i(z, r) = r_i - \frac{\epsilon}{2} \frac{\partial E}{\partial z_i} (z)$となり、他の関数は$g_k(z, r) = z_k$、$l_k(z, r) = r_k$となります。よって、この関数のヤコビ行列の対角成分は全て1で、それ以外の成分は下三角成分の$n+i$行目（$l_i$の$z_k$による偏微分）だけが非零になる可能性があります。ここで$l_i$が$r_i$と$z$の関数であることを使っています。このような行列の行列式は1なので、ヤコビアンが1であることが示されました。他の遷移についても同様に示すことができます。

棄却率

リープ・フロッグ法による遷移はヤコビアンが1になりますが、ハミルトニアンが保存されるとは限りません。このような場合には何もしないと詳細つりあい条件を満たさないことになります。式(\ref{acceptancerate})より、リープ・フロッグ法による遷移で得られたサンプルを以下の確率$A(z', z)$で採用すれば詳細つりあい条件を満たすということが分かります。

$$A(z', z) = \min \left( 1, exp(H(z, r)-H(z', r')) \right)$$

ハミルトニアン・モンテカルロ法は棄却率を小さくしたいという動機があったので、結局は棄却しなくてはならないことは残念ですが、普通のメトロポリス・ヘイスティング法とは全く状況が異なるということに注目してください。リープ・フロッグ法はハミルトニアンが保存される遷移を近似しているので、厳密に同じではないにせよ$H(z, r) \simeq H(z', r')$ということが期待されます。つまり、棄却率は常に非常に低くなるということです。リープ・フロッグ法を何度か行うという遷移をすることで、現在の状態$(z, r)$から大きく離れた状態に低い棄却率で移動することができます。

エルゴート性を満たすために

ハミルトニアンの微分に沿った移動は詳細つりあい条件を満たしていますが、一般にエルゴート性は満たしていません。メトロポリス・ヘイスティング法のように任意の値を取りうる提案分布を使ったサンプリングを追加すればエルゴート性は満たされますが、せっかく棄却率を下げた意味がなくなってしまいます。そこで、リープ・フロッグ法による遷移の後に$r$のみを真の分布からサンプリングするということを考えます。$r$の分布は自分で定義したものなので、サンプリングが簡単にできるような確率分布を持つように定義していました。真の分布からのサンプリングなので、得られたサンプルを棄却する必要なく詳細つりあい条件を満たしていることは明らかです。実際に、以下のように簡単に示されます。ここで、$z=z'$であることと$z$と$r$が独立であることに注意してください。

$$\begin{align*} p^{*}(z, r)p(z', r' \vert z, r) &= p(z)p(r)p(r') \\ &= p(z')p(r')p(r) \\ &= p^{*}(z', r')p(z, r \vert z', r') \end{align*}$$

実はこれでもエルゴート性を満たさない場合もあるのですが、それほど難しい話ではないので気になる人はNeal (2011) やBishop (2006) などを参照してください。

まとめ

ハミルトニアン・モンテカルロ法の1回の遷移は

1. リープ・フロッグ法でハミルトニアンが大きく変わらない遷移を行う
2. 新しく導入した変数$r$の真の分布からサンプリングを行う

という二段階で構成されます。ポイントとなるのはリープ・フロッグ法による遷移で、これによって低い棄却率を保ったまま現在の状態とは大きく異なる状態に遷移することができます。従来のメトロポリス・ヘイスティング法では、棄却率を小さくしようとすると大きく移動することができず、大きく移動すると棄却率が大きくなってしまうという問題があったので、ハミルトニアン・モンテカルロ法はこの大きな問題を解決したことになります。

さて、最後まで力学を考えずに説明することができました。この記事でハミルトニアン・モンテカルロ法を知った人は「どうして都合の良い遷移を思いつくことができたのか」という疑問を持ったかもしれないですが、これは解析力学でよく知られているハミルトニアンの性質から考えられたものです。気になる人はぜひ調べてもらうと良いと思いますが、知らなくてもハミルトニアン・モンテカルロ法の目的や手順を理解することができたのではないかと思います。

References

1. Duane, S., Kennedy, A. D., Pendleton, B. J., & Roweth, D. (1987). HYBRID MONTE CARLO. PHYSICS LETTERS B, 195(2).
2. Neal, R. M. (2011). MCMC using Hamiltonian dynamics. In Handbook of Markov Chain Monte Carlo.
3. Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer.

• 11. Introduction to Liouville’s Theorem
• Preservation of phase space volume and Liouville’s theorem

Appendix

確定的関数による遷移についての遷移確率と詳細つりあい条件

遷移確率

まずは一つの関数$g$による遷移について考えます。これについて遷移確率$p(z' \vert z)$は以下のように考えることができます。説明の際に記号を分かりやすくするためにヤコビアンを書き直しています。

$$\begin{align*} p(z' \vert z) = \frac{\delta(g^{-1}(z') - z)}{\vert J_g \vert} = \delta(g^{-1}(z') - z) \left\vert \det \frac{\partial g^{-1}(z')}{\partial x} \right\vert \end{align*}$$

本文の通りこれは直感的な定義ですが、これが$p(z' \vert z)$の定義として妥当だということを説明します。提案分布を用いた場合には元の確率密度関数$p^*$と遷移後の確率密度関数$p^{**}$に以下のような関係が成り立っていました。

$$p^{**}(z') = \int p^*(z)p(z' \vert z) dz$$

ここで、関数$g$による遷移というのは関数$g$による変数変換に他なりません。よって、変換後の確率密度は以下のように表されます。

$$p^{**}(z') = p^*(g^{-1}(z')) \left\vert \det \frac{\partial g^{-1}(z')}{\partial x} \right\vert$$

これは先ほど定義した遷移確率$p(z' \vert z)$を用いて表すことができ、以下のように提案分布を用いた遷移の形式と一致していることがわかります。どうやら確定的な関数による遷移についての$p(z' \vert z)$の定義は妥当なようです。

$$\begin{align*} p^{**}(z') &= p^*(g^{-1}(z)) \left\vert \det \frac{\partial g^{-1}(z')}{\partial x} \right\vert \\ &= \int p^*(z) \delta(g^{-1}(z') - z) \left\vert \det \frac{\partial g^{-1}(z')}{\partial x} \right\vert dz \\ &= \int p^*(z)p(z' \vert z) dz \end{align*}$$

詳細つりあい条件

遷移確率の定義について、詳細つりあい条件を考えた場合にも提案分布を用いた場合の一般化になっていることを確認します。上記の遷移確率を代入すると、サンプルの棄却を行わない場合の詳細つりあい条件は以下のように表されます。

$$p^*(z) \delta(g^{-1}(z') - z) \left\vert \det \frac{\partial g^{-1}(z')}{\partial x} \right\vert = p^*(z') \delta(g^{-1}(z) - z') \left\vert \det \frac{\partial g^{-1}(z)}{\partial x} \right\vert$$

これが成り立っているときに関数$g$が定常分布に準ずる遷移となっていること、つまり以下の関係が成り立っていることを示します。

$$p^{*}(z') = \int p^*(z)p(z' \vert z) dz$$

これは以下のように示されます。

$$\begin{align*} \int p^*(z)p(z' \vert z) dz &= p^*(g^{-1}(z)) \left\vert \det \frac{\partial g^{-1}(z')}{\partial x} \right\vert \\ &= \int p^*(z) \delta(g^{-1}(z') - z) \left\vert \det \frac{\partial g^{-1}(z')}{\partial x} \right\vert dz \\ &= \int p^*(z') \delta(g^{-1}(z) - z') \left\vert \det \frac{\partial g^{-1}(z)}{\partial x} \right\vert dz \\ &= p^*(z') \int \delta(g^{-1}(z) - z') \left\vert \det \frac{\partial g^{-1}(z)}{\partial x} \right\vert dz \\ &= p^*(z') \end{align*}$$

ここで、最後の等式についてはデルタ関数の以下の性質を利用しています。

$$\begin{align*} \delta(f(x)) = \left\{ \begin{array}{ll} \frac{1}{\left\vert \det \frac{\partial f(a)}{\partial x} \right\vert} & (f(a) = 0) \\ 0 & (otherwise) \end{array} \right. \end{align*}$$

よって、詳細つりあい条件が成り立っていれば関数$g$による遷移によって分布が変化しないということが分かりました。

以上のことが分かれば、棄却を行う場合や$n$個の関数$g_1, \cdots, g_n$の中から$g_i$を確率$p_i$で選ぶような場合についての一般化についても同様に考えることができます。