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の通り.