「Rによるモンテカルロ法入門」読書ノート：3章モンテカルロ積分その2

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

練習問題 3.1 正規・コーシー-ベイズ推定量のモンテカルロ積分を2通りの方法で計算する

正規・コーシー-ベイズ推定量について， $x = 0,2,4$ だとしてモンテカルロ積分を計算する．

$\delta(x) = \frac{\int_{-\infty}^{\infty} \frac{\theta}{1+\theta^2}e^{-(x-\theta)^2/2} d\theta}{\int_{-\infty}^{\infty} \frac{1}{1+\theta^2}e^{-(x-\theta)^2/2} d\theta}$

さて，この式だけだと何のことか全然分からないのだが，式中で消えている係数をしっかり書くと，

分子： $\int_{-\infty}^{\infty}\theta\times\frac{1}{\pi(1+\theta^2)}\frac{1}{\sqrt{2\pi}}e^{-(x-\theta)^2/2}d\theta$
分母： $\int_{-\infty}^{\infty}\frac{1}{\pi(1+\theta^2)}\frac{1}{\sqrt{2\pi}}e^{-(x-\theta)^2/2}d\theta$

に分解され，結局はコーシー分布と正規分布の掛け算をしていることになる．模式的に書くと以下のようになる．

$\delta(x) = \frac{\int_{-\infty}^{\infty} \theta \times Cauchy \times Normal}{\int_{-\infty}^{\infty} Cauchy \times Normal}$

このように，上と下でそれぞれ積分しているので，モンテカルロ積分の実装においても上下別々に推定値を計算した後に $\delta(x)$ の推定値を求める．また，コーシー分布と正規分布を掛けた形となっているため，どちらの分布からも乱数を生成してモンテカルロ積分ができる．

ちなみに，本問題は4章の例4.2（P.107）にも同様の推定量を扱った問題が出てくる．

a.1 被積分関数をプロットする

左から順に，

$\delta(x)$ の分子の分布
$\delta(x)$ の分母の分布

を， $x = 0,2,4$ についてプロットしたのが以下の図になる．

delta_num <- function(t){t/(1+t^2)*exp(-(x-t)^2/2)}
delta_den <- function(t){1/(1+t^2)*exp(-(x-t)^2/2)}

par(mfrow=c(3,2))
for(x in c(0,2,4)){
  curve(delta_num, from=-10, to=10, main=paste("numerator : x=",x))
  curve(delta_den, from=-10, to=10, main=paste("denominator : x=",x))
}

a.2 コーシー・シミュレーションにもとづくモンテカルロ積分を計算する

コーシー分布から乱数を生成してモンテカルロ積分をシミュレーションする．イテレーションごとの推定値と標準分散のプロットは例3.3のプロットと同じで，黒の線が推定値の推移，上下の黄色の線が推定値の誤差の範囲となっている．

モンテカルロ積分の結果としては

x=0のとき：-0.01243
x=2のとき：1.282
x=4のとき：3.389

となった．

Nsim <- 10^4
par(mfrow=c(1,3))
for(x in c(0,2,4)){
  theta <- rcauchy(Nsim)
  num <- theta * dnorm(theta, mean=x)
  den <- dnorm(theta, mean=x)

  n <- num[den != 0]
  d <- den[den != 0] # ゼロ除算を起こさないようにしている
  y <- n / d
  estint <- (cumsum(n)/(1:length(n))) / (cumsum(d)/(1:length(d)))
  esterr <- sqrt(cumsum((y - estint)^2))/(1:length(y))
  plot(estint, xlab="Iterations", main=paste("x = ",x,", Estimate= ",mean(estint)),
       type="l", lwd=2, xlim=c(0,10^4), ylim=c(x-1,x+1))
  lines(estint+2*esterr, col="gold", lwd=2)
  lines(estint-2*esterr, col="gold", lwd=2)
}

b. 収束する様子を推定値の標準誤差でモニタリングする

上図の黄色の幅で示してある．

c. 正規近似によるモンテカルロ積分を計算する

以下のグラフでは，コーシー乱数のシミュレーションのグラフを重ねてプロットすることで，コーシー乱数を使った場合と正規乱数を使った場合の標準誤差を比較をしている．正規乱数を使ったシミュレーションの結果は，赤とピンクの線で表している．このプロットを見ると，正規乱数を使った場合の方がより誤差が少なく推定値を求められることがわかる．

黒＆黄色：コーシー乱数（問題 a.）
赤＆ピンク：正規乱数（問題 c.）

Nsim <- 10^4
par(mfrow=c(1,3))
for(x in c(0,2,4)){
  theta <- rnorm(Nsim, mean=x)
  num <- theta * dcauchy(theta)
  den <- dcauchy(theta)

  n <- num[den != 0]
  d <- den[den != 0]
  y <- n / d
  estint <- (cumsum(n)/(1:length(n))) / (cumsum(d)/(1:length(d)))
  esterr <- sqrt(cumsum((y - estint)^2))/(1:length(y))
  plot(estint, xlab="Iterations", main=paste("x = ",x,", Estimate= ",mean(estint)),
       type="l", lwd=2, col="red", xlim=c(0,10^4), ylim=c(x-1,x+1))
  lines(estint+2*esterr, col="pink", lwd=2)
  lines(estint-2*esterr, col="pink", lwd=2)
}

練習問題 3.9

この問題は，練習問題3.1の続き．

a.コーシー候補にもとづく受理・棄却アルゴリズムから推定量を計算する

$X \sim \mathcal{N}(\theta,1)$ $\theta \sim \mathcal{C}(0,1)$ において，事後分布 $p(\theta|x)$ は以下のように表される．

$p(\theta|x)=\frac{p(x|\theta)p(\theta)}{p(x)}$

ここで， $p(x|\theta) = \frac{1}{\sqrt{2\pi}}e^{-\frac{(x-\theta)^2}{2}}$ より，

$\begin{align} p(\theta|x) =& \frac{p(\theta)}{p(x)}\frac{1}{\sqrt{2\pi}}e^{-\frac{(x-\theta)^2}{2}} \ \propto& p(\theta)e^{-\frac{(x-\theta)^2}{2}} \end{align}$

となる．受理・棄却アルゴリズムにおける目標分布と候補分布は

目標分布 $f(\theta)$ ： $p(\theta)e^{-(x-\theta)^2/2}$
候補分布 $g(\theta)$ ： $\frac{1}{\pi(1+\theta^2)}$
受理・棄却条件： $U\times M \leq \frac{f(Y)}{g(Y)}$ ，（ $Y \sim \mathcal{C}(0,1)$ ， $U \sim \mathcal{U}_{[0,1]}$ ）

のようになる．受理・棄却アルゴリズムの実験は，練習問題3.1と同様に $x = 0,2,4$ で実行している．以下の図は，受理された乱数のヒストグラムと，その乱数の平均値の値を赤い線で示してある．

Nsim <- 10^4
x <- rep(0, Nsim)
M <- pi
f <- function(theta){ exp(-(xx-theta)^2/2)/(1+theta^2)}

par(mfrow=c(3,1))
for(xx in c(0,2,4)){
  for(i in 1:Nsim){
    y <- rcauchy(1)
    while(runif(1)*M > f(y)/dcauchy(y)){
      y <- rcauchy(1)
    }
    x[i] <- y
  }
  hist(x, freq=F, nclass=150, main=paste("x = ",xx))
  abline(v=mean(x), col="red", lwd=2)
}

b. 正規・コーシー-ベイズ推定量のモンテカルロ積分における，分母と分子の乱数の選択

分母と分子で同じ乱数を使う場合と違う乱数を使う場合で，どのような誤差の影響が出るかを比較する．今回は正規乱数を用いたモンテカルロ積分で試してみる．

黒＆黄色：分母と分子で同じ乱数 $\theta_i$ を使う場合
青＆水色：分母と分子で違う乱数 $\theta_i$ を使う場合

以下の図を見ると，xの値が大きい時に違う乱数を使う場合の方が標準誤差が大きい．つまり，分母と分子で別々にシミュレーションを行う際にも，乱数から生成する値を統一することで推定値の分散を押さえることができる．逆に言えば，分母と分子で別々の乱数を使うと，分母と分子どちらかが極端な値を引いてしまったときには推定値が変な値を取ってしまう場合がある．

Nsim <- 10^4
par(mfrow=c(1,3))
for(x in c(0,2,4)){
  theta <- rnorm(Nsim, mean=x)
  theta_2 <- rnorm(Nsim, mean=x)

  num <- theta * dcauchy(theta)
  den <- dcauchy(theta)
  num_2 <- theta * dcauchy(theta)
  den_2 <- dcauchy(theta_2)

  n <- num[den != 0]
  d <- den[den != 0]
  y <- n / d
  n_2 <- num_2[den_2 != 0]
  d_2 <- den_2[den_2 != 0]
  y_2 <- n_2 / d_2

  estint <- (cumsum(n)/(1:length(n))) / (cumsum(d)/(1:length(d)))
  esterr <- sqrt(cumsum((y - estint)^2))/(1:length(y))
  estint_2 <- (cumsum(n_2)/(1:length(n_2))) / (cumsum(d_2)/(1:length(d_2)))
  esterr_2 <- sqrt(cumsum((y_2 - estint_2)^2))/(1:length(y_2))

  plot(estint, xlab="Iterations", ylab="" ,type="l", lwd=2, col="black", xlim=c(0,10^4), ylim=c(x-1,x+1))
  lines(estint+2*esterr, col="gold", lwd=2)
  lines(estint-2*esterr, col="gold", lwd=2)

  par(new=T)
  plot(estint_2, xlab="Iterations", ylab="", type="l", lwd=2, col="blue", xlim=c(0,10^4), ylim=c(x-1,x+1))
  lines(estint_2+2*esterr_2, col="skyblue", lwd=2)
  lines(estint_2-2*esterr_2, col="skyblue", lwd=2)
}