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)がある.