Markov Chains and Monte Carlo Methods 09日目

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

Gibbs samplerが終わったし後は流れで

Chapger 5. The Metropolis-Hastings algorithms

Gibbs samplerを使うにはfull conditionalsから効率よくサンプリングできなければならなかった． また,$X_1,..,X_p$に強い相関が有るとき,収束が遅くなるという欠陥が有る. これを克服するのがMetropolis-Hastings法である. rejection samplingのように,新しい$\mathbf{X}^{(t+1)}$を$\mathbf{X}^{(t)}$によって決まる局地的な分布に従って受理または棄却し,得られた$\{\mathbf{X}^{(t)}\}$をあるMarkov chainのpathと考える．

Algorithm 5.1 (Metropolis-Hastings)

$\mathbf{X}^{(0)} = (X_1^{(0)}, ..., X_p^{(0)})$を初期値として,$t=1,2,...$に
1. $\mathbf{X} \sim q(\cdot | \mathbf{X}^{(t-1)})$をとる.
2. $$\alpha(\mathbf{X}|\mathbf{X}^{(t-1)})=\min\left\{1, \frac{f(\mathbf{X})\cdot q(\mathbf{X^{(t-1)}|\mathbf{X}})}{f(\mathbf{X}^{(t-1)})\cdot q(\mathbf{X}|\mathbf{X}^{(t-1)})} \right\}$$
を計算する.
3. 確率$\alpha(\mathbf{X}|\mathbf{X^{(t-1)}})$で$\mathbf{X}^{(t)}=\mathbf{X}$とし,そうでなければ$\mathbf{X}^{(t)} = \mathbf{X}^{(t-1)}$とする.

Lemma 5.2

Metroplis-Hastingsのtransition kernelは
$$K(\mathbf{x}^{(t-1)}, \mathbf{x}^{(t)}) = \alpha(\mathbf{x}^{(t)}|\mathbf{x}^{(t-1)})q(\mathbf{x^{(t)}|\mathbf{x}^{(t-1)}})+(1-a(\mathbf{x^{(t-1)}}))\delta_[\mathbf{x}^{(t-1)}(\mathbf{x}^{(t)})$$
ここで$\delta_{\mathbf{x}^{(t-1)}}(\cdot)$はDirac-mass とする.

Proposition 5.3

Metropolis-Hastingsはdetailed balance
$$K(\mathbf{x}^{(t-1)}, \mathbf{x}^{(t)})f(\mathbf{x}^{(t-1)}) = K(\mathbf{x}^{(t)}, \mathbf{x}^{(t-1)})f(\mathbf{x}^{(t)})$$
をみたす.したがって$f(\mathbf{x})$は生成されるMarkov chainのstationary distributionであり,しかもMarkov chainはreversibleである.

さらに,chainがirreducibleかつaperiodicならばMarkov chainは任意のinitial distributionでstationary distributionに収束する.

Theorem 5.5 (Ergodic theorem)

Metropolis-Hastingsによって生成されるMarkov chainがirreducibleであるとき,可測な$h$に
$$\lim_{n \rightarrow \infty} \sum_{t=1}^n h(\mathbf{X}^{(t)})\rightarrow E_f[h(\mathbf{X})]$$
が任意の初期値に成立する.

5.3 The random walk Metropolis algorithm

Metropolis-Hastingsの特別な場合に,random walk Metropolisがある. Metropolis-Hastingsにおける提案分布$q$を，$\mathbf{X} = \mathbf{X}^{(t-1)} + \epsilon, \epsilon \sim g$に変えた場合である. ただし$g$は対称性をもつ分布とする(i.e. $g(x)=g(-x)$). このとき,$\mathbf{X}-\mathbf{X}^{(t-1)} \sim g \sim \mathbf{X}^{(t-1)}-\mathbf{X}$が対称性より言えるから,
$$\alpha(\mathbf{X}|\mathbf{X}^{(t-1)})=\min\left\{1, \frac{f(\mathbf{X})\cdot q(\mathbf{X^{(t-1)}|\mathbf{X}})}{f(\mathbf{X}^{(t-1)})\cdot q(\mathbf{X}|\mathbf{X}^{(t-1)})} \right\}=\min \left\{1, \frac{f(\mathbf{X})}{f(\mathbf{X}^{(t-1)})} \right\}$$
つまり
#### Algorithm 5.2 (Random walk Metropolis)

$\mathbf{X}^{(0)} = (X_1^{(0)}, ..., X_p^{(0)})$を初期値として,以下を$t=1,2,...$に繰り返す
1. $\epsilon \sim g$をとって,$\mathbf{X} = \mathbf{X}^{(t-1)} + \epsilon$
2. $$\alpha(\mathbf{X}|\mathbf{X}^{(t-1)}) =\min \left\{1, \frac{f(\mathbf{X})}{f(\mathbf{X}^{(t-1)})} \right\}$$
3. 確率$\alpha(\mathbf{X|X^{(t-1)}})$で$\mathbf{X}^{(t)}=\mathbf{X}$とし,そうでなければ$\mathbf{X}^{(t)}=\mathbf{X}^{(t-1)}$とする.

Example 5.2 (Bayesian probit model)

帝王切開による出産の際の感染の有無の調査(table 1)

$n_i$人の患者のうちの感染数$Y_i$を推測する.
$Y_i \sim Bin(n_i, \pi_i), \ \ \pi = \Phi(\mathbf{z}_i' \mathbf{\beta})$を仮定する.
ただし$\mathbf{z}_i=(1, z_{i1}, z_{i2}, z_{i3})$, $\Phi$は$N(0, 1)$のCDFとする.
$\mathbf{\beta}$のprior distributionには$N(0, I/\lambda)$を使う. $\beta$のposterior densityは
$$f(\beta|y_1, ...,y_n) \propto \left( \prod_{i=1}^n \Phi(z_i^T \beta )^{y_i} (1-\Phi(z_i^T \beta))^{n_i-y_i} \right) \exp(-\frac{\lambda}{2} \sum_{j=0}^3 \beta_j^2)$$
random walk Metropolis を使ってこのposterior からサンプリングを行う. 初期値$\beta^{(0)}$を適当に決めて,$t=1, 2, ...$で
1. $\epsilon \sim N(0, \Sigma), \beta = \beta^{(t-1)} + \epsilon$を計算する
2. $\alpha(\beta|\beta^{(t-1)})=\min \{ 1, f(\beta|Y_1,...,Y_n)/f(\beta^{(t-1)}|Y_1,...,Y_n)\}$を計算する
3. 確率$\alpha$で$\beta^{(t)} = \beta$, そうでなければ$\beta^{(t)} = \beta^{(t-1)}$とする.

を繰り返す. $\Sigma=0.08 I$とする.
table 5.2 とfig. 5.3が50000サンプルを取ったときの一つの結果である. ただし最初の10000サンプルは排除している.

5.4 Choosing the proposal distribution

Metropolis-Hastingsの効率はproposal distribution $q$の選び方に強く依存している. $q(\cdot|\mathbf{X}^{(t-1)})$は$\mathbf{X}^{(t-1)}$から離れたところに山をもっているのが望ましいし,また$\alpha(\mathbf{X}|\mathbf{X}^{(t-1)})$が大きいことが望ましい. この二つは相反する要求である.
経験的に,1か2次元のモデルでは受理率は1/2程度がよく,3次元以上のモデルでは1/4程度が良いと知られている.