「Rによるモンテカルロ法入門」読書ノート：2章受理・棄却法

読書ノートアーカイブ：http://yagays.github.com/blog/2012/10/20/archive-introducing-monte-carlo-methods-with-r/

内容：逆変換法や一般変換法で生成できない分布から乱数を生成する

逆変換法や一般変換法では，一様分布から生成される乱数に何らかの関数を通すことで任意の乱数を生成してきた．しかし，このような直接的な方法で乱数を生成できない分布の場合には，間接的な方法によって乱数を生成する必要がある．そのような場合には，本来求めたい乱数の分布とは違った分布から乱数を生成し，それが受理できる値か棄却できる値かを選り分けることで，間接的に求めたい分布から生成された乱数を表現する．

目標密度（target）: $f$
手段密度・候補密度（candidate）: $g$

この受理・棄却法では1つの乱数の値 $X$ につき，2つの乱数 $Y \sim g$ と $U \sim \mathcal{U}_{[0,1]}$ を考える．前者は候補密度から乱数を生成したもの，後者は受理・棄却に用いる一様乱数．この2つの乱数を，以下の式にあてはめて受理するか棄却するかを判断する．

$U \leq \frac{f(Y)}{M \times g(Y)}$

ここで登場するMは提案分布と候補密度から見積もった定数で，候補密度の制約の中で登場する値．基本的にある値より大きければ何でも良いのだが，シミュレーションの効率に関わってくる．これが小さい値であればあるほど，シミュレーションで棄却される乱数の候補が少なくて済むので，結果として効率が良くなる．最適なMは以下の式で与えられる（UserR!資料pdf p.51）．

$M = \sup_x \frac{f(x)}{g(x)}$

例 2.7, 例2.8 ベータ分布の乱数を生成する(2通り)

黒の点が一様分布/ベータ分布から生成した乱数，赤の点がアルゴリズムで受理された乱数（今回求めたかった乱数），青の線が極限 $M g(x)$ ．

受理・棄却の条件

例2.7: ，，，
- ここでgは何らかの分布ではなく「gは1に等しい」としているので，上式の1つ目の式の中にgが含まれていない．この場合，gで生成される乱数は一様乱数．
- 一様分布のグラフを考えてみるとわかりやすい．区間[0,1]の間で常に値1を取る関数なので，上では1に等しいと表している．
例2.8: $Ug(Y) < f(Y)$ ， $U \sim \mathcal{U}_{[0,M]}$ ， $f \sim \mathcal{Be(2.7,6.3)}$ ， $g \sim \mathcal{Be(2,6)}$

例2.8における「提案分布」という言葉は，新しく設定した候補分布のこと．

# 例2.7
Nsim <- 2500
a <- 2.7
b <- 6.3
M <- 2.67
u <- runif(Nsim, max=M)
y <- runif(Nsim)
x <- y[u < dbeta(y,a,b)]

y2 <- rbeta(Nsim,2,6)
u2 <- runif(Nsim, max=1.67) * dbeta(y2,2,6)
x2 <- y2[u2 < dbeta(y2,2.7,6.3)]

par(mfrow=c(1,2))
plot(y, u, xlim=c(0,1.0), ylim=c(0,2.8), cex=0.5, xlab="y", ylab="u.g(y)")
par(new=T)
plot(x, u[u < dbeta(y,a,b)], col="red", xlim=c(0,1.0), ylim=c(0,2.8), cex=0.5, xlab="y", ylab="u.g(y)")
curve(dbeta(x,a,b), add=T, col="red", lwd=2)
abline(h=M, lwd=2, col="blue")

plot(y2, u2, xlim=c(0,1.0), ylim=c(0,5), cex=0.5, xlab="y", ylab="u.g(y)")
par(new=T)
plot(x2, u2[u2 < dbeta(y2,2.7,6.3)], col="red", xlim=c(0,1.0), ylim=c(0,5), cex=0.5, xlab="y", ylab="u.g(y)")
curve(dbeta(x,2.7,6.3), add=T, col="red", lwd=2)
curve(dbeta(x,2,6)*1.67, add=T, col="blue", lwd=2)

上の図の解釈

乱数の値が赤の分布の線より中に入れば受理となる．候補として生成される乱数は青の分布の線の中でランダムに生成されるので，赤と青の線の間の部分が棄却された乱数になる．この赤と青の線の間の部分が少なければ少ないほど，この受理・棄却法の効率も良くなる．じゃあ効率を良くするにはどうすればいいかというと，いじることのできるポイントは2つ

提案分布gを良い感じにする（青の線の形を決める）
定数Mを出来るだけ小さくする（青の線の高さを決める，fとgの間の幅を決める）

この受理・棄却法の前提として目標分布（赤の線）から乱数は生成できないのだから，100％受理されるようなシミュレーションは理論上不可能．だから，乱数が生成できるような候補分布を使って，定数Mを定めて，間接的に乱数を発生させようということ．

練習問題 2.7 目標分布fと候補密度gがともにベータ分布だった場合に，どのような制約があるのか

$a \leq \alpha, b \leq \beta$ が必要なことの証明

$f(x) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} x^{\alpha-1} (1-x)^{\beta-1}$ ， $g(x) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} x^{a-1} (1-x)^{b-1}$ より

$\frac{f(x)}{g(x)} = \frac{\Gamma(\alpha+\beta)\Gamma(a)\Gamma(b)}{\Gamma(a+b)\Gamma(\alpha)\Gamma(\beta)} x^{\alpha-a} (1-x)^{\beta-b}$

となり，比 $f/g$ が有界になる（どっかのxの値が $\inf$ になったりしない）には，xや(1-x)のベキ数の部分が0以上でである必要がある．よって $a \leq \alpha, b \leq \beta$ が必要．

上の細かい証明は省略するが，要するにベータ分布 $\mathcal{Be}(a,b)$ において，aとbが1未満の時には0と1の両端で値が跳ね上がる．以下の図は，「パターン認識と機械学習上」のP.70を参考にベータ分布を作図したもの．これを見ると，aとbが1のときにちょうど一様分布のような形をとり，それより小さいと下に凸，それより大きいと上に凸のような分布の形になることがわかる．なので，ベータ分布が有界であるにはxや(1-x)のベキ数の部分がポイントになってくる．

# 練習問題 2.7 ベータ関数の分布を図示してみる
par(mfrow=c(2,2))
curve(dbeta(x,0.1,0.1),col="red",lwd=2,ylim=c(0,3),main="Be(0.1,0.1)")
curve(dbeta(x,1,1),col="red",lwd=2,ylim=c(0,3),main="Be(1,1)")
curve(dbeta(x,2,3),col="red",lwd=2,ylim=c(0,3),main="Be(2,3)")
curve(dbeta(x,8,4),col="red",lwd=2,ylim=c(0,3),main="Be(8,4)")

$a = \left\lfloor \alpha \right\rfloor$ ， $b = \left\lfloor \beta \right\rfloor$ の証明

$f(x)/g(x)$ を微分して最大値を取る時のxを求めると， $x_{max} = \frac{\alpha-a}{\alpha-a+\beta-b}$ となる．よって，上で証明した $a \leq \alpha, b \leq \beta$ を満たしつつ自然数aとbが最適な値を取るには， $a = \left\lfloor \alpha \right\rfloor , b = \left\lfloor \beta \right\rfloor$ となる．これもまあ厳密な証明はできなくても直感的に考えれば明らかだろう．

練習問題 2.8 二重指数分布から正規乱数を生成する

$f(x) = \frac{1}{\sqrt{2\pi}}\exp(-\frac{x^2}{2})$
$g(x|\alpha) = \frac{\alpha}{2}}\exp(-\alpha |x|)$

a. 定数Mの取りうる範囲と最小値を求める

$\frac{f(x)}{g(x|\alpha)} = \sqrt{\frac{2}{\pi}}\alpha^{-1}\exp(-\frac{x^2}{2} + \alpha|x|)$

$x = \alpha$ のとき， $\frac{f(\alpha)}{g(\alpha)} = \sqrt{\frac{2}{\pi}}\alpha^{-1}\exp(-\frac{\alpha^2}{2} + \alpha^2) = \sqrt{\frac{2}{\pi}}\alpha^{-1}e^{\alpha^2/2}$ となり，与式を満たす．

Mの最小値は，上式を微分すると $(-\alpha^2+1)\exp(-\frac{\alpha^2}{2})$ が出てくるので， $\alpha = \pm 1$ となる．候補分布に課された制約(P.59)より，fとgどちらか一方だけ負になることは無いので， $\alpha = 1$ が正解．（二階微分すれば-1のとき最小値，1のとき最大値だとわかる）