「Rによるモンテカルロ法入門」読書ノート：3章 重点サンプリング

内容：重点サンプリングにより推定量の分散と効率を改善する

例 3.5 Z \sim \mathcal{N}(0,1) でP(Z > 4.5) の確率を求める

確率を求めると書いてあるが，要するに正規分布の確率密度関数において4.5からInfまでの面積（\int_{4.5}^{\infty}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}} dx ）を求めるというもの．ただし，正規分布でxの値が4.5ともなると確率そのものがものすごく小さいので，普通のモンテカルロ積分で正規分布から乱数を生成したとしても，シミュレーション回数を非常に大きくしないと精度の良い値が求まらない．

1
2
3
4

pnorm(-4.5,log=T)
# [1] -12.59242
pnorm(-4.5)
# [1] 3.397673e-06

ということで，重点サンプリングでシミュレーション回数を抑える．g(y) = \frac{e^{-y}}{\int_{4.5}^{\infty}e^{-x}dx} = e^{-(y-4.5)} として，\frac{1}{m}\sum_{i=1}^{m}\frac{f(Y^{(i)})}{g(Y^{(i)})} を求める．ここで，gの選択の箇所は「gを，4.5で切り詰めた指数分布の密度\mathcal{Exp}(1) とする」ということ(P.81)．

1
2
3
4
5

Nsim <- 10^3
y <- rexp(Nsim)+4.5
weit <- dnorm(y)/dexp(y-4.5)
plot(cumsum(weit)/1:Nsim, type="l")
abline(h=pnorm(-4.5), col="red", lwd=2)

練習問題 3.5 切り詰められた指数分布\mathcal{Exp}(\lambda) をg(x)としたときの近似値の分散への影響を調べる

例3.5で\lambda = 1 だった指数分布の値を色々と変えてみて実験する．以下の図は，\lambda = 1,5,10,20 の時の推定値の変化を示している．左上の図が例3.5と同じ条件となり，そこから右上・左下・右下と\lambda の値を大きくしている．黒の線が推定値の推移，上下の黄色の線が推定値の誤差の範囲（＝分散）となっている． このプロットを見ると，\lambda の値が大きくなるにつれて，分散も大きくなることがわかる．

1
2
3
4
5
6
7
8
9
10
11
12

par(mfrow=c(2,2))
for(lambda in c(1, 5, 10, 20)){
  Nsim <- 10^4
  y <- rexp(Nsim)/lambda+4.5
  weit <- dnorm(y)/dexp(y-4.5, lambda)
  estint <- cumsum(weit)/1:Nsim
  esterr <- sqrt(cumsum((weit-estint)^2))/(1:Nsim)
  plot(estint, xlab="Mean and error range", ylab="prob", type="l", main=paste("lambda = ",lambda))
  lines(estint+2*esterr, col="gold", lwd=2)
  lines(estint-2*esterr, col="gold", lwd=2)
  abline(h=pnorm(-4.5), col="red")
}

また，がくっと値が変わっている部分は乱数でかなりレアな値を引いて全体的に値が引きずられたからだと思われる．weitの中から値の非常に大きいものを10個選び出して，その時の推定値の変化を見たのが以下の図になる．青の点線で示した箇所が，weitの値が非常に大きくなったポイントを示している．

1

abline(v=seq(1,Nsim)[rank(weit)>Nsim-10], col="blue", lty=2)