「Rによるモンテカルロ法入門」読書ノート：2章練習問題

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

練習問題 2.18

a. f(x)とMg(x)をプロットする

定数Mは例2.7のようにoptimize()関数で求めている．以下の図はfとgの分布を図示したもので，黒い曲線fに対して，青い曲線Mg(x)が全体を覆い被さるようになっていることがわかる．

f <- function(x){exp(-x^2/2) * (sin(6*x)^2 + 3*cos(x)^2*sin(4*x)^2 + 1)}
g <- function(x){exp(-x^2/2)/sqrt(2*pi)}
M <- optimize(f=function(x){f(x)/g(x)}, interval=c(0,1),maximum=T)$objective

par(mfrow=c(1,1))
curve(f(x),from=-5,to=5,ylim=c(0,5))
curve(M*g(x),add=T,col="blue")

b. 受理・棄却アルゴリズムを使ってfから2500個の乱数を生成する

ここで標準正規分布からの乱数を使うのは，g(x)が標準正規密度だから（P.59の受理・棄却法のアルゴリズムの1.にある $Y \sim g$ の部分に該当する）．

この問題では2500個の受理された乱数が欲しいので，P.61のコードを真似て作成している．ただし，このコードは少々効率が悪く，決められた数の乱数を作ったらループが終わりというものではない（受理確率を元にして多めに作ってる）ので，ヒストグラムとして図示する際には2500個のみを取り出して作図している．

Nsim <- 2500
x <- NULL
while(length(x) < Nsim){
  y <- rnorm(Nsim*M)
  x <- c(x, y[runif(Nsim*M,0,M) * g(y) < f(y) ])
}
hist(x[1:Nsim],freq=F,nclass=150,main="Nsim = 2500")
curve(f(x)/integrate(f,-Inf,Inf)$value, add=T, col="red")
# curveの中のintegrateに関しては後述(本問のc.)

概ねfの分布に沿った乱数が生成されているようだ．乱数が2500個だけでは少し不明瞭のため，Nsimの数を10⁴ にして再試したのが下の図．

c. シミュレーションから得られた受理率から正規化定数を求める

P.60にあるように，受理確率は基本的に1/Mだけれども，正規化されていない関数に関しては定数CがMに吸収されているので注意が必要となる．今回の場合，シミュレーションから求めた受理確率rを使うことで，以下の式から正規化定数を近似することができる．

$r = \frac{1}{C \times M}$

今回の場合，Mは10.94，rは0.54となったので，Cは0.17程度だと見積もった．

Nsim <- 10^4
x <- NULL # xは受理した乱数
z <- NULL # zは生成した乱数すべて
while(length(x) < Nsim){
  y <- rnorm(Nsim*M)
  x <- c(x, y[runif(Nsim*M,0,M) * g(y) < f(y) ])
  z <- c(z, y)
}
r <- length(x)/length(z)
C <- 1/(M*r)

> M; r; C
[1] 10.94031
[1] 0.5369048
[1] 0.1702445

さて，答え合わせ(?)だが，正規化定数は単純にfの $[-\infty,\infty]$ の積分を求めることによって計算することができる．これにはRのintegrate()関数を使って-infからinfまでを計算すればよい．確率分布の定義から，これが1になるように正規化定数を定めれば良いということで，以下のように計算した結果，0.1696543となった．

> 1/integrate(f,-Inf,Inf)$value
[1] 0.1696543

練習問題 2.19 二重指数分布が候補分布の受理・棄却アルゴリズムにおいて標準正規分布から乱数を生成する際のMとg(x)の最適化

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

$f/g = \frac{\sqrt{2/\pi}}{\alpha} \exp(\alpha|x|-\frac{x^2}{2})$ より，上界Mを取るときのx_maxの値を調べるには，とりあえずxで微分して0になるときの値を求めれば良い．expの部分を微分すると $(\alpha - x)$ みたいな項が出てくるので， $x_{max} = \alpha$ （厳密には $x_{max} = \pm \alpha$ ）．

次に，受理率を最適化するということは1/Mの値を大きくすれば良いので，結局はMの値の最大値を求めれば良い． $x_{max} = \alpha$ のときのMを $\alpha$ に関して微分して0になるときの値ということで計算していくと $\frac{(\alpha^2-1)\exp(\alpha^2/2)}{\alpha^2}$ みたいな項が出てくるので，つまり $\alpha = 1$ の時にMが最大値を取り，受理率が最適化される．

練習問題 2.22 切断正規分布から正規乱数を生成する

a. 切断正規分布から乱数を生成する

library(msm) # dtnorm()のため
mu <- 0
sigma <- 1

par(mfrow=c(4,2))
for(a in c(-2.32, -0.84, -0.67, 0, 0.67, 0.84, 2.32)){
  N <- 10^4
  X <- rep(0,N)
  for(i in 1:N){
    z <- rnorm(1, mean=mu, sd=sigma)
    while(z < a){
    z <- rnorm(1, mean=mu, sd=sigma)
    }
  X[i] <- z
  }
  hist(X, freq=F, nclass=150, main=paste("Truncated Normal Distribution: a = ",a))
  curve(dtnorm(x, mean=mu, sd=sigma, lower=a, upper=Inf), add=T, col="red", lwd=2)
}

b. シミュレーションので受理する確率を求める

上で示したアルゴリズムで考える．z < aならアタリ，それ以外ならハズレというようなときにどれだけ数を撃てばz<aに入るかということを考えると，正規分布のz < aの面積を求めればいいことがわかる．なお，ここで出てくる $\Phi(x)$ という関数は定義されていないが，以下のような累積標準正規分布（本書P.78に出てくる）のことだと思われる．

$\Phi(t) = \int_{-\infty}^{t} \frac{1}{\sqrt{2\pi}}e^{-y^2/2}dy$

この形式になるように，aからの累積正規分布 $F(x|\mu, \sigma^2, a) \sim \int_{a}^{\infty} \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{2\sigma^2}}dx$ を変形していけばよい．

ここで， $v = -\frac{x-\mu}{\sigma}$ とおくと， $dv = -\frac{dx}{\sigma}$ ， $x:a \rightarrow \infty$ のとき， $v:-\infty \rightarrow \frac{\mu-a}{\sigma}$ をとる．よって，

$F(v|\mu,\sigma^2,a) = \int_{\infty}^{\frac{\mu-a}{\sigma}} \frac{1}{\sqrt{2\pi}}e^{-\frac{v^2}{2}}dv = \Phi(\frac{\mu-a}{\sigma})$   となる．

また，aの値が裾にある（値が大きくなる）と， $\Phi(\frac{\mu-a}{\sigma})$ の値が小さくなるので，シミュレーション回数が多くなる．

c. 切断正規分布の乱数を正規分布から生成することを考える

目標分布； $\exp(-(x-\mu)/2\sigma^2)\mathbb{I}_{x \geq a}$
候補分布： $\frac{1}{\sqrt{2\pi}}\exp(-\frac{(x-\bar{\mu})^2}{2})$

より

$\frac{f(x)}{g(x)}=\frac{\exp(-(x-\mu)/2\sigma^2)\mathbb{I}_{x\geq z}}{\frac{1}{\sqrt{2\pi}}\exp(-\frac{(x-\bar{\mu})^2}{2})}\$ $= \sqrt{2\pi}\exp(-x\bar{\mu}+\frac{1}{2}\bar{\mu}^2)\mathbb{I}_{x \geq a}$

となり候補分布の制約を満たすので，正規候補にもとづく受理・棄却アルゴリズムは切断正規分布から乱数を生成するために使える．

次に，このときの最適なMを考える．x=aにおいて $\frac{f(a)}{g(a)} = \sqrt{2\pi}\exp(-a\bar{\mu}+\frac{1}{2}\bar{\mu}^2)\mathbb{I}_{x \geq a}$ より，この式を $\bar{\mu}$ で微分すると

$\Bigl(\frac{f(a)}{g(a)}\Bigr)' = \sqrt{2\pi}(-a+\bar{\mu})\exp(-a\bar{\mu}+\frac{1}{2}\bar{\mu}^2)\mathbb{I}_{x \geq a}$

となるので，Mは $\bar{\mu}=a$ のとき最適値をとる．

d. 候補分布に指数分布を使うことを考える

目標分布； $\exp(-(x-\mu)/2\sigma^2)\mathbb{I}_{x \geq a}$
候補分布： $g_{\alpha}(z)=\alpha e^{-\alpha(z-a)}\mathbb{I}_{x \geq a}$

$\begin{align}\frac{f(x)}{g(x)}&=\frac{\exp(-(x-\mu)/2\sigma^2)\mathbb{I}_{x\geq a}}{\alpha e^{-\alpha(z-a)}}\&=\frac{1}{\alpha}\exp(-\frac{x^2}{2}+\alpha(x-a))\mathbb{I}_{x \geq a}\end{align}$

より，

$\frac{f(x)}{g(x)} \propto e^{\alpha(x-a)}e^{-z^2/2}$

となる．上式を前と同じように微分すると $(-x+\alpha)$ が出てくるので， $x = \alpha$ のとき最大値を取り， $\exp(\alpha^2/2-\alpha a)$ で極限が与えられる．

次に， $a = \alpha$ の場合に正当な候補密度になることを導く．f/gの比が0より大きいことを示せば良いのだが，上式を見ると結局expの外側の $\frac{1}{\alpha}$ ，つまり $\alpha$ を調べないといけない．ということで， $\alpha$ で微分して最適値の時の値を求める．

$x = \alpha$ のとき $\frac{f(\alpha)}{g(\alpha)} = \frac{1}{\alpha}\exp(\frac{\alpha^2}{2}-\alpha a)$ なので，

$\Bigl(\frac{f(\alpha)}{g(\alpha)}\Bigr)' = -\frac{1}{\alpha^2}\exp(\frac{\alpha^2}{2}-\alpha a)+\frac{1}{\alpha}(\alpha-a)\exp(\frac{\alpha^2}{2}-\alpha a)$ $=\frac{1}{\alpha}\bigl(-\frac{1}{\alpha^2}+(\alpha-a)\bigr)\exp(\frac{\alpha^2}{2}-\alpha a)$