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

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が必要となる.

Markov Chains and Monte Carlo Methods 02日目

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

• • • 1.3 General state space Markov chains
      • Definition 1.23 (Markov chain)
      • Example 1.15 (Gaussian random walk on R\mathbb{R})
      • Definition 1.24 (Irreducibility)
      • Example 1.16 (Gaussian random walk (continued))
      • Definition 1.25 (Recurrence)
      • Definition 1.26 (Harris Recurrence)
      • Definition 1.27 (Stationary distribution)
      • Proposition 1.28 (証明略)
      • Definition 1.29 (Detailed balance)

1.3 General state space Markov chains

$|S| > |\mathbb{N}|$の場合の議論を始める. より一般的な場合にも定義できるのだが,ここでは$S=\mathbb{R}^d$とする.

Definition 1.23 (Markov chain)

$\{X_t\}$はdiscrete time stochastic processで, state spaec $S$は$S=\mathbb{R}^d$とする. $X$ がMarkov property を満たす
$\Leftrightarrow$
$$\forall A \in \mathcal{F}.\ P(X_{t+1}\in A| X_0=x_0,...,X_t=x_t)=P(X_{t+1}\in A|X_t=x_t)$$

ただし$\mathcal{F}$は$S$の可測集合族とする.
以後,Markov chainはhomogeneousとする.すなわち $(X_{t+1}\in A|X_t=x_t)$は$t$によらない. このときtransition kernel $K: S \times S \rightarrow \mathbb{R}_0^+$によって
$$P(X_{t+1}\in A| X_t=x_t) = \int_A K(x_t, x_{t+1})dx_{t+1}\ \ \ \ \ (1.3)$$
として得られる. $K(x, y)$というのは$X_t=x_t$が与えられたときの$X_{t+1}$のconditional probability densityである. def. 1.8における$K(i, j)=k_{ij}$というのはdiscrete spaceにおけるconting measureであって,(1.3)の式に合致する.
さらに
$$\begin{aligned}&P(X_{t+m}\in A|X_t=x_t) \\&= \int_A \int_S \cdots \int_S K(x_t, x_{t+1})K(x_{t+1},x_{t+2})\cdots K(x_{t+m-1},x_{t+m})dx_{t+1}\cdots dx_{t+m-1}dx_{t+m}\end{aligned}$$
だから,$m$-step transiton kernelは
$$K^{(m)}(x_0, x_m) = \int_S\cdots \int_SK(x_t,x_{t+1})\cdots K(x_{m-1},x_m)dx_{m-1}\cdots dx_1$$
であり,
$$P(X_{t+m}\in A| X_t=x_t) = \int_A k^{(m)} (x_t, x_{t+m})dx_{t+m}$$
と簡潔に書ける.

Example 1.15 (Gaussian random walk on $\mathbb{R}$)

$\mathbb{R}$上のrandom walkを
$$X_{t+1} = X_t + E_t, \ \ E_t \sim N(0, 1)$$
とすると,
$$X_{t+1}|X_{t} = x_t \sim N(x_t, 1)$$
と同値. $Et_t$は$X_0, E_1,...,E_{t-1}$と独立とし, $X_0 \sim N(0,1)$とすれば,
$$\begin{aligned} P(X_{t+1}\in A| X_t=x_t,...,X_0=x_0) &= P(E_t\in A-x_t|X_t=x_t,...,X_0=x_0) \\ &=P(E_t\in A-x_t)=P(X_{t+1}\in A|X_t=x_t) \end{aligned}$$
よって$X$はMarkov chainであり,しかも
$$P(X_{t+1}\in A| X_t=x_t) = P(E_t \in A - x_t) = \int_A \phi(x_{t+1}-x_t)dx_{t+1}$$
である.ただし$\phi(z)=\frac{1}{\sqrt{2\pi}}\exp(-\frac{1}{2}z^2)$である.したがって$K(x, x_{t+1}) = \phi(x_{t+1}-x_t)$.
$m$-step transition kernelは(1.3)によって計算できるが,非常に複雑になる. それよりも
$$X_{t+m} = X_t + E_t + ... + E_{t+m-1}$$
を利用すれば,$X_{t+m} |X_t = x_t \sim N(x_t, m)$が成立して,
$$P(X_{t+m}\in A| X_t=x_t) = P(X_{t+m}-X_t\in A-x_t) = \int_A \frac{1}{\sqrt{m}}\phi \left(\frac{x_{t+m}-x_t}{\sqrt{m}}\right)dx_{t+m}$$
したがって
$$K^{(m)}(x_t, x_{t+m}) = \frac{1}{\sqrt{m}}\phi \left(\frac{x_{t+m}-x_t}{\sqrt{m}} \right)$$
が$m$-step kernelとして得られる.

Definition 1.24 (Irreducibility)

$S$上の分布$\mu$が与えられて, Markov chain $X$が$\mu$-irreducibleである
$\Leftrightarrow$
$$P(X_{t+m}\in A| X_t=x) = \int_A K^{(m)}(x, y)dy > 0$$
が任意の$\mu (A) > 0$なる$A \in \mathcal{F}$と任意の$x \in S$に成立するような$m \in \mathbb{N}_0$が存在するとき$\mu$-irreducibleといい,特に任意の$A$に$m=1$で成立するときstrongly $\mu$-irreducibleという.

Example 1.16 (Gaussian random walk (continued))

ex. 1.15で$X_{t+1}|X_t = x_t \sim N(x_t, 1)$を見た. $P(X_{t+1}\in A|X_t=x_t)>0$が任意のnullでない$A$に成立するから,これは任意のcontinuous distributionにstrongly irreducibleである.

periodicity, recurrence, そしてtransienceのような概念を$S$が連続的なMakov chainに導入するため, atomsやsmall setsのような概念を導入する必要があって,これらはこのノートの範囲を超えるので,recurrenceのみを一般化する.
section 1.2.3で定義した$|S|\leq |\mathbb{N}|$の場合のrecurrenceとは,全てのstateがそれを初期のstateとしたとき,平均して無限回訪れられることであった. $S$が連続である場合には,あるstate一点ではなく,stateの集合たちを考えることになる. $V_A=\sum_{t=0}^\infty 1_{\{X_t \in A\}}|X_0=x$として,$A$が訪れられる回数をあらわす. expected valueを考えると
$$E[V_A|X_0=x] = E[\sum_{t=0}^\infty 1_{\{X_t \in A\}}|X_0=x] = \sum_{t=0}^\infty E(1_{\{X_t \in A\}}|X_0=x) =\sum_{t\geq 0} \int_A K^{(t)}(x, y)dy$$
である. recurrenceをMarkov chain全体に定義する前に,集合のrecurrenceを定義する.

Definition 1.25 (Recurrence)

(a) $A \subset S$がMarkov chain $X$においてrecurrentである
$\Leftrightarrow$ 任意の$x \in A$に
$$E(V_A|X_0=x) = \infty$$
(b) Markov chain $X$がrecurrent
$\Leftrightarrow$
(i) $X$はある分布$\mu$に対して$\mu$-irreducibleであって,かつ
(ii) 任意の$A \in \mathcal{F}, \mu(A) > 0$はrecurrent

定義より,recurrent setとは平均して無限回訪れられる集合であって,より強い命題に,その集合が無限回訪れられる確率が1であるというのがある. この強い性質によって定義できるrecurrenceをHarris Recurrenceという.

Definition 1.26 (Harris Recurrence)

(a) $A \subset S$が$X$にHarris-recurrentである
$\Leftrightarrow \forall x \in A P(V_A=\infty|X_0=x)=1$
(b) Markov chain $X$がHarris-recurrent
$\Leftrightarrow$
(i) $X$はある$\mu$に$\mu$-irreducible
(ii) 任意の$A \in \mathcal{F}, \mu(A) > 0$はHarris-recurrent

Harris recurrenceはrecurrenceを導くことは明らかで,discreteの場合2つは一致する.
どちらの概念も成立を証明することは非常に困難だが,あるMarkov chainがirreducibleで唯一のstationary distributionをもつならばrecurrentであるという命題を主張する. その前に,stationaryを定義する.

Definition 1.27 (Stationary distribution)

分布PDF$f_{\mu}$をもつ分布$\mu$がtransition kernel $K$をもつ$X$のstationary distibutionである
$\Leftrightarrow \forall y \in S.\ \ f_\mu(y) = \int_S f_\mu(x)K(x, y)dx$

Proposition 1.28 (証明略)

$X$が$\mu$-irreducible Markov chainで,$\mu$を唯一のstationary distibutionにもつなら,$X$はrecurrentである.

また,def. 1.27によってstationarityを確かめるのは困難なので,discreteの場合と同様により簡単な十分条件を定義する.

Definition 1.29 (Detailed balance)

transition kernel $K$がdistribution $\mu$にdetailed balanceである
$\Leftrightarrow \forall x,y \in S. \ \ f_\mu(x)K(x, y)=f_\mu(y)K(y,x)$

theorem 1.22と同様に,$\mu$によってdetailed balanceなMarkov chain $X$はtime-reversibleで$\mu$は$X$のstational distibutionである.

Markov Chains and Monte Carlo Methods 01日目

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 1. Markov Chains
  • 1.1 Stochastic processes
    • • Definition 1.1 (Stochastic process)
      • Definition 1.2 (Sample path)
      • Theorem 1.3 (Kolmogorov)
  • 1.2 Discrete Markov chains
    • 1.2.1 Introduction
      • Definition 1.4 (Discrete Markov chain)
      • Proposition 1.5
      • Example 1.1 (Phone line)
      • Example 1.2 (Random walk on $\mathbb{Z}$)
      • Proposition 1.6
      • Definition 1.7 (Homogeneous Markov Chain)
      • Definition (initial distribution)
      • Definition 1.8 (Transition kernel)
      • Example 1.3 (Phone line(continued))
      • Example 1.4 (Random walk on $\mathbb{Z}$ (continued))
      • Definition 1.9 (m-step transition kernel)
      • Proposition 1.10
      • Example 1.5 (phone line(continued))
    • 1.2.2 Classificaton of states
      • Definition 1.11 (classification of states)
      • Example 1.6
      • Definition 1.12 (Irreducibility)
      • Definition 1.13 (Period)
      • Example 1.7 (Ex. 1.6 continued)
      • Proposition 1.14
    • 1.2.3 Recurrence and transience
      • Definition 1.15 (Recurrence and transience)
      • Proposition 1.16
      • Proposition 1.17
      • Example 1.8 (ex.16, 1.7 continued)
    • 1.2.4 Invariant distribution and equilibrium
      • Definition 1.18 (Invariant distribution)
      • Example 1.9 (Phone line (continued))
      • Example 1.10 (Random walk on Z (continued))
      • Theorem 1.19 (convergence to equilibrium)
      • Example 1.11 (Phone line (continued))
      • Example 1.12
    • 1.2.5 Reversibility and detailed balance
      • definition 1.20 (Time-reversed chain)
      • Example 1.13(Photo line(cont.))
      • Definition 1.21 (Detailed balance, 詳細釣り合い条件)
      • Theorem 1.22
      • Example 1.14

Chapter 1. Markov Chains

1.1 Stochastic processes

Markov chain は無記憶性という特別な性質を持つ確率過程の一種である. Markov chainをよく学ぶため,まずはStochastic processの概念を形式的に定義する.

Definition 1.1 (Stochastic process)

$X$がstochastic processである
$\Leftrightarrow X$は $X=\{X_t: t \in T\}$という，$T$を添字集合とした確率変数の集合であって,domainとrange$X_t: \Omega \rightarrow S$は共通である. $t$は”time”(時刻)で,$S$は”state space”(状態空間)と呼ばれる.

$T$には様々な集合が考えられるが,我々が当面扱うのは$T$が離散的な集合である場合(stochastic processes in discrete time)で,例えば$T \subset \mathbb{N}$や$T\subset \mathbb{Z}$の場合である. ほかには$T=[0, \infty)$のような連続時間における過程や$T=\mathbb{R}^2$のような空間的な過程を考えることもある.
また,$S$がどのような集合かも問題で，$S$が離散的な集合であれば($X_t$がr.v.として離散的であれば)，このような過程を離散過程(discrete process)という.

Definition 1.2 (Sample path)

$\omega\in \Omega$について,$\{X_t(\omega); t \in T\}$を$X$の$\omega$におけるsample path という.

$T=\mathbb{N}_0$ならばsample pathは点列で，$T=\mathbb{R}$ならばsample pathは$\mathbb{R} \rightarrow S$なる関数である.
fig.1.1はsample pathの例である.

stochastic processは$X_t$のそれぞれの分布のみではなく,それらの依存関係によっても特徴づけられる. この依存関係の構造はprocessのfinite-dimentional distributionsによって表現できる．すなわち
$$P(X_{t_1}\in A_1, ..., X_{t_k} \in A_k)$$
という具合である. $S\subset \mathbb{R}$であればjoint distibution function(同時分布)は
$$F_{(t_1,...,t_k)}(x_1,...,x_k) = P(X_{t_1}\in (-\infty, x_1], ..., X_{t_k} \in (-\infty, x_k])$$
と記述できる.

$X$が
$$F_{(t_1,...,t_{j-1},t_j,t_{j+1},...,t_k)}(x_1,...,x_{j-1},+\infty,x_{j+1},...,x_k) = F_{t_1,...,t_{j-1},t_{j+1},...,t_k)}(x_1,...,x_{j-1},x_{j+1},...,x_k)$$
を満たすとき,$X$のfinite dimentional distribution functionはconsistentであるという.
stochastic process $X$について,$X$がfinite dimentional distributionsによって完全に記述できるか否かについての部分的な答えが以下の定理である.

Theorem 1.3 (Kolmogorov)

$F_{(t_1,...,t_k)}$はconsistent なfinite-dimensional distribution functionの族とする. このとき,以下を満たすprobability spaceとstochastic process $X$が存在する.
$$F_{(t_1,..,t_k)}(x_1,...,x_k) = P(X_{t_1}\in(-\infty, x_1], ..., X_{t_k} \in (-\infty, x_k])$$

この定理から,あるprocessのfinite-dimensional distributionsを与えれば，そのprocessを特徴づけられる(本当か?). ただし,あるfinite-dimensional distributionsによって特徴づけられる$X$は一意ではない. しかし,そのdistributionsはたかだか可算個のr.v.によるeventの全てに確率を一意に割り当てることができて，この講義の範囲ではそれで十分である.

1.2 Discrete Markov chains

1.2.1 Introduction

この節ではMarkov chainのうち，とくに$|S|= |\mathbb{N}|$であるときを考える. これをdiscrete Markov processと呼ぶことはすでに述べた($|S|<\infty$の時も含むが,深くは考えない). discrete Markvo processでは$S$を$\mathbb{N}$として一般性を失わない.

Definition 1.4 (Discrete Markov chain)

$X$はdiscrete stochastic processで，しかも時間についてもdiscreteとする.
$X$がMarkov chain (with discrete state space)
$\Leftrightarrow P(X_{t+1}=x_{t+1}|X_t=x_t,...,X_0=x_0) = P(X_{t+1}=x_{t+1}|X_t=x_t)$
またこの性質をMarkov propertyという.

この定義は，ある時刻における状態がその直前の時刻における状態のみによって決まる(確率的)ということの定式化である.

Proposition 1.5

Markov propertyが成立する$\Leftrightarrow$
任意の$k \in \mathbb{N}$と$t_1<...<t_k\leq t$について
$$P(X_{t+1}=x_{t+1}|X_{t_k}=x_{t_k},...,X_{t_1}=x_{t_1})=P(X_{t+1}=x_{t+1}|X_{t_k}=x_{t_k})$$

Example 1.1 (Phone line)

電話線が使われている状態(1とする)と使われていない状態(0とする)があって，毎分この電話線を監監視する過程$\{X_t|t \in \mathbb{N}_0\}$のstochastic processを考える. $\{X_t\}$がMarkov chainと仮定する. すなわち,これまでどれほど長く電話をしていても1分間後にその電話が切れている確率は変わらず，同様にどれほど長く電話がかかってこなかったとしても1分後に電話がかかってきている確率は変わらない仮定する. Markov assumptionは$\{X_t\}$が同じ分布であることを要求しないので,$P(X_{t+1}=1|X_t=0)=0.8$($t$は昼頃), $P(X_{t+1}=1|X_t=0)=0.1$($t$は深夜)というふうに，時刻によって利用のパターンが異なるモデルにも適用できる.

Example 1.2 (Random walk on $\mathbb{Z}$)

$X_0=0$から始まるrandom walkという確率過程を考える. 全ての時刻で，次の時刻に今あるstateにとどまるか,+1進むか,-1進むかを確率的に選ぶ過程である. 現在あるstateにかかわらず,そのstateにとどまる確率を$1-\alpha-\beta$, -1進む確率を$\alpha$, $+1$進む確率を$\beta$とする.$\alpha, \beta \geq 0, \alpha + \beta \leq 1$である.
$X_{t+1}$を,$P(E_t=-1)=\alpha, P(E_t=0)=1-\alpha-\beta, P(E_t=1)=\beta$任意の$t$に成立するr.v. $E_t$によって
$$X_{t+1} = X_t + E_t$$
と記述する.このとき
$$\begin{aligned} &P(X_{t+1}=x_t-1|X_t=x_t)=\alpha, \\ &P(X_{t+1}=x_t|X_t=x_t)=1-\alpha-\beta, \\ &P(X_{t+1}=x_{t}-1|X_t=x_t)=\beta\end{aligned}$$
であることは明らかである.さらに
$$\begin{aligned} &P(X_{t+1}=x_{t+1}|X_t=x_t,...,X_0=x_0) \\ =&P(E_t=x_{t+1}-x_t|E_{t-1}=x_t-x_{t-1},..,E_0=x_1-x_0,X_0=x_o)\\ =&P(E_t=x_{t+1}-x_t) \ \ \ (\because \{E_t\}\text{の独立性})\\ =&P(X_{t+1}=x_{t+1}|X_t=x_t) \end{aligned}$$
が成り立つから,$\{X_t|t\in \mathbb{N}_0$はMarkov chainである.

Markov chainの分布は初期分布$P(X_0=x_0)$によって完全に定まる.さらに*transition probabilitiesを $P(X_{t+1}=x_{t+1}|X_t=x_t)$と定めると,以下の命題が成立する.

Proposition 1.6

discrete Markov chain $\{X_t | t \in \mathbb{N}_0\}$について，
$$P(X_t=x_t,...,X_0=x_0)=P(X_0=x_0)\cdot \prod_{\tau=0}^{t-1} P(X_{\tau+1}=x_{\tau+1}|X_\tau = x_\tau)$$

proof.

$$\begin{aligned} P(X_t=x_t,...,X_0=x_0) =&P(X_0=x_0) \\ &P(X_1=x_1|X_0=x_0)\\ &P(X_2=x_2|X_1=x_1,X_0=x_0) \\ &\cdots \\ &P(X_t=x_t|X_{t-1}=x_{t-1},...,X_0=x_0) \\ &=P(X_0=x_0)\prod_{\tau=0}^{t-1}P(X_{\tau+1}=x_{\tau+1}|X_\tau =x_\tau) \end{aligned}$$

この証明の1つめの等号は全てのr.v.の組に成立するが,2つ目の等号が成立するのはMarkov chain特有である.
Homogeneous Markov chainという更に特別なクラスは$t$によって$X_t$が変化せず，非常に扱いやすい.以下,全てのMarkov chainはhomogeneousとする.

Definition 1.7 (Homogeneous Markov Chain)

Markov chain $\{X_t|t \in \mathbb{N}_0$がhomogeneous
$\Leftrightarrow P(X_{t+1}=j|X_t=i)=p_{ij} \in [0, 1]$という$t$によらない実数$p_{ij}$が,任意の$i,j \in S$に存在する.

Definition (initial distribution)

initial distributionを$\mathbf{\lambda_0}=(P(X_0=i))_{i \in S}$と書く．$\mathbf{K}$と$\mathbf{\lambda_0}$の組によってhomogeneous Markov chainの分布は完全に定まる(後述).

Definition 1.8 (Transition kernel)

$$\mathbf{K}=(k_{ij})_{ij}, \ \ k_{ij}=P(X_{t+1=j}|X_t=i)$$
という行列$\mathbf{K}$をhomogeneous Markov chain $X$のtransition kernelとかtransition matrixという．
$\sum_{j}k_{ij}=\sum_j P(X_{t+1}=j|X_t=i)=P(X_{t+1}\in S|X_t=i)=1$が成立する.

Example 1.3 (Phone line(continued))

$P(X_{t+1}=0|X_t=0)=0.9, \ P(X_{t+1}=1|X_t=0)=0.1$
$P(X_{t+1}=0|X_t=1)=0.3, \ P(X_{t+1}=1|X_t=1)=0.7$
とするとき,transition kernelは
$$\mathbf{K} = \left(\begin{array}{} 0.9 & 0.1 \\ 0.3 & 0.7 \end{array} \right)$$
となる. transition probabilityを有向グラフを使って表現することが有る. Markov graphという. この例のMarkov graphはfig. 1.4のようになる.

Example 1.4 (Random walk on $\mathbb{Z}$ (continued))

前に挙げたrandom walkのhomogeneous Markov chainのtransition kernelは行，列ともに無限大のToeplitz matrix(テープリッツ行列)で，具体的には
$$\left(\begin{array}{} \ddots & \ddots & \ddots & \ddots & \ddots \\ \ddots & \alpha & 1-\alpha-\beta & \beta & 0 \\ \ddots & \ddots & \alpha &1-\alpha-\beta & \beta &\ddots \\ \ddots & \ddots & \ddots & \ddots & \ddots \end{array} \right)$$
という形をしている.

Definition 1.9 (m-step transition kernel)

homogeneous Markov chain $\{X_t\}$について,$\mathbf{K}^{(m)} = (k^{(m)}_{ij})_{ij}, \ k^{(m)}_{ij} = P(X_{t+m}=j|X_t=i)$を$m$-step transition kernelという.

Proposition 1.10

$X$をhomogeneous Markov chainとすると,
i. $\mathbf{K}^{(m)} = \mathbf{K}^m$
ii. $P(X_m=j) = (\mathbf{\lambda_0}^T \mathbf{K}^{(m)})_j$
が成立する.

proof.

i. $\mathbf{K}^{(m_1+m_2)} = \mathbf{K}^{(m_1)} \cdot \mathbf{K}^{(m_2)}$を示す.
$$\begin{aligned} P(X_{t+m_1+m_2}=k|X_t=i) &=\sum_j P(X_{t+m_1+m_2}=k,X_{t+m_1}=j|X_t=i)\\&=\sum_j P(X_{t+m_1+m_2}=k|X_{t+m_1}=j,X_t=i)P(X_{t+m_1}=j|X_t=i) \\ &= \sum_j P(X_{t+m_2}=k|X_t=j)P(X_{t+m_1}=j|X_t=i) \\&=\sum_j \mathbf{K}_{ij}^{(m_1)} \mathbf{K}_{jk}^{(m_2)} = (\mathbf{K}^{(m_1)} \mathbf{K}^{(m_2)})_{i,k} \end{aligned}$$
ゆえに$\mathbf{K}^{(2)}=\mathbf{K}^2$がたしかに成立し,帰納法により$\mathbf{K}^{(m)}=\mathbf{K}^m$である.
ii. $$P(X_m=j) = \sum_i P(X_m=j,X_0=i) =\sum_i P(X_m=j|X_0=i)P(X_0=i)=(\lambda_0^T\mathbf{K}^m)_j$$

Example 1.5 (phone line(continued))

$$\mathbf{K}=\left(\begin{array}{} 0.9 & 0.1 \\ 0.3 & 0.7 \end{array} \right)$$
のm-step transition kernelは
$$\mathbf{K}^{(m)} = \mathbf{K}^m = \frac{1}{4}\left(\begin{array}{} 3+(\frac{3}{5})^m& 1-(\frac{3}{5})^m \\ 1+(\frac{3}{5})^m & 3-(\frac{3}{5})^m \end{array} \right)$$

1.2.2 Classificaton of states

Definition 1.11 (classification of states)

(a) ある$m \geq 0$があって,$k^{(m)}_{ij} = P(X_{t+m}=j|X_t=i)>0$が成り立つとき$i$は$j$を導くといい,$i \leadsto j$と書く.
(b) $i \leadsto j,\ j \leadsto i$であるとき$i,j$はcommunicateするといい,$i \sim j$と書く.このとき$\cdot \sim \cdot$は同値関係である.

$S$に同値関係を入れられたから,$S$は$\sim$で同値類別できる. ある同値類$C$の全ての元で他の$S$の元に出るpathが無いとき($\Leftrightarrow [\forall i \in C i \leadsto j \Rightarrow j \in C]$), $C$はclosedであるという. ある同値類の元は様々な性質を共有していることをあとで示す.

Example 1.6

$$\mathbf{K} = \left(\begin{array}{}\frac{1}{2} & \frac{1}{4} & 0 & \frac{1}{4} & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & \frac{3}{4} & 0 & 0 & \frac{1}{4} \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & \frac{3}{4} & 0 & 0 & 0 & \frac{1}{4} \\ 0 & 0 & \frac{1}{2} & 0 & 0 & \frac{1}{2} \end{array} \right)$$
をtransition kernelにもつMarkov chainのMarkov graphはfig. 1.5のとおりである. $2 \sim 4, 2 \sim 5, 3 \sim 6, 4\sim 5$が見て取れて,同値類別は$S\backslash \sim = \{\{1\},\{2,4,5\},\{3,6\}\}$であり，特に$\{3, 6\}$のみがclosedである.

Definition 1.12 (Irreducibility)

$S$の全ての元が互いにcommunicateするとき$S$はirreducibleであるという.

phone lineの例とrandom walkの例はreducibleであり,example 1.6はirreducibleである. ところで,ex. 1.6において,$2$から再び$2$に至るのは,$2$から$4$,$4$から$5$,$5$から$2$というステップを踏まなければならないので,$P(X_t=2|X_\tau=2) = \begin{cases} = 0 \ \ &(t \notin \tau + 3\mathbb{N}) \\ >0 &\text{otherwise} \end{cases}$
である. このような振る舞いをperiodicityという.

Definition 1.13 (Period)

(a) $i \in S$のperiod $d(i)$を
$$d(i) = gcd\{m\geq 1| K^{(m)}_{ii} > 0\}$$
と定める. periodを持たないstateも存在する.
(b) $d(i) = 1$なら,$i$はaperiodicという.
(c) $d(i) > 1$なら,$i$はperiodicという.

Example 1.7 (Ex. 1.6 continued)

すでに述べたように$d(2)=3$である. 同様に$d(4)=d(5)=3$である.
また$\mathbf{K}_{1,1},\mathbf{K}_{3,3}, \mathbf{K}_{6,6}>0$だから,$d(1)=d(3)=d(6)=1$.すなわち$3,6$はaperiodicである.
ある同値類(以後,communicateによる同値類をcommunicating classか,単にclassという)においてperiodを共有しているのは偶然ではない.

Proposition 1.14

(a) あるclassの全ての元はperiodを共有する.
(b) irreducible chainでは全ての元はperiodを共有する.

1.2.3 Recurrence and transience

ex. 1.6のMarkov chainを辿り続けると,そのうち$3, 6$を往復するだけになる. このような振る舞いを定式化するため, number of visits in state $i$:
$$V_i = \sum_{t=1}^\infty 1_{\{X_t =i\}}$$
を導入する. 初期値が$i$であるときの条件付き期待値は
$$E[V_i|X_0=i] = E\left[\sum_t 1_{\{X_t=i\}} | X_0=i \right] = \sum_{t=0}^\infty E[1_{\{X_t=i\}}|X_0=i] = \sum_t P(X_t=i|X_0=i) = \sum_{t} k^{(t)}_{ii}$$
この値が有限か無限かによってstateを分類する.

Definition 1.15 (Recurrence and transience)

(a) $i$がrucurrent
$\Leftrightarrow\ E(V_i|X_0=i)=\infty$
(b) $i$がtransient
$\Leftrightarrow\ E(V_i|X_0=i)<\infty$

すなわち,$i$がa.s.無限回訪れられるというのがrecurrentで，a.s.有限回訪れられるというのがtransientである.
prop. 1.14から,あるcommunication classの元たちはrecurrentであるか否かを共有する.

Proposition 1.16

あるcommunicating classにおいて，全ての元がrecurrentであるか，全ての元がtransientであるかのどちらかが成立する.

proof.

$i \sim j$ならば,$i$から$j$に至る長さ$m_{ij}$のpathがあり,$j$から$i$に一有る長さ$m_{ji}$のpathがある. すなわち$k^{(m_{ij})}_{ij}, k_{ji}^{(m_{ji})}>0$である.
$E(V_i|X_0=i) = \sum_t k^{(t)}_{ii} < \infty$とすると，
$$\begin{aligned} EV(V_j|X_0=j) &=\sum_{t} k^{(t)}_{jj} = \frac{1}{{k_{ij}^{(m_{ij})} k^{(m_{ji})}_{ji}}}\sum_t \underline{k_{ij}^{(m_{ij})}k_{jj}^{(t)}k_{ji}^{(m_{ji})}}_{(1)} \\ &\leq \frac{1}{{k_{ij}^{(m_{ij})} k^{(m_{ji})}_{ji}}} \sum_t k^{(m_{ij}+t+m_{ji})}_{ii} \\ &\leq \frac{1}{{k_{ij}^{(m_{ij})} k^{(m_{ji})}_{ji}}} \sum_{s\geq0} k_{ii}^{(s)} < \infty\end{aligned}$$
((1): $j$から$i$に行って,さらに$t$後にまた$i$に戻って,そこから$j$に戻るという確率)
よって$j$もまたtransientである.
また$i$がrecurrentであるとき,
$$\begin{aligned}E[V_j|X_0=j] = \sum_t k^{(t)}_{jj} &\geq \sum_{t=0}^{m_{jj}+ m_{ij}} k^{(t)}_{jj} + \sum_{\tau \geq 0}k^{(m_{ji})}_{ji} k^{(\tau)}_{ii} k^{(m_{ij})}_{ij} \\ &\geq \sum_{t=0}^{m_{jj}+ m_{ij}} k^{(t)}_{jj} + k^{(m_{ji})}_{ji} k^{(m_{ij})}_{ij}\sum_{\tau \geq 0} k^{(\tau)}_{ii} \geq \infty \end{aligned}$$
から,$j$もrecurrent.

Proposition 1.17

(a) closedでないclassはtransient
(b) 有限かつclosedなclassはrecurrent

Example 1.8 (ex.16, 1.7 continued)

ex.1.6において,$\{1\}, \{2,4,5\}$はclosedでないからtransient.
一方$\{3,6\}$は有限かつclosedなのでrecurrent.

1.2.4 Invariant distribution and equilibrium

invariant distributionを導入して,Markov chainの長期的な振る舞いを調べる.

Definition 1.18 (Invariant distribution)

$\mathbf{\mu} = (\mu_i)_{i \in S}$は$S$上のprobablility distributionとする. また$X$がMarkov chainでtransition kernel $\mathbf{K}$をもつとする. $\mathbf{\mu}$が$X$のinvariant distribution (stationary distibution)
$\Leftrightarrow \mathbf{\mu}^T \mathbf{K} = \mathbf{\mu}^T$

さらにこのとき,右から$\mathbf{K}$を掛けることで,
$$\forall i\in \mathbb{N}. \mu^T = \mu^T \mathbf{K}^{n}=\mu^T \mathbf{K}^{(n)}$$
が成立する.したがってprop. 1.10より
$$P(X_m=j) = (\mu \mathbf{K}^{(m)}_j) = (\mu)_j$$
が任意の$m$に成立する. つまり,$X$の分布は時刻によって変化しない.

Example 1.9 (Phone line (continued))

$$\mathbf{K} = \left(\begin{array}{} 0.9 & 0.1 \\ 0.3 & 0.7 \end{array} \right)$$
のstationary distributionを見つける.
$\mu^T \mathbf{K}=\mu^T$を変形して,$\mathbf{K}^T\mu = \mu$.よって$\mu$は$\mathbf{K}^T$のeigenvector(固有ベクトル)であって,ただし確率の公理から$\mu$のそれぞれの要素は非負で，総和は1である.これを解くと,$\mu=[\frac{3}{4}, \frac{1}{4}]^T$である.
Markov chainは必ずしもstationary distributionを持たない.$\mathbb{Z}$上のrandom walkがその例である.

Example 1.10 (Random walk on Z (continued))

$$\mathbf{K}= \left(\begin{array}{} \ddots & \ddots & \ddots & \ddots & \ddots \\ \ddots & \alpha & 1-\alpha-\beta & \beta & 0 \\ \ddots & \ddots & \alpha &1-\alpha-\beta & \beta &\ddots \\ \ddots & \ddots & \ddots & \ddots & \ddots \end{array} \right)$$
であることはすでに言った.$\mu = [1, 1,...]^T$は$\mu^T \mathbf{K}=\mu^T$の唯一の解だが,$\mu$は無限次元のベクトルなので正規化できない.

ある種のMarkov chainは長期的にはstationary distributionに至る.

Theorem 1.19 (convergence to equilibrium)

$X$がirreducibleかつaperiodicなMarkov chainで,stationary distribution $\mu$をもつとする.このとき
$$\lim_{t \rightarrow \infty}P(X_t = i) \rightarrow \mu_i$$
が任意の$i \in S$に成立する.

proof. (sketch)

$X$のinitial distribution, transition kernelをそれぞれ$\lambda$, $\mathbf{K}$とする. initial distribution $\mu$(stational)とtransition matrix $\mathbf{K}$をもつ新しいMarkov chain $Y$を定める. また$T$を$X, Y$が初めて$i\in S$に同時に到達する時刻の確率変数を$T$とする.すなわち
$$T = \min \{t\geq 0| X_t=Y_t=i\}$$
さらに$P(T<\infty)=1$であり，また新しいprocess $Z$を
$$Z_t = \begin{cases} X_t \ \ &(t \leq T) \\ Y_t &(t > T)\end{cases}$$
によって定める. $Z$の概略はfig.1.6のようになる. $Z$はinitial distribution $\lambda$をもち,transition kernel $\mathbf{K}$である. したがって$X$と$Z$は常に同じ分布を持つ.すなわち$\forall j. \ P(X_t=j)=P(Z_t=j)$である.
$Y$のinitial distributionはstationaryなので,$\forall t\forall j.\ P(Y_t = j) = \mu_j$である. $t \rightarrow \infty$において$P(\{Y_t = Z_t\})=1$であって,ゆえに
$$P(X_t=j)=P(Z_t=j)\rightarrow P(Y_t=j)=\mu_j$$

Example 1.11 (Phone line (continued))

phone lineの例で,$\mu=(3/4, 1/4)$だから,$P(X_t=0) \rightarrow 3/4, P(X_t=1) \rightarrow 1/4$.

Example 1.12

Theorem 1.19におけるaperiodicityの仮定が必須であることを示す.
$S=\{1, 2\}$, $\mathbf{K} = (\begin{array}{} 0 & 1 \\ 1 & 0\end{array})$とする. これは明らかにirreducibleだがperiod 2である. stationary distibutionは$(1/2, 1/2)$だが,これは決定論的な過程だから,明らかに$X$はこれに収束しない.

1.2.5 Reversibility and detailed balance

$P(X_{t+1}|X_t)$のように,現在(あるいは過去)を条件にした未来の状態の確率を論じてきたが，今度は逆に未来の状態を条件にした過去あるいは現在の状態の確率を議論する.
$$P(X_t=j|X_{t+1}=i) = P(X_{t+1}=i|X_t=j)\cdot \frac{P(X_t=j)}{P(X_{t+1}=i)}$$
のように条件付き確率の前後が交換できるのはBayesの定理の教えるところである.
chainを逆に辿っていくような新しいMarkov chainの定義を動機づける.

definition 1.20 (Time-reversed chain)

$\tau \in \mathbb{N}$について,$\{X_t | t = 0,...,\tau\}$をMarkov chainとする. $\{Y_t:t=0,..,\tau\}$を$Y_t = X_{\tau - t}$とすると,$Y$を$X$のtime-reversed chainという.
$$\begin{aligned} P(Y_t = j|Y_{t-1}=i)&=P(X_{\tau-t}=j|X_{\tau-t+1}=i)=P(X_s=j|X_{s+1}=i)\\ &=P(X_{s+1}=i|X_s=j)\cdot \frac{P(X_s=j)}{P(X_{s+1}=i)} = k_{ji}\frac{P(X_s=j)}{P(X_{s+1}=i)} \end{aligned}$$
である.

$X$がhomogeneousでも$Y$がhomogeneousとは限らない.
しかし,$X$のinitial distributionがstational $\mu$であれば,$P(X_{s+1}=i)=\mu_i, P(X_s=j)=\mu_j$が任意の$s$に実数として決まって,
$$P(Y_t=j|X_{t-1}=i)=k_{ji}\frac{\mu_j}{\mu_i} \ \ \ \ \ \ (1.2)$$
が成立するから,$\{Y_t\}$はhomogeneousである.

Example 1.13(Photo line(cont.))

すでに挙げた例で,
$\mathbf{K} = (\begin{array}{} 0.9 & 0.1 \\ 0.3 & 0.7 \end{array}), \mu = (3/4, 1/4)^T$だから, 式(1.2)により
$$\begin{aligned} &P(Y_t=0|Y_{t-1}=0) = k_{00}\frac{\mu_0}{\mu_0}=k_{00}=P(X_t=0|X_{t-1}=]1) \\ &P(Y_t=0|Y_{t-1}=1)=k_{01}\frac{\mu_0}{\mu_1} =0.3 =k_{10}=P(X_t=0|X_{t-1}=1) \\ &P(Y_t=1|Y_{t-1}=0)=k_{10}\frac{\mu_1}{\mu_0}=0.1 = k_{01}=P(X_t=1|X_{t-1}=0) \\ &P(Y_t=1|Y_{t-1}=1) =k_{11}\frac{\mu_1}{\mu_1}=k_{11}=P(X_t=1|X_{t-1}=1) \end{aligned}$$
以上よりこの例では$X$と$Y$は同じtransition probabilityをもつ. このようなchainはtime-reversibleであるという.

time-reversible であるか否かを判別する基準を導入する.

Definition 1.21 (Detailed balance, 詳細釣り合い条件)

transition kernel $\mathbf{K}$がdistribution $\mu$によってdetailed balanceを満たす
$\Leftrightarrow \forall i,j\in S. \ \ \mu_ik_{ij}=\mu_jk_{ji}$

detailed balanceは後で学ぶMarkov Chain Monte Carlo(MCMC)の議論でも極めて重要な役割を果たす. detailed balanceを満たす$\mu$はstationary distributionであり，これはdef. 1.19の定義よりも扱いやすいことが多い.

Theorem 1.22

$X$はMarkov chainで,そのtransition matrix を$\mathbf{K}$はdetailed balanceを$\mu$によって満たすとする.このとき
(i) $\mu$は$X$のstationary distributionである.
(ii) $\mu$がinitial distributionであれば,$X$はtime-reverisbleである．

proof.

(i)仮定より
$$(\mu^T K)_i = \sum_j \mu_j k_{ji} =_{(2)}\sum_j \mu_i k_{ij}= \mu_i \sum_j k_{ij} =_{(1)}\mu_i$$
((1): distributionの定義 (2): detailed balanceの定義)
よって確かに$\mu^T\mathbf{K}=\mu^T$だから,$\mu$はstationary distributionである.
(ii) $Y$を$X$のtime-reversalとする.
$$P(Y_t=j|Y_{t-1}=i)=_{(1)}k_{ji}\frac{\mu_j}{\mu_i}=_{(2)}\frac{k_{ij}\mu_i}{\mu_i}=k_{ij}=P(X_t=j|X_{t-1}=i)$$
(1): 式1.2 (2): detailed balanceの定義
よって確かに$X,Y$はtransition matrixを共有する

一方,stationary distributionを持つからと行ってtime-reversibleであるとは限らない.

Example 1.14

$S=\{1, 2, 3\}$.
$$\mathbf{K} = \left( \begin{array}{} 0 & 0.8 & 0.2 \\ 0.2 & 0 & 0.8 \\ 0.8 & 0.2 & 0\end{array} \right)$$
というMarkov chainを考えると,stationary distributionは$\mu = (1/3, 1/3, 1/3)^T$. Markov graphはfig.1.7の通り.

(1.2)からtime-reversed chain $Y$のtransition matrixが
$$\left( \begin{array}{} 0 & 0.2 & 0.8 \\ 0.8 & 0 & 0.2 \\ 0.2 & 0.8 & 0\end{array} \right)$$
と得られるが,これは$\mathbf{K}$と異なる行列である.