Markov Chains and Monte Carlo Methods 10日目

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が終わったし後は流れで

• Chapter 7. State-space models and the Kalman filter algorithm
  • 7.1 Motivation
  • 7.2 State-space models
    • 7.2.1 Inference problems in SSMs

Chapter 7. State-space models and the Kalman filter algorithm

7.1 Motivation

現実世界では,観測は時間的に離散的な列で行われて,その観測が行われるたびに以前の観測たちから興味有る対象の量($x$)を推測することになる. これをon-line inferenceという. 観測データが$y_1, y_2, ..., y_t, ...$と時刻$t$に沿って与えられて,各時刻で興味有る対象$x_t$を推測することを考える. $x_t$の事前分布$p(x_t)$が$t$によって変化するようなモデルをDynamic Modelという.
Dynamic modelの例には,レーダーによる観測からの飛行機の監視(場所,速度を推測する)やノイズの有る音声データからの発話の認識(発音された単語を推測する) などがある. これらの問題を扱うのに適した方法の一つがSequential Monte Carlo (SMC)である. SMCはiterativeではないMCMCの一種で,$t-1$における$x$の分布をapproximate sampleで表現し,それを$t$における分布の表現で再利用する.

7.2 State-space models

SSMでは,根底となり,観測できないMarkov process(state process) $\{X_t\}$と,観測される過程(observation procss)$\{Y_t\}$からなる.
observation: $y_t = a(x_t, u_t) \sim g(\cdot| x_t, \phi)$
hidden state: $x_t = b(x_{t-1}, v_t) \sim f(\cdot|x_{t-1}, \theta)$
$u,v$はノイズを表す変数で,$\phi, \theta$は既知である. $p(x_1)$が初期状態$x_1$の分布とする. state processは$p(x_t|x_1,...,x_{t-1})=p(x_t|x_{t-1})=f(x_t|x_{t-1},\theta)$だからMarkov chainである. さらに$p(y_t|x_{1:t}, y_{1:t-1})=p(y_t|x_t)=g(y_t|x_t, \phi)$すなわち$y_t$の分布は$x_t$それ以前の観測値と独立である(fig. 7.1).

$x_{1:t}$で$x_1,...,x_t$を表し,$y_{1:t}$も同様である. 簡単のため$\theta$や$\phi$との依存関係の記述を省略して$f(\cdot|x_{t-1}), g(\cdot|x_t)$などと書く.

7.2.1 Inference problems in SSMs

$$\begin{aligned} p(x_{1:t}, y_{1:t}) &= p(x_1)g(y_1|x_1) \prod_{i=2}^t p(x_i,y_i|x_{1:t-1},y_{1:t-1}) \\ &=p(x_1)g(y_1|x_1)\prod_{i=2}^t f(x_i|x_{i-1})g(y_i|x_i) \end{aligned}$$
またBayesの定理から,
$$p(x_{1:t}|y_{1:t}) \propto p(x_{1:t}|y_{1:t-1})g(y_t|x_t) = p(x_{1:t-1}|y_{1:t-1})f(x_t|x_{t-1})g(y_t|x_t)$$
が成立する.
- $p(x_t|y_{1:t})$からのサンプルを取ることをfilgeringという
- $p(x_{1:t}|y_{1:t})$からのサンプルを取ることをsmoothingという
これらについて議論する.

filtering, smoothingはともに$p(x_{1:t}|y_{1:t})$の扱いやすさが問題となる. ほとんどのSSMでは,この分布はnormalizing constantしかわからない.その例外が
- transitionとobservationが離散的である場合(Hidden Markov modelという)には,再帰的なアルゴリズムが使える
- 関数$a, b$が線形であり,ノイズ$u_t, v_t$が正規分布である時(linear Gaussian SSMという),これを解くアルゴリズムをKalman filterという.

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程度が良いと知られている.

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,からわかるとおり,収束が非常に遅い.

Markov Chains and Monte Carlo Methods 07日目

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

3.3 Importance sampling

rejection samplingでは,target $f(x)$のかわりにinstrument $g(x)$からサンプリングし,$f(x)$に合致しなそうなサンプルを棄却することで$f(x)$からのサンプリングを行った. importance samplingでは$g(x)$からのサンプルを重み付けして$f(x)$からのサンプリングを実現する. impotrance samplingの最も重要な基礎は
$$P(X\in A) = \int_A f(x)dx = \int_A g(x) \underline{\frac{f(x)}{g(x)}}_{:w(x)}dx = \int_A g(x)w(x)dx$$
が$f(x)>0 \Rightarrow g(x)>0, a.e.$なる全ての$g$に成立することである. これはまた，任意の可測関数$h$に,
$$E_f[h(X)]=\int_S f(x)h(x)dx = \int_S g(x)\frac{f(x)}{g(x)}h(x)dx=\int_S g(x)w(x)h(x)=E_g[w(X)h(X)]$$
と一般化出来る.
$X_1,...,X_n \sim g$があって,$E_g|w(X)\cdot h(X)|$が存在するとき
$$\frac{1}{n} \sum_{i=1}^n w(X_i)h(X_i) \rightarrow^{\text{a.s.}} E_g[w(X)\cdot h(X)]$$
が大数の強法則から言える. $E_g[w(X)h(X)]=E_f[h(X)]$だから
$$\frac{1}{n}\sum_1^n w(X_i)h(X_i)\rightarrow^{\text{a.s.}} E_f[h(X)]$$
つまり$\mu = E_f[h(X)]$は
$$\tilde{\mu} =\frac{1}{n}\sum_{1}^n w(X_i)h(X_i)$$
で近似できる.
$E_g(w(X))=\int_S \frac{f(x)}{g(x)}g(x)dx = \int_S f=1$だが,$w(X_1),...,w(X_n)$の総和は必ずしも$n$ではないので,self-normalized版
$$\hat{\mu} = \frac{1}{\sum_{i=1}^n w(X_i)} \sum_{i=1}^n w(X_i)h(X_i)$$
を正当化でき，以下のアルゴリズムが導かれる.

Algorithm 3.2 (Impotrance sampling)

$\text{supp}(f\cdot h) \subset \text{supp}(g)$なる$g$を選んで,
1. $i = 1,...,n$に
(i) $X_i \sim g$を生成する
(ii) $w(X_i)= f(X_i)/g(X_i)$とする
2.
$$\hat{\mu} = \frac{\sum_{i=1}^n w(X_i)h(X_i)}{\sum_{i=1}^n w(X_i)}$$
あるいは
$$\tilde{\mu} = \frac{\sum_{i=1}^n w(X_i)h(X_i)}{n}$$
を返す.

Theorem 3.3 (Bias and Variance of Importance Sampling)

(a) $E_g(\tilde{\mu})=\mu$
(b) $var_g[\tilde{\mu}] = \frac{var_g[w(X)h(X)]}{n}$
(c) $$E_g(\hat{\mu})=\mu+\frac{\mu var_g[w(X)] - cov_g[w(X), w(X)h(X)]}{n}+O(n^{-2})$$
(d) $$var_g[\hat{\mu}] = \frac{var_g[w(X)h(X)]-2\mu cov_g [w(X), w(X)h(X))+\mu^2 var_g[w(X)]]}{n}+O(n^{-2})$$

proof.

(a) $$E_g[\frac{1}{n} \sum_{i=1}^n w(X_i)h(X_i))] = \frac{1}{n}\sum_i E_g[w(X_i)h(X_i)]=E_f[h(X)]$$
(b) $$var_g[\frac{1}{n}\sum_{i=1}^n w(X_i)h(X_i)] = \frac{1}{n^2}\sum_i var_g(w(X_i)h(X_i))=\frac{var_g [w(X)h(X)]}{n}$$
(c, d) 略

この定理から,$\tilde{\mu}$は不偏だが分散が大きく,$\hat{\mu}$は不偏でないが分散が$\tilde{\mu}$より小さいことがわかる. さらに,$f(x) = C\pi(x)$とすると
$$\hat{\mu} = \frac{\sum w(X_i)h(X_i)}{\sum w(X_i)} = \frac{\sum \frac{f(X_i)}{g(X_i)}h(X_i)}{\sum \frac{f(X_i)}{g(X_i)}} = \frac{\sum \frac{C\pi(X_i)}{g(X_i)}h(X_i)}{\sum \frac{C\pi(X_i)}{g(X_i)}} = \frac{\sum \frac{\pi(X_i)}{g(X_i)}h(X_i)}{\sum \frac{\pi(X_i)}{g(X_i)}}$$
だから,$C$がわからなくとも$\hat{\mu}$は計算できる.
$g$はsupportの条件を満たせば良いが,普通$\tilde{\mu}$の分散を有限にするように選ぶ.これは以下の２つの条件のどちらかが成立すればよい.
- $f(x) < Mg(x) \text{ and } var_f[h(X)]<\infty$ ・・・・ $g$はrejection samplingにも使える
- $S$がコンパクトで,$f$が$S$上有界

さらに$g$が最良である,すなわち$var[\tilde{\mu}]$が最小になるような$g$の選び方を考える.

Theorem 3.4 (Optimal proposal) 証明略

$var[\tilde{\mu}]$を最小にする$g^*$は
$$g^* (x) = \frac{|h(x)|f(x)}{\int_S |h(t)|f(t)dt}$$
で与えられる.

Corollary

importance samplingはsuper-efficientである. すなわちTh. 3.4 による$g^*$を使うと,$\tilde{\mu}$は$f$から直接サンプリングしたときの分散よりも小さくなる.
$\because$
$$\begin{aligned} n\cdot var_f [\frac{h(X_1)+\cdots +h(X_n)}{n}] &= E_f(h(X)^2)-\mu^2 \geq_{\text{Jensen's inequality}} (E_f[|h(X)|])^2 -\mu^2 \\ &=(\int_S |h(x)| f(x)dx )^2-\mu^2 =n \cdot var_{g^*}[\tilde{\mu}] \end{aligned}$$

$g^*$のnormalisaton constantを知らなければならず,また$g^*$からのサンプリングが難しいことも有るので,$g^*$に近い別の$g$をinstrumental として使うことが有る．

Example 3.5 (Computing $E_f|X| \text{ for } X \sim t_3$)

$X$は自由度3のt分布($t_3$とする)に従うとして,$E_f[X]$をMonte Carlo methodで計算する. 以下の３つの方法が考えられる.
- X_1,…,X_nを$t_3$から直接サンプリングし, $\frac{1}{n} \sum_{i=1}^n |X_i|$で推測する
- $t_1$(Cauchy分布に同値)をinstrumentalにしてimportance samplingする.
- $N(0, 1)$をinstrumental にしてimportance samplingする.このとき$var[\tilde{\mu}]=\infty$

2つのinstrumentalとtargetのグラフはfig. 3.4の通り.

Markov Chains and Monte Carlo Methods 06日目

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

assignment 1

Q and A

2.

答案.

(a) 数学的帰納法を使う.$m=1$で確かに成立.
m=mで成立を仮定し,$m=m+1$での成立を示すのは$K^{(m+1)}=K^{m+1}=K K^{m}=K K^{(m)}$を使えばたやすい.
(b) stationary distributionを$\mu$とすると,$\mu^T K = \mu^T$が成り立つから,
$\mu^T(K-I) =0$を解く.
その解は$span\{(\beta, \alpha)^T\}$. $\mu$の要素は全て非負で和が1になることを考えれば,求めるdistibutionは$\mu = \frac{1}{\alpha + \beta}(\beta, \alpha)^T$.
(c) detailed-balanceを満たしていることを見れば良い.
$\mu_1 K_{1, 2} = \frac{1}{\alpha + \beta} \beta \alpha = \frac{1}{\alpha + \beta} \alpha \beta = \mu_2 K_{2,1}$だから,たしかに成立.

3.

(a, b) $\{\{1, 2, 3, 6, 7, 8 \}, \{ 4, 9\}, \{5\}, \{10\}\}$
それぞれ(recurrent, period=2), (transient, aperiodic), (transient, aperiodic), (transient, aperiodic)

4.

(a) 任意の$i < j$について,$P(X_{t+m_1}=j|X_i=i)>0, P(X_{t+m_2}=i|X_{t}=j)$を満たすような$m_1,m_2$が有ることを示せば良い. 定義より$P(X_{t+(j-i)}=j|X_t=i)=2^{-(j-i)}>0, P(X_{t+i}=i|X_{t}=j)=\underline{(1/2)}_{(1)}\cdot \underline{2^{-i}}_{(2)}>0$から示せた.
((1): $j$から$0$に遷移する確率, (2): $0$から$i$に遷移する確率)
(b) (i)recurrence

$S$が一つのcommunicating classだから,$0$がrecurrentであることを示せば十分である.
$k_{00}^{(t)} = 1/2 + 1/2 \cdot 1/2 + 1/2 \cdot (1/2 \cdot 1/2) + ... + 1/2 \cdot \underline{(1/2 \cdot ... \cdot 1/2)}_{t\text{個}} = 1-2^{-(1+t)}$であって,
$\lim_t k^{(t)}_{00}=1$ だから$\sum_t k_{00}^{(t)} \rightarrow \infty$すなわちrecurrent. 以上より示せた.

(ii) aperiodicity

$K_{ii}^{(m)} \begin{cases} = 0 \ \ \ (m < i) \\ \geq 2^{-(i+1)} \ \ (m \geq i) \end{cases}$から,$\{m \geq 1| K_{ii}^{(m)}>0\}= \{i, i+1, ...\}$.GCDは1であり,aperiodic. これが全ての$i$に成立する.

(c) transition matrixは
$$K = \left(\begin{array}{} 1/2 & 1/2 & 0 & 0 &\cdots \\ 1/2 & 0 & 1/2 & 0 &\cdots \\ 1/2 & 0 & 0 & 1/2 & \cdots \\ \vdots \end{array} \right)$$
である. stationaryな$\mu = (\mu_0 ,\mu_1, ...)^T$を計算する.
$\mu^T K = \mu^T$を解けばよい.
$\mu_0 = (\mu_0 + \mu_1 + \cdots) / 2$
$\mu_1 = \mu_0 /2$
$\mu_2 = \mu_1 / 2$
$\vdots$
であって,$\mu_i = \mu_{i-1}/2 \ \ ( i \geq 1)$ゆえに$\mu_i = \mu_0 \cdot 2^{-i}$.
$\sum_{i \geq 0} \mu_i = 2\mu_0 = 1$だから,$\mu_0=1/2$.以上より$\mu_i = 2^{-(i+1)}$.
stationary $\mu=(1/2, 1/4, 1/8, ...)^T$が示せた.

6.

(a) $$p(s|X_t = x_t) = \underline{(1- k_{x_tx_t})}_{(1)}\ \underline{(k_{x_tx_t})^{s-1}}_{(2)}$$
(1): $x_t$から他のstateに遷移する確率
(2): $s-1$回,$x_t$から$x_t$に遷移する確率
(b) $P(S_t > s_0 + s| S_t > s_0)=P(S_t>s)$を示す.
$lhs = P(S_t > s_0 + s, S_t > s_0)/P(S_t > s_0) = P(S_t > s_0 + s)/P(S_t > s_0)$
$$=\frac{1-\sum_{r=1}^{s_0+s}(1-k)k^{r-1}}{1-\sum_{r=1}^{s_0}(1-k)k^{r-1}}=k^{s_0+s}/k^{s_0}=k^s=1-P(s_t\leq s_0) = P(S_t > s)=rhs$$
よって示せた.

8.

def Wright_Fisher(a, b):
    path = [a]
    tn = a + b #2N
    xt = a
    for i in range(1000):
        xt = np.random.binomial(tn, xt/tn)
        path.append(xt)

    return path

MIT OCW, Machine Learning 11日目 宿題3

Rohit Singh, Tommi Jaakkola, and Ali Mohammad. 6.867 Machine Learning. Fall 2006. Massachusetts Institute of Technology: MIT OpenCourseWare, https://ocw.mit.edu. License: Creative Commons BY-NC-SA.

assigmnemt 3
Q and A

1.

答案.

(b) $E[error_{LOOCV}(S_n)]=E[\frac{1}{n}\sum_t error_t(S_n)] = \frac{n}{n} \cdot E[(y-\hat{f}_{S_{n-1}}(\mathbf{x}))^2]=E[(y-\hat{f}_{S_{n-1}}(\mathbf{x}))^2]$
(c) $$\begin{aligned}var[error_{LOOCV}(S_n)] &= var[\frac{1}{n} (error_1(S_n) + \cdots+error_n(S_n))] \\ &=\frac{1}{n^2}var(\sum_1^n error_i(S_n)) <_{(1)} var(error_1(S_n)) \\ &= E[(y-\hat{f}_{S_{n-1}}(\mathbf{x}))^2] \end{aligned}$$
(1)の不等号を示せなければならなのだが,模範解答でも定性的に言及されただけだからもう明らかでいいと思う
(d)
$f^k_{keep}$がtraining errorを0にするというのは，任意の$t\in\{1,..,n\}$で$y_t = x^{(t)}_r$が成り立つということ. $y_t$を生成する$Y$と$x_r$を生成する$X_r$は独立で,$P(Y=1)=P(X_r=1)=0.5, P(Y=-1)=P(X_r=-1)=0.5$だから,$f^{k}_{keep}$がtraining errorを0にする確率は$P(\forall t. Y=X_r)=_{(1)}\prod_t P(Y=X_r)=1/2^n$. (1): 各サンプルの生成の独立性
$f^k_{flip}$も同様で,足し合わせると求める確率$2^{-(n-1)}$が得られる.
(e)
$\mathcal{M}_r$において$f^r_{keep}$のtraining error $\epsilon << 1/2$とする.
$f^r_{keep}$の$error_i(S_n)$を考える($i = 1,..., n$)
$\epsilon << 1/2$だから，training setから$i$番目を引いても$\mathcal{M}_r$から選ばれるestimatorは変わらず.$f^r_{keep}$のまま. したがって$\hat{f^r}_{-i} = f^r_{keep}$が任意の$i$で成立.
よって$error_{LOOCV}(S_n) = (1/n) \cdot \sum_i (y_i - f^r_{-i}(x_i))^2 = (1/n)\cdot\sum_i (y_i-f^r_{keep}(x_i))^2$
これは$f^r_{flip}$のtraining errorが$\epsilon << 1/2$の時も同じ. よって示せた.

模範解答.

(a)
r.v. $A, B$が同じdistributionをもつとき，$E[f(A)]=E[f(B)]$であるのを利用する.これは
$$E[f(A)] = \int f(x)p_A(x)dx = \int f(x)p_B(x)dx = E[f(B)]$$
からわかる.
$A,B$がそれぞれr.v.の集合であっても成立する. $A=\{S^{-1}_n, (\mathbf{x}_1, y_1)\},\ B=\{S_{n-1}, (\mathbf{x}, y)\}$とする.ただし$S_n^{-1}$でtrainし，$(\mathbf{x}_1, y_1)$を識別子, $S_{n-1}$は$n-1$個のtraining dataで$(\mathbf{x}, y)$を識別する. $A,B$は同じdistributionを持つから,与えられた四季が成立する.
(f)
training error を$\delta$とすると,classifierは$n\delta/4$回間違える. ある次元$i$においてtraining errorが$\epsilon$以下である時,すなわち間違いがせいぜい$floor(n\epsilon/4)$であるとする.間違いの回数を$k$とおくと,間違いの起こる場合の数は$nCk$通りで,まさにそこで間違いが起こる確率は$2^{1-n}$. よってtraining errorが$\epsilon$未満である確率は
$$p = \sum_{k=0}^{floor(n\epsilon/4)} nCk \frac{1}{2^{n-1}}$$
errorが$\epsilon$以上である確率は$1-p$で,$d$次元全てがそうである確率は$(1-p)^d$.これが$1/2$以下であれば少なくとも1つの次元でerrorが$\epsilon$未満となる. よって
$(1-p)^d \leq 1/2$を解くと,
$$d \geq \frac{1}{\log_2 \frac{1}{1-p}}$$

2.

答案.

(a)
$$P(S_n|\{l\}) = \sum_{\theta \in \{-1, 1\}} 2^{-1} \prod_{t=1}^n \left[ \frac{1+y_t\theta_l x_{t_l}}{2} \right]=\frac{1}{2^{n+1}} \left[\prod_{t=1}^n (1+y_tx_{t_l})+\prod_{t=1}^n(1-y_tx_{t_l}) \right]$$
(d) marginal likelihoodが減少し始めるとき

模範解答.

(b) 与えられた式は間違いに対して確率0を割り当ててしまう.
$$P(S_n|\mathcal{J})=\sum_{\theta \in \{1, -1\}^{|\mathcal{J}|}}2^{-|\mathcal{J}|}\prod \left[\frac{1}{|\mathcal{J}|}\sum_{j \in \mathcal{J}} f(y_t\theta_j , x_{tj}) \right]$$
とする.

Markov Chains and Monte Carlo Methods Computer Practical 1. Markov chains

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
Q and A

task1

def m_trans(K, m):
    # A, Bがnp.array型のときA * BやA**nはアダマール積として扱われるので,
    # np.matrix型に変換する.
    K = np.matrix(K)
    return K**m

m_trans(np.array([[0.9, 0.1], [0.3, 0.7]]), 3)

->

matrix([[ 0.804,  0.196],
        [ 0.588,  0.412]])

task2

def dist_t(K, t, lambda_0):
    return np.transpose(lambda_0).dot(m_trans(K, t))

# phone lineの例で,stationary distributionを試すと
dist_t(np.array([[0.9, 0.1], [0.3, 0.7]]), 100, np.array([0.75, 0.25]))

->

matrix([[ 0.75,  0.25]])

task3

def stationary(K):
    K = np.matrix(K)
    dim = len(K)
    b = np.zeros(dim+1)
    b[dim] = 1

    up = np.transpose(K - np.matrix(np.identity(dim)))
    mu = np.linalg.lstsq(np.array(np.concatenate((up, [np.ones(dim)]))), b)
    print(mu[0])


stationary(K)

->

[ 0.75  0.25]

task4

t = np.linspace(0, 25)
a = np.random.uniform()
init = np.array([a, 1-a])
def pow_generator(K): # Kのべき乗を生成するジェネレータ
    K = np.matrix(K)
    tmp = np.identity(len(K))
    while True:
        tmp = K.dot(tmp)
        yield tmp

gen = pow_generator(K)
process = [np.array(init.dot(next(gen)))[0] for t in range(25)]

plt.plot(process, marker=".")
plt.ylim([0, 1])
plt.show()

->

task5

def sample_path(K, init, tau):
    states = list(range(len(init))) # ここではstateは0から番号を付ける
    x0 = np.random.choice(states, p=init)
    path = [x0]
    tmp = x0

    for i in range(tau):
        suc = np.random.choice(states, p=K[tmp])
        tmp = suc
        path.append(suc)

    return path

plt.plot(sample_path(K, [0.5, 0.5], 100))
plt.show()

->

task6

path = sample_path(K, [0.5,0.5],20000)
avg = []
for i in range(1, 100):
    avg.append(np.average(path[: 200 * i]))
plt.plot(avg)
plt.show()

->

Markov Chains and Monte Carlo Methods 05日目

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

Chapter 3. Fundamental Concepts: Transformation, Rejectino, and Reweighting

3.1 Transformation methods

$U[0, 1]$の現れを生成する(サンプリングする)方法はすでに見た. CDF $F$をもつ分布からサンプリングする方法を考える. transformation methodはそのようなアルゴリズムのひとつのクラスであって,transformation methodの最も単純なアルゴリズムがInversion Methodで,generalized inverse(一般化逆関数) $F^-(u)=\inf \{x|F(x)\geq u\}$を用いる.

Theorem 3.1 (Inversion Method)

$U \sim U[0, 1]$として,$F$はあるCDFとする. $F^-(U)$のCDFはまた$F$である.

proof.

$F^-(u) \leq x\Leftrightarrow u \leq F(x)$だから,$U\sim U[0, 1]$に
$$P(F^-(U)\leq x)=P(U\leq F(x)) = F(x)$$

Example 3.1 (Exponential distribution)

パラメータ$\lambda$のexponential distribution($exp(\lambda)$)のCDFは$F_\lambda(x)=1-\exp(\lambda x)\ \ \ \ (x\geq 0)$であって,$F_\lambda^-(u)=\log(1-u)/\lambda$. inversion methodから,これで$U[0, 1]$からの現れを写像すれば$exp(\lambda)$からのサンプリングを行える.

Inversion Methodはinverse CDFが効率的に計算できる分布に対してのみ効率が良いアルゴリズムである. 例えば正規分布はCDFもその逆関数も解析的に書けない. しかし,generalised inverseでない変換によって欲しい分布を実現する方法も有る.

Example 3.2 (Box-Muller Method for Sampling from Gaussian)

$X_1, X_2 \sim N(0, 1)$,IIDとする. この2つの実数の組を平面上の点と考えるとその極座標$(R, \theta)$について,$R, \theta$は独立で,$\theta \sim U[0, 2\pi], \ R^2 \sim exp(1/2)$である.
$$X_1 = \sqrt{R^2} \cos (\theta), X_2 = \sqrt{R^2} \sin (\theta)$$
が成立するから,$U_1, U_2 \sim U[0,1]$を使って
$$X_1 = \sqrt{-2\log(U_1)} \cos(2\pi U_2), X_2 = \sqrt{-2\log (U_1)}\sin (2\pi U_2)$$
で$X_1, X_2の現れが得られる.

transformation methodは，目的とする分布以外の,扱いやすい分布からサンプリングを行い,そのサンプルたちを目的とする分布のサンプルとなるように変換する技術である. 多くの場合，そのような変換をclosed formで得ることはできず，そのような場合，目的とする分布に似ているが実は異なる分布からサンプリングを行い，不合理なサンプルを棄却することで目的とする分布のサンプリングを行う方法が有る. これをrejection samplingといい，次節であつかう.

3.2 Rejection Sampling

rejection samplingは,instrumental distributionからサンプリングし,目的の分布の点ではなさそうなサンプルを棄却する. 目的分布のPDF $f$は既知とする. rejection samplingの根底には,
$$f(x) = \int^{f(x)}_0 1 du = \int^1_0 1_{0 <u<f(x)}du$$
がある. $f(x)$を,$\{(x, u) | 0 \leq u \leq f(x)\}$における一様分布の，$x$による周辺分布と考えるのである. fig. 3.2はその概略図である.

Example 3.3 (Sampling from a Beta distribution)

$Beta(a, b)$は
$$f(x) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}x^{a-1}(1-x)^{b-1}\ \ \ \ (0 < x < 1)$$
ただし,$\Gamma(a)=\int^\infty_0 t^{a-1}\exp(-t)dt$はGamma関数である. $Beta(a, b)$のPDFは$(a-1)/(a+b-2)$をmodeとする単峰なグラフを持つ(fig. 3.2).
fig.3.2の影の部分からサンプルを取るには, ex.2.1と2.2でみたのと同じ技術を使う. つまり,明るいグレーの四角形に一様にサンプルの候補を置き,影になっている部分のみをサンプルとして保存するのである.
形式的には,$X \sim U[0,1], U \sim U[0, 2.4]$から独立にサンプリングし,$U < f(X)$となるような$(X, U)$の組のみをサンプルとする.
$$P(U<f(X)|X=x)=P(U<f(x))=f(x)/2.4$$
は,$(X, U)$の組が,$X=x$という条件のもとでサンプルになる条件付き確率である.

ex.3.3の例では,BetaのPDFが短径に覆われることを利用したが,PDFが性の値を取るrange(support, 台という)が非有界な分布にはそのまま適用できない. しかしそのような$f(x)$を,より簡単な$g(x)$によって$M\cdot g(x)$として抑えることでrejection samplingを実現できる.$g(x)$をproposal distribution(提案分布)という.

Algorithm 3.1 (Rejection sampling)

任意の$x$に$f(x)<M g(x)$が成り立つような$M\in \mathbb{R}$と$g$を与えられたとき,$f$からのsampleを以下のようにして得る.
1. $X \sim g$を得る.
2. $X$を,確率
$$\frac{f(X)}{M g(X)}$$
で受理して,受理しないときには１にもどる

proof.

$\mathcal{X}$を,棄却を考えずに$g$から得た$X$の集合とする.
$$P(X \in \mathcal{X} \text{ and is accepted}) = \int_\mathcal{X} \underline{g(x)}_{x \text{is from }g} \underline{\frac{f(x)}{Mg(x)}}_{P(X \text{ is accepted}|X=x)} dx = \frac{\int_\mathcal{X} f(x)dx}{M}$$
さらに, $S$を$X$が取りうる値全ての集合とすると$(\int_\mathcal{X} f(x)dx) \leq \int_S f(x)dx = 1$で,
$$P(X\text{ is accepted}) = P(X\in S \text{ and is accepted}) = 1/M$$ を代入すれば
$$P(X \in \mathcal{X}|X\text{ is accepted})=\frac{P(X\in \mathcal{X} \text{ and is accepted})}{P(X \text{ is accepted})} = \frac{\int_\mathcal{X}f(x)dx/M}{1/M}=\int_\mathcal{X}f(x)dx$$
よってこのアルゴリズムで生成された値たちの密度は($\mathcal{X}$が一様なら)$f$.

Remark 3.2

$f(x)=C\cdot \pi(x)$について,$\pi(x)$しかわかっていないときには
$$\frac{\pi(X)}{M\cdot g(X)}$$
によってrejection samplingを行える.

Example 3.4 (Rejection sampling from the $N(0,1)$ using a Cauchy proposal)

$N(0, 1)$とCauchy distributionのPDFはそれぞれ
$$\begin{aligned}f(x) = \frac{1}{\sqrt{2\pi}}\exp (-\frac{x^2}{2})\\ g(x) = \frac{1}{\pi(1+x^2)} \end{aligned}$$
であって,$M=\sqrt{2\pi}\exp(-1/2)$とすれば,$f(x) \leq Mg(x)$が言える. (fig.3.3)
一方で,$N(0, 1)$をproposal distributionとしてCauchy distributionをrejection samplingすることはできない. $g(x) < Mf(x)$なる$M$が存在しないためである.

Markov Chains and Monte Carlo Methods 04日目

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

• • 2.1 What are Monte Carlo Methods?
  • 2.2 Introductory examples
    • • Example 2.1 (A raindrop experiment for computing $\pi$)
      • Example 2.2 (Monte Carlo Integration)
      • Example 2.3 (Buffon’s needle)
    • 2.4 Pseudo-random numbers
      • Algorithm 2.1 (Conguruential pseudo-RNG)

2.1 What are Monte Carlo Methods?

• Stochastic integration
  積分をシミュレーションで近似する
• Monte Carlo tests
  p値をシミュレーションで近似する
• Markov Chain Monte Carlo(MCMC)
  興味有る分布に収束するMarkov chainを構成する

2.2 Introductory examples

Example 2.1 (A raindrop experiment for computing $\pi$)

$\pi$をMonte Carloによって推測する. ある雨粒が落ちる位置の確率変数を,$x$軸$X$,$y$軸$Y$とする. ある正方形$R=[-1, 1]\times[-1,1]$に一様に雨粒が落ちると仮定して,その中の単位円にも一様に雨粒が落ちる. $X, Y$がiidでuniformally distribution ,$U[-1, 1]$に従うとする.
$$P(\text{drop within circle}) = \frac{\text{area of the unit circe}}{\text{area of the square}}=\frac{\int \int_{x^2+y^2\leq 1}1 dxdy}{\int\int_{-1\leq x,y \leq 1}1 dxdy}=\pi/4$$
これは$\pi = 4 \cdot P(\text{drop within circle})$と同じ.
$n$個のraindropに対して,単位円に落ちる個数のr.v.を$Z$とすると$Z$はbinomialである.つまり
$$Z \sim B(n,p), \ \ \ p = P(\text{drop within circle})$$
$p$を最尤法での推定値は$\hat{p} = Z / n$. よって$\hat{\pi}=4\hat{p}=4\cdot \frac{Z}{n}$.
law of large numbersによって,$\hat{\pi}$がほとんど必ず$\pi$に収束する. 中心極限定理によって,例えば$n=100$として$Z \sim B(100, p)$とすれば,$Z \sim N(100p, 100p(1-p))$で近似できる. よって$\hat{p} =\hat{Z}/100 \sim(p, p(1-p)/100)$であって,$p$の95%信頼区間は
$$\left[ 0.77-1.96\sqrt{\frac{0.77(1-0.77)}{100}}, 0.77 + 1.96\sqrt{\frac{0.77(1-0.77)}{100}} \right]=[0.6875, 0.8525]$$
さらに$\pi$の95%信頼区間は$[2.750, 3.410]$.
以上やってきたことは
- $\pi$をある期待値として表現した
- 代数的な表現を,それのsample approximationに書き換えた. そのsample approximationが収束することを大数の法則で保証し,CLTによって収束の測度を議論した.

Example 2.2 (Monte Carlo Integration)

$$\int^1_0 f(x) dx \text{ with } f(x)=\frac{1}{27} (-65536x^8+262144x^7-409600x^6+311296x^5-114688x^4+16384x^3)$$
をMonte Carlo integrationすることを考える.$f([0, 1]) = [0,1]$だから,$[0,1]$上の$f$のグラフは$[0, 1] \times [0, 1]$に収まる. またraindrop experimentを考える. $f(x)=\int^{f}_0 1 dt$だから
$$\int^1_0 f(x) dx = \int^1_0 \int^{f(x)}_0 1 dt dx = \int \int _{\{(x, t): t \leq f(x)\}} 1dt dx = \frac{\int \int _{\{(x, t): t \leq f(x)}1dtdx}{\int \int _{0\leq x, t \leq 1}1 dtdx}$$
分子は$f(x)\leq y$のグラフの面積で,分母は$[0,1]\times[0,1]$の面積である. $n$個の雨粒を落として$f(x) \leq y$に落ちる確率が$\hat{p}_n$なとき,$(1-2\alpha)$信頼区間は
$$\left[\hat{p_n} - z_{1-\alpha}\sqrt{\frac{\hat{p_n}(1-\hat{p_n})}{n}},\hat{p_n} + z_{1-\alpha}\sqrt{\frac{\hat{p_n}(1-\hat{p_n})}{n}} \right]$$
だから,収束の早さは$O_P(n^{-1/2})$. 一方Riemann積分の速度は$O(n^{-1})$.
Monte Carloの場合の収束の早さは次元に依存しない一方で,他の決定論的な積分評価の場合は次元の増加とともに収束が遅くなっていくので,高次元な関数の積分でMonte Carloは威力を発揮する.

Example 2.3 (Buffon’s needle)

3本の間隔$\delta$の平行な直線で平面が区切られていて,長さ$l < \delta$の針を落とすとき,その針が直線と交わる確率はどれほどだろうか?

解答 (Buffon, 1777)

針が直線との角度$\theta$で着地したとき,針が直線と交わる$\Leftrightarrow$針の一端と直線の距離が$l\sin \theta$以下(fig. 2.5(a)). したがって
$$P(\text{intersect}|\theta) = \frac{l\sin \theta}{\delta}$$
さらに$\theta$は$[0, \pi)$上一様分布していると仮定すると
$$P(\text{intersect})=\int^\pi_0 P(\text{intersect}|\theta)\cdot \frac{1}{\pi}d\theta = \int^\pi_0 \frac{l\sin \theta}{\delta}\frac{1}{\pi}d\theta=\frac{l}{\pi\delta}\int^\pi_0 \sin\theta d\theta=\frac{2l}{\pi \delta}$$

Lazzarini,1901は$l=2.5cm, d=3cm$の場合に,1808本の針を使って$\pi \sim 3.14159292035$を算出した. これは非常に良い近似である. 力学的にMonte Carlo法を行うのは非常に時間がかかるが,電子計算機の到来によってこの欠点は克服された. しかし,例からわかるように,それぞれの実験での確率変数の現れがたしかにもとの分布から生成されていなければならないので,乱数の生成が重要になってくる.

2.4 Pseudo-random numbers

ここでは$U[0, 1]$の現れを生成するpseudo-random number generator(RNG)を考える. これには以下の性質が必要である.
- RNGの生成する値は独立である
- RNGの生成する値は$[0, 1]$にまんべんなく分布する

以下にlinear congruential generator(線形合同法)の概要を述べる. linear congruential generatorは上で述べた性質をあまり満たしていないので実践すべきではない.

Algorithm 2.1 (Conguruential pseudo-RNG)

1. $M \in \mathbb{N}, c \in \mathbb{N}_0, Z_0 \in \{1,...,M-1\}$を選ぶ
2. $i = 1,2,...$に
   $Z_i = (aZ_{i-1}+c) \mod M, X_i = Z_i / M$とする.

これは明らかに決定論的なアルゴリズムで,それぞれのパラメータを一致させれば完全に一致する出力をおこなう. また, 生成される値$\{X_i\}$は,$(X_{nk+1},...,X_{n(k+1)-1})$を$n$次元空間の点と考えることで,$n$次元立方体のテント見ることが出来る. これらの点は有限の-しばしばごく小さい数の-超平面に乗っていて,したがってまんべんなく分布していると見ることができない(fig. 2.6, fig. 2.7).

よりよいpseudo-RNGには,例えばMarsaglia and Zaman(1991)やMatsumoto and Nishimura(1998)がある.

Markov Chains and Monte Carlo Methods 03日目

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

• • • 1.4 Ergodic theorems
      • Theorem 1.30 (Ergodic Theorem) 証明略
      • Example 1.17

1.4 Ergodic theorems

Makrov chainを観測して,そのstationary distributionを推測する方法とその条件を考える. IIDな確率変数の列では大数の法則が観測した値の平均によって期待値を推測することが正当化される一方,Markov chainではErgodic theorems (エルゴード性に関する諸定理)によって似たような主張が出来る. これらはMarkov Chain Monte Carlo(MCMC) を正当化する根拠でも有る.

Theorem 1.30 (Ergodic Theorem) 証明略

$X$は$\mu$-irreducibleかつrecurrentな$\mathbb{R}^d$上の Markov chainで,stationary distribution $\mu$をもつとする. このとき$g: \mathbb{R}^d \rightarrow \mathbb{R}$があって,確率1で
$$\lim_{t\rightarrow \infty} \frac{1}{t}\sum_{i=1}^t g(X_i)\rightarrow E_\mu[g(X)]=\int_S g(x)f_\mu(x)dx$$
がほとんどすべての初期値$X_0=x$に成立する. また,$X$がHarris-recurrentであればこれは任意の初期値$X_0=x$で成立する.

左辺は$g$の時間的な平均であり,右辺は$g$の空間的な期待値と考えることが出来る.
ex. 1.17はrecurrenceやirreducible性がth. 1.30の仮定に必要であることの例である.

Example 1.17

$S= \{1, 2\}$でtransition matrix $\mathbf{K} = (\begin{array}{} 1 & 0 \\ 0 & 1\end{array})$をもつdiscrete Markov chainを考える.
任意の$\mu$がstationary distributionであり,さらにこのchainはirreducibleでもrecurrentでもない. $\mu = (\alpha, 1-\alpha)^T$ならば任意の$t$に
$$P(X_t=1)=\alpha, P(X_t=2) = 1-\alpha$$
が成立する. 一方sample pathはただ(1, 1, 1,…)か(2, 2, 2,…)しかありえず,どちらを一つだけ得たとしても$\alpha$や$X_t$についての推測はできない. このchainは$S$をexploreできないためであり,様々な推測をするには複数のsample pathが必要となる.