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