サンプリング法3〜MCMC〜

補足.
・Monte Carlo法とは、乱数を用いた数値計算全般を指します。つまり、MCMCではMarcov chainで乱数を取得→何かしらの計算(Monte Carlo)をします。


・マルコフ連鎖では離散的な状態のみを扱いますが、MCMCでは連続的な状態も扱います。

説明. (次元の大きいベクトルについて受容棄却法ではうまくサンプリングができないことについて)

目標分布$p(\vec z)$：$\N_d(\vec 0, \sigma_p^2 I)$
提案分布$q(\vec z)$：$\N_d(\vec 0, \sigma_q^2 I)$ （$\sigma_p^2 \gt \sigma_q^2$）
とします。




$q(z),k$は、
$$\begin{aligned} k q(\vec z) \geqq p(\vec z) ~~~~~ \small{(\forall \vec z \in \mathrm{\vec R}^d)} \end{aligned}$$を満たす必要がありますが、これを満たす最小の$k$は、
$$\begin{aligned} k = \frac{p(0)}{q(0)} = (\frac{\sigma_q}{\sigma_p})^d \end{aligned}$$となりますが、この時のサンプリングの受容率は、
$$\begin{aligned} \int \frac{p(\vec z)}{k q(\vec z)} q(\vec z) d \vec z &= \frac{1}{k} \int p(\vec z) d \vec z \\ &= \frac{1}{k} \\[10px] &= (\frac{\sigma_p}{\sigma_q})^d \\[10px] &\overset{ d \to \infty}\longrightarrow 0 \end{aligned}$$となります。
（参照：C.M. ビショップ『パターン認識と機械学習』丸善出版, 2012.）




つまり、次元の大きいベクトルをサンプリングしようとした時、その受容率が0に近くなってしまい、うまく（正確には『効率的な』）サンプリングができないことになります。

すたどく

まとめると、受容棄却法の様に独立に（1つずつ別個に）乱数を得るアプローチは高次元の状況で困難を生じてしまうわけです。そこで系列的（$\neq$独立）にサンプリングする方法が求められますが、MCMCはその様なサンプリング法の1つです。

for i=0,1,…,do
$~~~~$# 提案分布$q(\vec z| \vec z^{(i)})$からサンプルする
$~~~~$sample $\vec z^* \sim q(\vec z| \vec z^{(i)})$

$~~~~$# 下記の通り$A(\vec z^*, \vec z^{(i)})$を定める
$~~~~$# ただし$p()$は目標分布を表す
$~~~~$$A(\vec z^*, \vec z^{(i)}) = \min\{ 1, \frac{p(\vec z^*) q(\vec z^{(i)} | \vec z^*)}{p(\vec z^{(i)} ) q(\vec z^* | \vec z^{(i)})} \}$

$~~~~$# 下記の通り$\vec z^{(i+1)}$を更新する
$~~~~$・$\vec z^{(i+1)}$←$\vec z^*$ with probability $A(\vec z^*, \vec z^{(i)})$
$~~~~$・$\vec z^{(i+1)}$←$\vec z^{(i)}$ otherwise (probability $1-A(\vec z^*,\vec z^{(i)})$)

$~~~~$end

最終的にiが一定以降(1000以降など)のサンプルを採用する。

すたどく

この様なアルゴリズムの下では$\vec z$の定常分布は目標分布$p(\vec z)$となります。
（この証明に興味がある方は＜補足. MCMC＞を参照ください）

補足.
原点対称な分布（正規分布、$\t$分布、一様分布）$K()$に対して、提案分布$q()$を、
$$\begin{aligned} q(\vec z^{*} | \vec z^{(i)}) = K(\vec z^*-\vec z^{(i)}) \end{aligned}$$と定めることが多いです。




この時、
$$\begin{aligned} q(\vec z^{(i)} | \vec z^*) = q(\vec z^* | \vec z^{(i)}) \end{aligned}$$が成立するため、前出のアルゴリズムにおいて、
$$\begin{aligned} A(\vec z^*, \vec z^{(i)}) = \min\{ 1, \frac{p(\vec z^*)}{p(\vec z^{(i)} )} \} \end{aligned}$$となります。




よって、提案分布$q()$からランダムサンプルした$z^*$について、$p(\vec z^*) \geqq p(\vec z^{(i)})$、となれば、$A(\vec z^*,\vec z^{(i)}) = 1$、より$\vec z^{(i)}$→$\vec z^*$と必ず推移します。一方で、$p(\vec z^*) \lt p(\vec z^{(i)})$、となれば、$A(\vec z^*, \vec z^{(i)}) = \frac{p(\vec z^*)}{p(\vec z^{(i)} )}$、より$\vec z^{(i)}$→$\vec z^*$と必ず推移するわけではありません。（確率$\frac{p(\vec z^*)}{p(\vec z^{(i)} )} (\leqq 1)$で推移する）




言い換えると、目標分布$p()$において確率が高い点をより取りたがり、目標分布$p()$において確率が低い点をより取りたがらないことになります。その結果、ゆくゆくは乱数が目標分布$p()$からサンプリングされることになります。

for i=0,1,…,do
$~~~~$# Gibbs samplingではベクトルを構成する要素を1つずつ更新する
$~~~~$# $\vec z^{(i)} = (z_1^{(i)},\cdots,z_K^{(i)})$、と$K$要素からなるものとする
$~~~~$# $p()$は目標分布を表す
$~~~~$sample $z_1^{(i+1)} \sim p(z_1 | z_2^{(i)},\cdots,z_K^{(i)})$
$~~~~$sample $z_2^{(i+1)} \sim p(z_2 | z_1^{(\color{red}{i+1})},z_3^{(i)},\cdots,z_K^{(i)})$
$~~~~$・・・（以下同様）・・・
$~~~~$sample $z_K^{(i+1)} \sim p(z_K | z_1^{(\color{red}{i+1})},\cdots,z_{K-1}^{(\color{red}{i+1})})$

$~~~~$# 上記プロセスにおいては、『目標分布$p()$において、ベクトル$\vec z$を構成する要素の1つは、残りの要素が条件づけられていればsample可能（＊）』という条件下で話を進めています

$~~~~$end

最終的にiが一定以降(1000以降など)のサンプルを採用する。

すたどく

この様なアルゴリズムの下では$\vec z$の定常分布は目標分布$p(\vec z)$となります。
（この証明に興味がある方は＜補足. MCMC＞を参照ください）

目標:
特定の目標分布$p()$からsampleする。


MCMCによる解決策：
特定の目標分布$p()$を定常分布に持つようなマルコフ連鎖を構成して、十分長い時間シミュレーションする。マルコフ連鎖は定常分布に近づくので（ほぼ）特定の目標分布$p()$からのサンプリングができる。


適切なマルコフ連鎖の構成の方法：
マルコフ連鎖は次の状態への推移確率$R(\vec z^{(i+1)}\mid \vec z^{(i)})$で特徴づけられる。この$R$をマルコフ連鎖の定常分布が特定の目標分布$p()$に一致するように定めたい。


この条件を数式で書くと、
$$\begin{aligned} p(\vec z^{(i+1)}) = \int R(\vec z^{(i+1)} \mid \vec z^{(i)}) p(\vec z^{(i)})d \vec z^{(i)} \end{aligned}$$
この条件を満たす$R$の構成方法として、Metropolis-HastingsやGibbs Samplingがありました。

すたどく

最後の式は$p()$が定常分布であることの定義式です。なかなかこの式の意味を直感的に理解することは難しいので説明してみます。まず、$\vec z^{(i+1)}$を$\vec z^{(i)}$と推移確率$R(\vec z^{(i+1)}\mid \vec z^{(i)})$によって得るとします。一般に『（$a,b$の同時分布）＝（$a$の分布）$\times$（$a$を条件づけた時の$b$の条件付き分布）』であることから、$(\vec z^{(i)}, \vec z^{(i+1)})$の同時分布は$p(\vec z^{(i)}) R(\vec z^{(i+1)} \mid \vec z^{(i)})$となり、上式右辺の被積分項となります。上式右辺ではこれを$\vec z^{(i)}$について積分してるので、右辺は$\vec z^{(i+1)}$の従う分布を表します。それが左辺、つまり、$p(\vec z^{(i+1)})$に一致します。以上より、$\vec z^{(i)}, \vec z^{(i+1)}$いずれの従う分布も$p()$となるわけです。

まとめ.

• MCMCとはマルコフ性を利用した系列的なサンプリング法である。
  
  
• MCMCの代表的アルゴリズムとしてMetropolis-HastingsとGibbs samplingがあり、サンプリング方式や適用条件は異なるが、「定常分布が目標分布に一致する様にアルゴリズムが組まれている」という点で共通している。