『ベイズ深層学習』勉強会資料 4章の1

記事の概要

『ベイズ深層学習』の輪読会で使用する資料です。
記事では『ベイズ深層学習』を持っていることを前提として、一部の数式や説明を省略しています。同書を持っている人が自習する時の参考になることを目的としています。

記事中の誤りは全て私の理解不足に起因し、輪読会参加者および本の著者、引用記事の作成者とは無関係になります。
誤りや不明点がございましたら、コメント欄などでご指摘いただけると助かります。

サンプリング法

分布$p(\mathbf{z})$に対する関数$f(\mathbf{z})$の期待値$\int f(\mathbf{z})p(\mathbf{z})d\mathbf{z}$を積分計算するのは困難なので、近似計算をする必要がある。
そのための手段の１つがサンプリング法である。

単純モンテカルロ法

分布$p(\mathbf{z})$から任意にT個のサンプル$(\mathbf{z}^{1}, \cdots \mathbf{z}^{T}) \sim p(\mathbf{z})$を選ぶことで、期待値を関数の和に近似できる。

$$\int f(\mathbf{z})p(\mathbf{z})d\mathbf{z} \approx \frac{1}{T} \sum_{t=1}^T f(\mathbf{z}^{(t)}) \equiv \hat{f}$$

この近似式$\hat{f}$の期待値は、関数$f(\mathbf{z})$の期待値に等しい。よって、この近似式は正しい平均を持つと言える。

\begin{eqnarray}
\frac{1}{T} \sum_{t=1}^T \int f(\mathbf{z}^{(t)})p(\mathbf{z}^{(t)})d\mathbf{z}^{(t)}
&=& \frac{1}{T} \sum_{t=1}^T \mathbb{E}[f] = \mathbb{E}[f]
\end{eqnarray}

また、$\hat{f}$の分散も関数$f(\mathbf{z})$の分布$p(\mathbf{z})$に対する分散と等しくなる。
実際に以下に計算してみる。

\begin{eqnarray}
\mathrm{var}[\hat{f}] &=& \mathbb{E} \Bigl[ \Bigl( \hat{f} - \mathbb{E}[ \hat{f}] \Bigl)^2 \Bigl] \\
&=& \mathbb{E}[\hat{f}^2] -  \mathbb{E}[ \hat{f}]^2
\end{eqnarray}

第1項は以下のように展開できる。

\begin{eqnarray}
\mathbb{E}[\hat{f}^2]  &=& 
\mathbb{E} \Bigl[ \frac{1}{T} \sum_{n=1}^T f(\mathbf{z}^{(n)}) \frac{1}{T} \sum_{m=1}^T f(\mathbf{z}^{(m)})  \Bigl] \\

&=& \frac{1}{T^2} \sum_{n=1}^T \sum_{m=1}^T \{ \mathbb{E}[f^2] + \delta_{nm} \mathrm{var}[f]  \} 
\end{eqnarray}

元の式に代入して以下を得る。

\begin{eqnarray}
\mathrm{var}[\hat{f}] 
&=& \frac{1}{T^2} \sum_{n=1}^T \sum_{m=1}^T \{ \mathbb{E}[f^2] + \delta_{nm} \mathrm{var}[f]  \}  -  \mathbb{E}[ \hat{f}]^2 \\

&=&　\frac{1}{T} \mathrm{var}[f] \\

&=&　\frac{1}{T} \mathbb{E} \bigl[  (f - \mathbb{E}[f] )^2 \bigl]

\end{eqnarray}

これより、一部のサンプルをサンプリングしたものから、原理的には十分な精度の期待値を得られることが分かる。

欠点

計算効率が悪いため、この手法は実際に用いられることは少ない。

例えば、本書では式(4.2)において、パラメータ$\theta$と事前分布$p(\theta)$を持つN個の観測データ（事後分布）の確率$p(\mathbf{X})$を単純モンテカルロ法で近似している。

現実の問題では、この近似式は特定のパラメータ$\theta$の範囲でしか大きな値を取らないことが分かっている。
つまりほとんどの$\theta$について、観測データと一致しないので棄却しなくてはならず、サンプルとして受容できる率は低い。
よって受容できるサンプルを十分な数だけ得るには、広範囲な$p(\theta)$を取らないといけなくなり、計算効率が悪くなる。

棄却サンプリング

分布$p(\mathbf{z})$からのサンプリングが難しい場合に、サンプル$(\mathbf{z}^{1}, \cdots \mathbf{z}^{T}) \sim p(\mathbf{z})$を取得する手法。

• 仮定

正規化されていない$\tilde{p}(\mathbf{z})$が分かっている。
$$p(\mathbf{z}) = \frac{1}{Z_p}\tilde{p}(\mathbf{z})$$

• 手法

1. 仮の分布として提案分布$q(\mathbf{z})$を作成する
2. 定数kを任意の$\mathbf{z}$に対して$kq(\mathbf{z}) \geqslant \tilde{p}(\mathbf{z})$が成立するように定める
3. サンプル$\mathbf{z}^{(t)}$を提案分布$q(\mathbf{z})$から取得する
4. サンプル$\tilde{u}$を区間$[0,kq(\mathbf{z}^{(t)})]$の一様分布$\mathrm{Uni}(u|0,kq(\mathbf{z}^{(t)}))$から取得する
5. $\tilde{u} > \tilde{p}(\mathbf{z}^{(t)})$ならばサンプル$\mathbf{z}^{(t)}$を棄却する。$\tilde{u} \leq \tilde{p}(\mathbf{z}^{(t)})$ならばサンプル$\mathbf{z}^{(t)}$を受容する

例：2次式からのサンプリング

簡単な例として、opabinia2さんが機械学習に詳しくなりたいブログ　棄却サンプリングにおいて$p(z)=\sin^2 0.5 \pi z$からのサンプリングを解説されている。
先のアルゴリズムに従い、以下のことを行っている。

1. 提案分布として区間$[-1, 1]$の一様分布$q(z)=0.5$を採用
2. 定数kを2に設定。q(z)の最大値は1なので$kq(z)=1$は常に$1 \geqslant \sin^2 0.5 \pi z$を満たす。
3. 提案分布からサンプリングする
4. 区間$[0,kq(z)]$の一様分布からサンプリングする
5. 受容するかどうかを判定

紹介したサイトでは、提案分布から乱数を生成し、また区間$[0,kq(z)]$の一様分布から乱数を生成し、両者の大小を比較して受容するかどうかを判定するプログラムを作成している。
そして、受容されたサンプルの集合の棒グラフは、$p(z)=\sin^2 0.5 \pi z$の実線グラフとよく一致しており、正しくサンプリングできていることが分かる。

例：ガンマ分布からのサンプリング

PRMLのP243では、棄却サンプリング法を用いて、ガンマ分布

$$
p(z) = \mathrm{Gam}(z|a,b) = \frac{b^a z^{a-1} \exp(-bz)}{\Gamma(a)}
$$
からサンプリングしている。
提案分布にはコーシー分布を用いる。
$$
\frac{1}{\pi} \frac{1}{1+y^2}
$$
ここでyを区間$[0,1]$上の一様分布として、$z=b\tan y + c$でコーシー分布の変数変換を行ったものを$q(z)$とする。

\begin{eqnarray}
q(z) &=& p(y) \Bigl| \frac{dy}{dz} \Bigl|
\end{eqnarray}

ヤコビアンは以下のようにに計算できる。

$$y = \tan^{-1} \bigl( \frac{z-c}{b} \bigl)$$

より$u \equiv z-c/b$として

\begin{eqnarray}
\frac{dy}{dz}  &=& \frac{d}{du} tan^{-1}(u) \frac{du}{dz} \\
&=& \frac{1}{1+u^2} \frac{du}{dz} \\
&=& \frac{1}{1+(z-c)^2/b^2}\frac{1}{b}
\end{eqnarray}

これを元の式に代入して以下を得る。$c=a-1$、$b^2=2a-1$に選んだものを提案分布q(z)とする

$$
q(z) = \frac{k}{1+(z-c)^2/b^2}
$$

具体的にはysekkyさんがサンプリング手法について(自分なりに)わかりやすく書くよ！で計算されている。

1. $a=3$、$b=1.0$に選んだものを提案分布q(z)とする
2. $k$をなるべく小さく選ぶ
3. 提案分布$q(z)$から乱数yを生成
4. $[0,1]$の一様分布に従う乱数uを生成
5. $k*q(y)*u \leq p(y)$を満たすuを受容する

ここで区間$[0,kq(z)]$の一様分布を求める代わりに、区間$[0,1]$の一様分布に$kq(z)$を乗じたものを等価に扱っているのが計算手法として興味深い。

結果のグラフを見ると、サンプリング結果がガンマ分布の確率密度関数と同じような曲線を描いていることが分かる。

また、PRMLの図11.4では、棄却される領域は、正規化されていない$\tilde{p}$と$kq$の2つのグラフの差分であることが示されている。

サンプル受容確率

サンプル値$\mathbf{z}$が受容される確率求める。

まずはサンプル値$\mathbf{z}$の時に受容される確率を求める。
区間$[0,kq(z)]$の一様分布の確率$1/kq(\mathbf{z})$をサンプルuについて積分する。積分範囲については、$kq(\mathbf{z}) \geqslant \tilde{p}(\mathbf{z})$の条件より、$\tilde{p}$より大きい値域は考えなくてもいい。

\begin{eqnarray}
p(受容|\mathbf{z}) = \int_0^{\tilde{p}} \frac{1}{kq(\mathbf{z})} du = \frac{\tilde{p}}{kq(\mathbf{z})}
\end{eqnarray}

よって提案分布$q$がサンプル値$\mathbf{z}$を与える時に受容する確率は

\begin{eqnarray}
q(\mathbf{z}) \frac{\tilde{p}(\mathbf{z})}{k q(\mathbf{z})}
\end{eqnarray}

であり、これを全てのサンプル値の範囲について積分したものがサンプル受容確率になる。

\begin{eqnarray}
\int q(\mathbf{z}) \frac{\tilde{p}(\mathbf{z})}{k q(\mathbf{z})} d\mathbf{z} = \frac{1}{k} \int \tilde{p} d\mathbf{z}(\mathbf{z}) 
\end{eqnarray}

欠点

高次元の変数のサンプリングが必要な場合、受容率が極めて低くなる。

自己正規化重点サンプリング

棄却サンプリングと同様に、分布$p(\mathbf{z})$からのサンプリングが難しい場合に、サンプル$(\mathbf{z}^{1}, \cdots \mathbf{z}^{T}) \sim p(\mathbf{z})$を取得する手法。

$f(\mathbf{z})$が大きい領域からサンプリングしても、その領域の$p(\mathbf{z})$が小さければ期待値には寄与しない。
また、$p(\mathbf{z})$が大きい領域からサンプリングしても、その領域の$f(\mathbf{z})$が小さければ期待値には寄与しない。
そこで最初から$f(\mathbf{z})p(\mathbf{z})$が大きくなる領域を重点的にサンプリングすることで、計算効率を高くするのが自己正規化重点サンプリング法である。

• 仮定

正規化されていない$\tilde{p}(\mathbf{z})$が分かっている。
$$p(\mathbf{z}) = \frac{1}{Z_p}\tilde{p}(\mathbf{z})$$

また、提案分布$q(\mathbf{z})$についても正規化されていない値しか取得できないとする。
$$q(\mathbf{z}) = \frac{1}{Z_q}\tilde{q}(\mathbf{z})$$

• 手法

上記の仮定の下で期待値を計算する

\begin{eqnarray}
\int f(\mathbf{z})p(\mathbf{z})d\mathbf{z}
&=& \int f(\mathbf{z}) \frac{p(\mathbf{z})}{q(\mathbf{z})} p(\mathbf{z}) d\mathbf{z} \\
&=& \frac{Z_q}{Z_p} \mathbb{E}_{q(\mathbf{z})} \Biggl[ f(\mathbf{z}) \frac{\tilde{p}(\mathbf{z})}{\tilde{q}(\mathbf{z})} \Biggl]
\end{eqnarray}

ここでT個のサンプルT個のサンプル$(\mathbf{z}^{1}, \cdots \mathbf{z}^{T})$を選んで、単純モンテカルロ法と同じように和で近似する。ここで$\tilde{p}(\mathbf{z}) / \tilde{q}(\mathbf{z}) \equiv \mathbf{w}^{(t)}$として

\begin{eqnarray}
\int f(\mathbf{z})p(\mathbf{z})d\mathbf{z}
&\approx& \frac{Z_q}{Z_p} \frac{1}{T} \sum_{t=1}^T f(\mathbf{z}^{(t)}) \mathbf{w}^{(t)}
\end{eqnarray}

を得る。

この$\mathbf{w}^{(t)}$は重要度重みと呼ばれ、意図しない分布からのサンプリングの影響を補正する。
棄却サンプリングと異なり、サンプリングしたサンプルの全てを保持しているのが、自己正規化重点サンプリング法の特徴である。

正規化項についても近似的に求めることができる。

\begin{eqnarray}
\frac{Z_q}{Z_p}
&=& \int \frac{\tilde{p}(\mathbf{z})}{Z_q} d\mathbf{z} \\
&=& \int \frac{\tilde{p}(\mathbf{z})}{\tilde{q}(\mathbf{z})} q(\mathbf{z}) d\mathbf{z} \\
&=& \mathbb{E}_{q(\mathbf{z})} \Biggl[  \frac{\tilde{p}(\mathbf{z})}{\tilde{q}(\mathbf{z})} \Biggl] \\
&\approx& \frac{1}{T} \sum_{t=1}^T \mathbf{w}^{(t)}
\end{eqnarray}

マルコフ連鎖モンテカルロ法（MCMC）

マルコフ連鎖モンテカルロ法（MCMC）は、棄却サンプリングや自己正規化重点サンプリングと異なり、高次元空間においてもサンプリングを行える手法である。

MCMCとはサンプルを取り出したい事後分布が定常分布に収束するように遷移確立を設計する手法である。
各用語について、以下で説明する。

1次マルコフ連鎖

確率には、各試行が独立で過去の試行に影響されないものと、独立ではないものがある。
例えば、サイコロを振るのは独立試行である。今回のサイコロの目の出る確率が、前回のサイコロの目によって影響を受けることはない。

独立試行の場合、N個の観測値 $\mathbf{Z}=(z^{(1)}, z^{(2)}, \cdots, z^{(n)})$ に対して以下が成立する。

$$p(\mathbf{Z}) = p(z^{(1)})p(z^{(2)})\cdots p(z^{(n)})$$
$$p(z^{(i)}|z^{(j)}, z^{(k)}, \cdots ,z^{(l)})=p(z^{(i)})$$

それに対して、今回の状態が何になるかの確率が、前回の状態に依存するのをマルコフ連鎖と呼ぶ。この場合、先の等式は成立しない。

$$p(\mathbf{Z}) \neq p(z^{(1)})p(z^{(2)})\cdots p(z^{(n)})$$

今回の状態がj回前までの状態の影響を受ける場合、j重マルコフ過程と呼ぶ。
今回の状態が、前回の状態のみの影響をうけるならば1次マルコフ連鎖と呼び、以下が成立する。

$$p(z^{(m+1)}|z^{(1)}, \cdots ,z^{(m)}) = p(z^{(m+1)}|z^{(m)})$$

定常分布（不変分布）

1次マルコフ連鎖に従う変数の周辺確率は、1つ前の変数の周辺確率を用いて表現できる。

$$
p(\mathbf{z}^{(t+1)}) = \sum_{\mathbf{z}^{(t)}}p(\mathbf{z}^{(t-1)}| \mathbf{z}^{(t)})p(\mathbf{z}^{(t)})
$$

ここで遷移確率$\mathcal{T}(\mathbf{z}^{(t-1)}, \mathbf{z}^{(t)})$を以下のように定義しておく。

$$\mathcal{T}(\mathbf{z}^{(t-1)}, \mathbf{z}^{(t)}) \equiv p(\mathbf{z}^{(t)}| \mathbf{z}^{(t-1)})$$

1次マルコフ連鎖の全てのステップでの周辺確率が等しい場合、以下の定常分布$p_*(\mathbf{z})$の式が成り立つ。

\begin{eqnarray}
p_*(\mathbf{z}) = \int \mathcal{T}(\mathbf{z}', \mathbf{z}) p_*(\mathbf{z}') d \mathbf{z}'
\end{eqnarray}

詳細釣り合い条件

釣り合い条件は、$p_*(\mathbf{z})$が定常分布になるための十分条件である。

\begin{eqnarray}
p_*(\mathbf{z}) \mathcal{T}(\mathbf{z}, \mathbf{z}') = p_*(\mathbf{z}') \mathcal{T}(\mathbf{z}', \mathbf{z})
\end{eqnarray}

$p_*(\mathbf{z})$が釣り合い条件式を満たすならば、定常分布となることは、上式の両辺を$\mathbf{z}'$で積分すれば分かる。

\begin{eqnarray}
(lhs) &=& \int p_*(\mathbf{z}) p(\mathbf{z}' | \mathbf{z}) d \mathbf{z}' \\
&=& \int p_*(\mathbf{z}) \frac{p(\mathbf{z}' , \mathbf{z})}{p(\mathbf{z})} d \mathbf{z}' \\
&=& p_*(\mathbf{z}) \frac{1}{p(\mathbf{z})} \int p(\mathbf{z}' , \mathbf{z}) d \mathbf{z}' \\
&=& p_*(\mathbf{z})
\end{eqnarray}

\begin{eqnarray}
(rhs) &=& \int \mathcal{T}(\mathbf{z}', \mathbf{z}) p_*(\mathbf{z}') d \mathbf{z}'
\end{eqnarray}

エルゴード性

本書P85参照

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

メトロポリス・ヘイスティングス法はMCMCの基本的な手法である。

• 仮定

正規化されていない$\tilde{p}(\mathbf{z})$が分かっている。
$$p(\mathbf{z}) = \frac{1}{Z_p}\tilde{p}(\mathbf{z})$$

• 手法

1. 提案分布$q(\mathbf{z}|\mathbf{z}^{(t)})$を作成する。一般的に提案分布にはサンプリングしやすいものを選ぶ
2. 提案分布から次のサンプル候補として$\mathbf{z}_*$をサンプリングする
3. 以下の確率を計算する

\begin{eqnarray}
\mathrm{min} \Biggl( 1, \frac{\tilde{p}(\mathbf{z}_*) q(\mathbf{z}^{(t)}|\mathbf{z}_*)}{\tilde{p}(\mathbf{z}^{(t)}) q(\mathbf{z}_*|\mathbf{z}^{(t)})} \Biggl)
\end{eqnarray}

4.区間[0,1]の一様乱数$u$を求め、上記で計算した確率の値の方が$u$を上回っていれば受容し、次のサンプル$\mathbf{z}^{(t+1)} = \mathbf{z}_*$とする。それ以外の場合は、$\mathbf{z}^{(t+1)} = \mathbf{z}^{(t)}$として、再度サンプル候補のサンプリングを行い、同じ工程を繰り返す

ここで、確率式の以下の部分を比率$r$と呼ぶ。

\begin{eqnarray}
\frac{\tilde{p}(\mathbf{z}_*) q(\mathbf{z}^{(t)}|\mathbf{z}_*)}{\tilde{p}(\mathbf{z}^{(t)}) q(\mathbf{z}_*|\mathbf{z}^{(t)})}
\end{eqnarray}

分母は、次のサンプル点候補が $z_{*}$ になる確率と、次のサンプル点候補が$z_{*}$の時に現在のサンプル点が$z^{(t)}$になる確率の積なので、「現在のサンプル点における分布$p$の確率密度」を表わす。
分子は、現在のサンプル点が$z^{(t)}$になる確率と、現在のサンプル点が$z^{(t)}$の時に次のサンプル点候補が$z_{*}$になる確率の積なので、「次のサンプル点候補における分布$p$の確率密度」を表わす。

（追記）
この記載は誤りで、分母が「次のサンプル点候補における分布$p$の確率密度」を表わし、分子が「現在のサンプル点における分布$p$の確率密度」を表わしていると考えるべきでした。
それなら次のサンプル点候補の確率が高い程、それを受容する確率が増えるという自然な解釈ができます。

最初に誤った解釈をしたのは、

\tilde{p}(\mathbf{z}^{(t)}) q(\mathbf{z}_*|\mathbf{z}^{(t)})

を

p(\mathbf{z}_* )=p(\mathbf{z}^{(t)})p(\mathbf{z}_*|\mathbf{z}^{(t)})

と同様に考えたことによるものでした。（追記終わり）

実際の計算においては、この2つの確率密度の比を求めることが大事であり、分母や分子を正確に計算する必要はない。
後述する実例でも、だいたいの比率が分かれば近似計算には十分であり、上記の数式通りの厳密な計算で$r$を求めてはいない。

数値計算を行う際は、pythonにおけるmcmcライブラリを利用すると簡単だが、ここでは解説しない。

提案分布が対象である場合は、メトロポリス法と呼ばれるものになる。

遷移確率

メトロポリス・ヘイスティングス法では、遷移確率は以下のようになる。

\begin{eqnarray}
\mathcal{T}(\mathbf{z}, \mathbf{z}') 
&=& q(\mathbf{z} | \mathbf{z}') \mathrm{min} \Biggl( 1, 
\frac{\tilde{p}(\mathbf{z}') q(\mathbf{z} | \mathbf{z}')}
{\tilde{p}(\mathbf{z}) q(\mathbf{z}' | \mathbf{z})} 
\Biggl)
\end{eqnarray}

詳細釣り合い条件

MCMCとは分布が定常分布に収束するように遷移確率を設計する手法である。

実際に、定常分布が成立するならば、上で求めた遷移確率は詳細釣り合い条件を満たしているはずである。
以下にそれを証明する。

\begin{eqnarray}
p(\mathbf{z}) \mathcal{T}(\mathbf{z}, \mathbf{z}') 
&=& p(\mathbf{z}) q(\mathbf{z} | \mathbf{z}') \mathrm{min} \Biggl( 1, \frac{\tilde{p}(\mathbf{z}') q(\mathbf{z}|\mathbf{z}')}{\tilde{p}(\mathbf{z}) q(\mathbf{z}'|\mathbf{z})} \Biggl) \\
&=& p(\mathbf{z}) q(\mathbf{z} | \mathbf{z}') \mathrm{min} \Biggl( 1, \frac{p(\mathbf{z}') q(\mathbf{z}|\mathbf{z}')}{p(\mathbf{z}) q(\mathbf{z}'|\mathbf{z})} \Biggl) \\
&=& \mathrm{min} ( p(\mathbf{z}) q(\mathbf{z}'|\mathbf{z}) , p(\mathbf{z}') q(\mathbf{z}|\mathbf{z}') ) \\
&=& \mathrm{min} ( p(\mathbf{z}') q(\mathbf{z}|\mathbf{z}') , p(\mathbf{z}) q(\mathbf{z}'|\mathbf{z}) ) \\
&=& p(\mathbf{z}') q(\mathbf{z}|\mathbf{z}') \mathrm{min} \Biggl(1, \frac{p(\mathbf{z}) q(\mathbf{z}'|\mathbf{z})}{p(\mathbf{z}') q(\mathbf{z}|\mathbf{z}')} \Biggl) \\
&=& p(\mathbf{z}') q(\mathbf{z}|\mathbf{z}') \mathrm{min} \Biggl(1, \frac{\tilde{p}(\mathbf{z}) q(\mathbf{z}'|\mathbf{z})}{\tilde{p}(\mathbf{z}') q(\mathbf{z}|\mathbf{z}')} \Biggl) \\
&=& p(\mathbf{z}') \mathcal{T}(\mathbf{z}', \mathbf{z}) 
\end{eqnarray}

例：2次元ガウス分布からのサンプリング

kenmatsu4さんが【統計学】マルコフ連鎖モンテカルロ法(MCMC)によるサンプリングをアニメーションで解説してみる。において、実際に2次元ガウス分布からのサンプリングを計算されている。
数式などの説明は[基礎からのベイズ統計学 輪読会資料 第4章 メトロポリス・ヘイスティングス法(https://www.slideshare.net/matsukenbook/4-56002293)に詳しい。

記事を見ると、サンプリングを繰り返したサンプル点が2次元ガウス分布になっていることが分かる。

プログラムのアルゴリズムは以下のようになっている。

1. 提案分布には、2次元ガウス分布と一様分布のどちらかを選択。提案分布から次のサンプル点候補（次の候補位置）をサンプリングする。
2. 現在位置における目標分布の確率密度(に比例した数値)と次の候補位置における目標分布の確率密度(に比例した数値)を求め、比率$r$の計算をする。比率の時に説明した通り、両者の比がおおよそ分かればいいので『(に比例した数値)』の計算で十分である
3. 比率が $r \geq 1$の場合は必ず受容。$0 \leq r < 1$ の場合は区間$[0, 1]$の一様分布から乱数$u$を生成し、$u \leq r$ ならば受容する

例：2次元ガウス分布からのサンプリング

棄却サンプリングの例でも紹介したopabinia2さんの機械学習に詳しくなりたいブログ　マルコフ連鎖モンテカルロ法 - メトロポリス・ヘイスティングス法（２）において多変量正規分布からのサンプリングを解説されている。
先のアルゴリズムに従い、以下のことを行っている。提案分布が対称であるため、メトロポリス法で計算を行っている。

1. 提案分布として2次元ガウス分布を採用し、次のサンプル点候補をサンプリング
2. サンプル点候補の分布$p(z)$と目標分布の現在点の分布$p(x)$の比$p(z)/p(x)$を計算し、区間$[0,1]$の一様分布からの乱数と比較して、受容するかどうかを判定

ハミルトンモンテカルロ法

ハミルトンモンテカルロ法（もしくはハイブリッドモンテカルロ法）は詳細釣り合い条件を満たす遷移確率を作成する手法の1つである。

ハミルトン方程式

解析力学においてハミルトニアンは運動エネルギーと位置エネルギーの和、つまりシステムの全エネルギーを表わす。

$$
\mathrm{H}(\mathbf{z}, \mathbf{p})
= \mathrm{U}(\mathbf{z}) + \mathrm{K}(\mathbf{p})
$$

ハミルトン方程式はハミルトニアンの変分を取ることで求まる。

\begin{eqnarray}
\delta \mathrm{H} 
&=& \frac{\partial \mathrm{H}}{\partial \mathbf{z}} \delta \mathbf{z} + \frac{\partial \mathrm{H}}{\partial \mathbf{p}} \delta \mathbf{p} 
\end{eqnarray}

ラグランジアン$\mathcal{L} = \mathrm{K}(\mathbf{p}) - \mathrm{U}(\mathbf{z})$を用いるとハミルトニアンは

$$
\mathrm{H}(\mathbf{z}, \mathbf{p})
= \mathbf{p}\frac{\partial \mathbf{z}}{\partial \tau} - \mathcal{L}
$$
となるので、これも変分を取り

\begin{eqnarray}
\delta \mathrm{H} 
&=& \delta \mathbf{p} \frac{\partial \mathbf{z}}{\partial \tau} + \mathbf{p} \delta \frac{\partial \mathbf{z}}{\partial \tau}
- \frac{\partial \mathcal{L}}{\partial \mathbf{z}} \delta \mathbf{z} - \frac{\partial \mathcal{L}}{\partial \frac{\partial \mathbf{z}}{\partial \tau}} \delta \frac{\partial \mathbf{z}}{\partial \tau}
\end{eqnarray}

これに

\begin{eqnarray}
\frac{\partial \mathcal{L}}{\partial \frac{\partial \mathbf{z}}{\partial \tau}} = \mathbf{p}
\end{eqnarray}

\begin{eqnarray}
\frac{\partial \mathcal{L}}{\partial \mathbf{z}} 
= \frac{\partial \mathbf{p}}{\partial \tau}
\end{eqnarray}

を代入して

\begin{eqnarray}
\delta \mathrm{H} 
&=& \frac{\partial \mathbf{z}}{\partial \tau} \delta  \mathbf{p}
- \frac{\partial \mathbf{p}}{\partial \tau} \delta \mathbf{z}
\end{eqnarray}

となるので、最初の変分を取った式と比較してハミルトン方程式が求まる。

リープフロッグ（カエル跳躍）法

時間$t=1$をN等分して展開するという特徴がある。
詳細は本書p88参照

リウヴィルの定理

ハミルトンの方程式に従う位相空間は体積が不変になる。
ここでの位相空間とは$\mathbf{z}$軸と$\mathbf{p}$軸からなる空間のことを指す。

位相空間の流れ場を以下で与える。

$$
\mathbf{V} = \Bigl( \frac{d \mathbf{z}}{d \tau}, \frac{d \mathbf{p}}{d \tau} \Bigl)
$$

この流れ場の発散は0になる。ここでの式展開にはハミルトン方程式を代入した。

\begin{eqnarray}
\mathrm{div} \mathbf{V} 
&=& \sum_i \Bigl\{  \frac{\partial }{\partial z_i} \frac{dz_i}{d \tau} + \frac{\partial }{\partial p_i} \frac{dp_i}{d \tau} \Big\} \\
&=& \sum_i \Bigl\{  \frac{\partial }{\partial z_i} \frac{\partial \mathrm{H}}{\partial p_i} + \frac{\partial }{\partial p_i} \frac{\partial \mathrm{H}}{\partial z_i} \Big\} \\
&=& 0
\end{eqnarray}

よって位相空間の体積は変化しないことが分かる。

体積が不変であると、変数変換に伴うヤコビアンの計算が省略できるという利点がある。
つまり変数$X=(x_1, x_2, \cdots, x_n)$を変数$Y=(y_a, y_2, \cdots, y_n)$に変換した時、以下を満たす。

$$
dx_1 dx_2 \cdots dx_n = dy_1 dy_2 \cdots dy_n
$$

サンプリングアルゴリズムの適用

p89 (4.21)より、全エネルギーがハミルトニアンである位相空間上での同時分布は以下で与えられることが分かる。

$$
p(\mathbf{z}, \mathbf{p}) = \exp (- \mathrm{H}(\mathbf{z}, \mathbf{p}))
$$

以下の手順でメトロポリス法を適用できる

1. ガウス分布$\mathcal{N}(\mathbf{p}|\mathbf{0}, \mathbf{I})$から運動量$\mathbf{p}$をサンプリングする
2. 初期状態$(\mathbf{z}, \mathbf{p})$にリープフロッグ法を適用して$(\mathbf{z}_* , \mathbf{p} _* )$を求める　
3. 比率rを計算する。比率rからは確率$\mathrm{min}(1, r)$が求まる。

\begin{eqnarray}
r &=& \frac{p(\mathbf{z}_* , \mathbf{p}_*)}{p(\mathbf{z}, \mathbf{p})} \\
&=& \exp (- \mathrm{H}(\mathbf{z}_* , \mathbf{p}_*) + \mathrm{H}(\mathbf{z}, \mathbf{p}))
\end{eqnarray}

4.区間[0,1]の一様乱数$u$を求め、上記で計算した確率の値の方が$u$を上回っていれば受容し、次のサンプル$\mathbf{z}^{(t+1)} = \mathbf{z}_*$とする。それ以外の場合は、$\mathbf{z}^{(t+1)} = \mathbf{z}^{(t)}$として、再度サンプル候補のサンプリングを行い、同じ工程を繰り返す

ランジュバン動力学法

分子動力学法での実際のハミルトンモンテカルロ法を適用してみる。
ここでは『格子状の場の理論 $12.2』の解説を参照した。

分子動力学を以下のハミルトニアンで表す。$S(\mathbf{z})$は状態$\mathbf{z}$の確率を決める作用とする。

\begin{eqnarray}
\mathrm{H} &=& \frac{1}{2} \mathbf{p}^2 + S(\mathbf{z})
\end{eqnarray}

1. 初期状態を$(\mathbf{z}, \mathbf{p})$とする

2. リープフロッグ法でN分割したステップを計算する。途中式は『格子状の場の理論』を参照されたい。最終ステップのみを以下に示す。

\begin{eqnarray}
\mathbf{p}' &=& \mathbf{p}(N) = \mathbf{p}(N-1) + \epsilon \mathbf{z}(N - 1/2) \
\mathbf{z}' &=& \mathbf{z}(N) = \mathbf{z}(N-1/2) -\frac{1}{2} \epsilon \frac{\partial S(\mathbf{p}(N))}{\partial \mathbf{p}(N)}
\end{eqnarray}
```

1. 比率$r=\exp(-(\mathrm{H}(\mathbf{z}', \mathbf{p}') - \mathrm{H}(\mathbf{z}, \mathbf{p}))$を計算する
2. $\mathrm{H}(\mathbf{z}', \mathbf{p}') - \mathrm{H}(\mathbf{z}, \mathbf{p})<0$ならば受容し、状態を$(\mathbf{z}, \mathbf{p})$から$(\mathbf{z}', \mathbf{p}')$に更新する。受容しない場合は、前の状態$(\mathbf{z}, \mathbf{p})$をそのまま保持する
3. それ以外の場合は、確率$\mathrm{min}(1, r)$を計算し、区間$[0,1]$の一様分布からの乱数$u$と比較し、確率が乱数を上回っていれば受容する

例：2次元イジングスピン系の熱力学量

MCMCは統計力学で使用される事例が多い。
統計力学への適用例としてsci_Haruさんの[Pythonによる科学・技術計算] ２次元イジングスピン系の熱力学量のメトロポリス法によるモンテカルロシミュレーションを紹介する

2次元イジングモデルのハミルトニアンの統計力学的な意味には立ち入らず、これを『2値変数であるスピン$s_i$の分布からハミルトンモンテカルロ法でサンプリングを行う問題』と解釈しておく。

プログラムは以下のアルゴリズムに従っている。

1. 任意のスピン配列$A_k$を作成して、そこから変数$s_i$をサンプリングする
2. $s_i$の値を反転させたものを、新しい$A_{k+1}$の候補$s_{trial}$とする
3. $s_{trial}$のエネルギー$E_{trial}$と$s_i$のエネルギー$E_i$を求め、$E_{trial}-E_i<0$ならば$A_{k+1}=s_{trial}$とする
4. それ以外の場合は、比率$\exp(- \beta (E_{trial} - E_i))$を計算し、区間$[0,1]$の一様分布からの乱数$u$と比較し、確率が乱数を上回っていれば受容する。受容するなら$A_{k+1}=s_{trial}$として、棄却するなら何も更新しない

例：クォーク色力学での実例

クォーク色力学（以降「QCD」）では、偶数個のクォークがある場合の格子QCDでハミルトンモンテカルロ法を用いる。
『格子状の場の理論 $12.3』を参照されたい。

ギブスサンプリング

確率分布$p(\mathbf{Z})$の集合$\mathbf{Z}$が大きすぎる場合に、M個の部分集合$\mathbf{Z}= { \mathbf{Z}_1, \mathbf{Z}_2, \cdots, \mathbf{Z}_M }$に分割する手法をギブスサンプリングと呼ぶ。

\begin{eqnarray}
\mathbf{Z}_1 &\sim& p(\mathbf{Z}_1 | \mathbf{Z}_1, \mathbf{Z}_2, \mathbf{Z}_{M-1}, \cdots, \mathbf{Z}_M) \\
\mathbf{Z}_2 &\sim& p(\mathbf{Z}_2 | \mathbf{Z}_2, \mathbf{Z}_3, \mathbf{Z}_{M-1}, \cdots, \mathbf{Z}_M) \\
&.& \\
&.& \\
&.& \\
\mathbf{Z}_M &\sim& p(\mathbf{Z}_M | \mathbf{Z}_1, \mathbf{Z}_2, \cdots, \mathbf{Z}_{M-2}, \mathbf{Z}_{M-1}) \\
\end{eqnarray}

本書のP93では、ギブスサンプリングの比率が常に１になるので、サンプリングしたサンプル点が必ず受容されることを説明している。

アルゴリズムは以下になる。

1. $\mathbf{Z}$を適当な初期値に設定する
2. $\mathbf{Z}$から$\mathbf{Z}_i$を選ぶ
3. 条件付き確率密度に従い、新しい$\mathbf{Z} _* $を発生させる
4. $\mathbf{Z}_i$を新しい$\mathbf{Z} _*$で置き換える（必ず受容されるので判定を行う必要はない）
5. 2に戻り同じことを繰り返す

例：2次元ガウス分布を分割してのサンプリング

『続わかりやすいパターン認識』のA.5章は、2次元ガウス分布にギブスサンプリングを適用しているので参照されたい。

ManInMLさんが2次元正規分布でギブスサンプリングする【Python】において、上記を数値計算しているので併せて参照されたい。

参考

• Pythonによるモンテカルロ法入門

モンテカルロ法についてpythonで計算しながら詳細に解説している。
棄却サンプリングと重点サンプリングについても解説があり、他のどの教科書よりも分かりやすいかもしれない。

関連記事

『ベイズ深層学習』勉強会資料　4章の2.1~2.4
期待値伝播法
『ベイズ深層学習』勉強会資料　5章の1