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$が存在しないためである.