Markov Chains and Monte Carlo Methods 08日目

Ioana A. Cosma and Ludger Evers, Markov Chains and Monte Carlo Methods
http://users.aims.ac.za/~ioana/notes.pdf
CC-by-nc-sa 2.5
http://creativecommons.org/licenses/by-nc-sa/2.5/za/legalcode

もしかして： chap.3 いらない？

Chapter 4. The Gibbs Sampler

4.1 introduction

importance samplingでは$f$から直接サンプリングせずに$E_f[h(X)]$を求めたが,性質の良いinstrumental ditributionをみつけるのは特に高次元に置いて困難になる． この章で議論するサンプリング法は,$f$がstationary distributionであるようなMarkov chainを設計することが最終目標である. こうした技術を総称してMarkov Chain Monte Carlo (MCMC)とよぶ. $\mathbf{x} = (x_1,..., x_p), f(\mathbf{x})$をサンプルを生成したいdistributionとして,$\{\mathbf{X}^{(i)}\}_{i=1}^n$が,$f$をstationary distributionにもつMarkov chainであるようにする. このとき$\{\mathbf{X}^{(t)}\}$は従属で,$\mathbf{X}^{(t)}$は$t$の極限で$f(\mathbf{x})$の正確なサンプルとなる.　
$f(\mathbf{x})$からサンプリングをすることが困難でも,full conditional distributions
$$f_{X_j|X_{-j}}(\cdot | x_1,..., x_{j-1}, x_{j+1}, ..., x_p) \text{ where } X_{-j} = (x_1, .., x_{j-1}, x_{j+1},...,x_p)$$
が全ての$j$について効率的にサンプリングできるとき,Gibbs sampler が使える.
記述がまちまちだが，Gibbs samplerによって生成される列はあるMarkov chainの一つのrealizationつまりsample pathである

4.2 Algorithm

Algorithm 4.1 ((Systematic sweep) Gibbs sampler)

$(X_1^{(0)}, ..., X_p^{(0)})$から初めて, $t=1, 2, ...$に
1, $X_1^{(t)} \sim f_{X_1|X_{-1}} (\cdot | X_2^{(t-1)},..., X_p^{(t-1)})$ を取る
$\vdots$
j, $X_j^{(t)} \sim f_{X_j|X_{-j}} (\cdot | X_1^{(t)},..., X_{j-1}^{(t)}, X_{j+1}^{(t-1)}, ..., X_p^{(t-1)})$ をとる
$\vdots$
p. $X_p^{(t)} \sim f_{X_p|X_{-p}} (\cdot | X_1^{(t)},...,X_{p-1}^{(t)})$　をとる
を繰り返す.

Gibbs samplerはreversible でない. Liu et al.(1995)はreversibleなchainを返すアルゴリズムを開発した.

Example 4.2 (Random sweep Gibbs sampler)

$(X_1^{(0)}, ..., X_p^{(0)})$から初めて, $t=1, 2, ...$に
1. $\{1, ..., p\}$から,(例えばuniformで) $j$を選んで,
2. $X_j^{(t)} \sim f_{X_j|X_{-j}} (\cdot | X_1^{(t-1)},...,X_{j-1}^{t-1},X_{j+1}^{t-1},...,X_p^{(t-1)})$をとって,すべての$i \neq j$に$X_i^{(t)}=X_i^{(t-1)}$とする.

4.3 The Hammersley-Clifford Theorem

Gibbs samplerの基礎であるfull conditionalはjoint distributionを一意に決定するという著しい特徴が有る(Hammersley and Cliford).

Definition 4.1 (Positivity condition)

density $f(x_1,..., x_p)$とmarginal density $f_{X_i}(x_i)$をもつ分布とがpisitivityをもつ
$\Leftrightarrow [\forall x_1,…,x_p.\ (f_{X_i}(x_i)>0 \Rightarrow f(x_1,…,x_p)>0)]$

positivityは,joint distribution $f$の台が$f_{X_i}$の台たちのデカルト積であるということである.

#### Theorem 4.2 (Hammersley-Clifford)

$(X_1,..., X_p)$がpositivityをみたし,joint densityは$f(x_1,..,x_p)$とする. このとき任意の$(\xi_1, ..., \xi_p) \in supp(f)$に,
$$f(x_1,...,x_p) \propto \prod_{j=1}^p \frac{f_{X_j|X_{-j}}(x_j|x_1,..,x_{j-1}, \xi_{j+1},...,\xi_p)}{f_{X_j|X_{-j}}(\xi_j|x_1,...,x_{j-1},\xi_{j+1},...,\xi_{p})}$$

proof.

$$f(x_1,...,x_p) = f_{X_p|X_{-p}}(x_p|x_1,...,x_{p-1})f(x_1,..,x_{p-1})$$
であって,$x_p$を$\xi_p$に置き換えても成立する.
$$f(x_1,...,x_{p-1}, \xi_p) = f_{X_p|X_{-p}}(\xi_p|x_1,...,x_{p-1})f(X_1,...,x_{p-1})$$
したがって
$$\begin{aligned} f(x_1,..., x_p) &= f(x_1,...,x_{p-1}) f_{X_p|X_{-p}}(x_p|x_1,...,x_{p-1}) \\ &=\underline{f(x_1,...,x_{p-1})}_{=f(x_1,...,x_{p-1},\xi_p)/f_{x_p|x_{-p}}(\xi_p|x_1,..,x_{p-1})} \ \ \ {f_{X_p|X_{-p}}(x_p|x_1,...,x_{p-1})} \\ \vdots \\ &= f(x_1,...,x_{p-1},\xi_p)\frac{f_{X_1|X_{-1}}(x_1|\xi_2,...,\xi_p)}{f_{X_1|X_{-1}}(\xi_1|\xi_2,..,\xi_2)} \cdots \frac{f_{X_p|X_{-p}}(x_p|x_1,...,x_{p-1})}{f_{X_p|X_{-p}}(\xi_p|x_1,...,x_{p-1})} \end{aligned}$$
よって成立. positivity conditionが分母が$0$でないことを保証する.

Hammersley-Cliffford theoremはjoint probability distributionの存在を,任意のconditionの選び方にも保証するわけではない. このような問題はBayesian modelingで，prior distributionの設定に問題が有る時によく起きる.例えば
$X_1|X_2 \sim expo(\lambda X_2), X_2|X_1 \sim expo(\lambda X_1)$とする. Hammersley-Cliffordから
$$f(x_1, x_2) \propto\frac{f_{X_1|X_2}(x_1|\xi_2)f_{X_2|X_1}(x_2|x_1)}{f_{X_1|X_2}(\xi_1|\xi_2)f_{X_2|X_1}(\xi_2|x_1)} \propto \exp(-\lambda x_1 x_2)$$
しかし$\int \int \exp(-\lambda x_1 x_2)dx_1dx_2$は無限であって,$f(x_1, x_2)$がPDFとなるような分布は存在しない.

4.4 Convergence of the Gibbs sampler

$f(x_1,...,x_p)$が実際にGibbs sampler(この節ではalg. 4.1とする)で生成されるMarkov chainのstationary distributionであることを確かめる. まず,Gibbs samplerによって生成されるtransition kernelを議論する.

Lemma 4.3

Gibbs samplerのtransition kernelは
$\begin{aligned}&K(\mathbf{x}^{(t-1)}, \mathbf{x}^{(t)}) = \\ &f_{X_1|X_{-1}}(x_1^{(t)}|X_2^{(t-1)},..,x_p^{(t-1)}) \cdot \\ & f_{X_2|X_{-2}}(x_2^{(t)}|X_1^{(t)}, x_3^{(t-1)},...,x_p^{(t-1)}) \cdots \\ & f_{X_p|X_{-p}}(x_p^{(t)}|x_1^{(t)},...,x_{p-1}^{(t)})\end{aligned}$

proof.

$$\begin{aligned} &P(\mathbf{X}^{(t)} \in \mathcal{X} | \mathbf{X}^{(t-1)}=x^{(t-1)}) = \int_\mathcal{X} f_{(\mathbf{X}^{(t)}|\mathbf{X}^{(t-1)})}(\mathbf{x}^{(t)}|\mathbf{x}^{(t-1)})d\mathbf{x}^{(t)} \\ &=f_{X_p|X_{-p}}(x^{(t)}_1|x_2^{(t-1)}, ..., x_p^{(t-1)}) f_{X_2|X_{-2}}(x_2^{(t)}|x_1^{(t)}, x_3^{(t-1)}, ..., x_p^{(t-1)}) \cdots \\ &f_{X_p|X_{-p}}(x_p^{(t)}|x_1^{(t)}, ..., x_{p-1}^{(t)}) d\mathbf{x}^{(t)} \end{aligned}$$

Proposition 4.4 証明略

$f(x_1, ..., x_p)$はたしかに生成されるMarkov chain$(\mathbf{X^{0}, X^{(1)}, ...})$のstationary distributionである.

以上, Gibbs samplerが生成するMarkov chainは$f$をstationary distributionにもつことが言えた. Theorem 1.19では,Markov chain がstationary distributionに収束する十分条件がirreducibleかつaperiodicであることを見たが,Gibbs samplerが生成するMarkov chainがこれを満たすかは議論の余地が有るし,実際満たさない.

Example 4.3 (Reducible Gibbs sampler)

$C_1 := \{(x_1, x_2)| \|(x_1, x_2)-(1, 1)\|\leq 1\}, C_2:=\{(x_1, x_2)|\|(x_1, x_2)-(-1, -1)\| \leq 1\}$
とし,$f$は$C_1 \cup C_2$上一様分布のPDFとする.このとき,$X_1^{(0)} < 0$なる初期値から開始したGibbs samplerはfig. 4.2のように,$C_2$の点のみを取り出してしまう.

これは生成されたMarkov chainがirreducibleでないために起きる. 次の命題はGibbs samplerの生成するMarkov chain のirreducibilityの十分条件を与える. より弱い条件の十分な命題もある(Robert and Casella, 2004, Lemma 10.11)

Proposition 4.5

$f(x_1, ..., x_p)$がpositivity conditionを満たすとき,Gibbs samplerはirreducibleかつrecurrentなMarkov chainを生成する.

proof.

$\mathcal{X} \subset supp(f)$は$\int_\mathcal{X} f(x_1^{(t)}, ..., x_p^{(t)})d(x_1^{(t)},..,x_p^{(t)})>0$を満たすとする.
$$\int_\mathcal{X} K(\mathbf{x}^{(t-1)}, \mathbf{x}^{(t)})d\mathbf{x}^{(t)}=\int_\mathcal{X} \underline{f_{X_1|X_{-1}}(x_1^{(t)}|x_2^{(t-1)},...,x_p^{(t-1)})}_{>0} \cdots \underline{f_{X_p|X_{-p}}(x_p^{(t)},x_1^{(t)}, ..., x_{p-1}^{(t)})}_{>0} d\mathbf{x}^{(t)}>0$$
が,positivityよりconditional densityが正の値であることから言える.よって$\{\mathbf{X}^{(t)}\}_t$はstrongly f-irreducibleで, prop. 1.28から,Markov chainはまたrecurrentである.

さらに,エルゴード性の帰結としてTh. 4.6が得られる.

Theorem 4.6

Gibbs sampler によって生成されるMarkov chainがirreducibleかつrecurrentであるとき,可積な$h:E \rightarrow \mathbb{R}$について
$$\frac{1}{n}\lim_{n\rightarrow \infty} \sum_{t=1}^n h(\mathbf{X}^{(t)}) \rightarrow E_f [h(\mathbf{X})]$$
がほとんどすべての初期値$\mathbf{X}^{(0)}$で成立する.

これが$E_f[h(\mathbf{X})]$を,生成した一つのMarkov chainの平均によって推測することを正当化する.

Example 4.6

$(X_1, X_2) \sim N_2 \left( \left(\begin{array}{} \mu_1 \\ \mu_2 \end{array} \right) , \left(\begin{array}{} \sigma_1^2 & \sigma_{12} \\ \sigma_{12} & \sigma_{2}^2 \end{array} \right)\right)$
について,$P(X_1 \geq 0, X_2 \geq 0)$をGibbs samplerによって計算する.
marginal distributionは$X_1 \sim N(\mu_1, \sigma_1^2), X_2 \sim N(\mu_2, \sigma_2^2)$である
conditional distibution $X_1|X_2 = x_2$と$X_2|X_1=x_1$は正規分布の多項式表現から
$$X_1|X_2=x_2 \sim N(\mu_1+\sigma_{12}/\sigma_2^2(x_2-\mu_2), \sigma_1^2-(\sigma_{12})^2/\sigma_2^2)$$
$$X_2|X_1=x_1 \sim N(\mu_2 +\sigma_{21}/\sigma_1^2(x_1-\mu_1), \sigma_2^2 -(\sigma_{21})^2 /\sigma_1^2)$$
よってGibbs samplerが,$t=1,2...$に
1. $X_1^{(t)} \sim N(\mu_1+\sigma_{12}/\sigma_2^2(X_2^{(t-1)}-\mu_2), \sigma_1^2-(\sigma_{12})^2/\sigma_2^2)$を取る.
2. $X_2^{(t)} \sim N(\mu_2 +\sigma_{21}/\sigma_1^2(X_1^{(t)}-\mu_1), \sigma_2^2 -(\sigma_{21})^2 /\sigma_1^2)$を取る.

を繰り返してMarkov chain $\{X^{(t)}=(X^{(t)}_1, X^{(t)}_2)\}_t$を生成する.
$\mu_1 = \mu_2 = 0, \sigma_1^2 = \sigma_2^2 = 1, \sigma_{12}=0.3$とするとき,fig.4.4はひとつのsample pathの例である. さらにTh. 4.6により,$P(X_1\geq 1, X_2 \geq 0)$が$(X_1^{(t)}\geq 0, X_2^{(t)}\geq 0)$の$1$から$t$までの平均によって推測できる. $t$を横軸として平均をプロットしたのがfig. 4.3である.

Markov性から$(X^{(0)}, ...)$は従属であり，普通は正の相関を持つ. $\{X^{(t)}\}$の相関が大きいほどMarkov chainはゆっくりと動く(slowly mixingという). Gibbs samplerにおいても, $X_j$が正であれ負であれ強く相関しているときにはそのような現象が見られる. ex.4.5はその例である.

Example 4.5 (Sampling from a highly correlated bivariate Gaussian)

4.4 の例で,ただ$\sigma_{12}=0.99$にした場合に,$\rho(X_1, X_2)=0.99$である. このときGibbs samplerはslower mixingで,fig. 4.5,からわかるとおり,収束が非常に遅い.