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

内容：乱数を利用して積分を数値的に求める

モンテカルロ積分とは，乱数を利用して単変量や多変数の積分を近似し，解析的に解けない問題を数値的に解くという手法である．これは2章で示したように，任意の分布に従う乱数を理論上無限に生成できるからこそ成り立つ手法であり，本手法はi.i.dを前提とした大数の法則や中心極限定理のような確率論に落としこんで解析している．具体的なモンテカルロ法の利用例としては確率分布を積分といったことが挙げられ，ベイズなどに限らず確率分布を仮定するような統計手法などにおいて必要となる．

例 3.3 モンテカルロ積分の具体例

図の上半分はh(x) = [\cos(50x) + \sin(20x) ]^2 を0から1の間で図示したもの．モンテカルロ積分ではこの曲線の下側の面積を求めることになる． 図の下半分は，実際にモンテカルロ積分のシミュレーションをしたもの．横軸がシミュレーション回数（サンプル数m）で，縦軸が経験平均\bar{h}_m = \frac{1}{m}\sum_{j=1}^{m}h(x_j) を示したもの．黒の線が平均値の推移，上下の黄色の線が推定標準誤差に基づく信頼幅，赤の直線がRのintegrate関数で求めた積分値を示している．この図から，シミュレーション回数を重ねるごとに数値的に求めた値が実際の積分値に収束していく様子がわかる．

1
2
3
4
5
6
7
8
9
10
11

h <- function(x){(cos(50*x) + sin(20*x))^2}
par(mar=c(2,2,1,1), mfrow=c(2,1))
curve(h, xlab="Function", ylab="", lwd=2)
integrate(h,0,1)
x <- h(runif(10^4))
estint <- cumsum(x)/(1:10^4)
esterr <- sqrt(cumsum((x-estint)^2))/(1:10^4)
plot(estint, xlab="Mean and error range", type="l", lwd=2, ylim=mean(x)+20*c(-esterr[10^4], esterr[10^4]))
lines(estint+2*esterr, col="gold", lwd=2)
lines(estint-2*esterr, col="gold", lwd=2)
abline(h=integrate(h,0,1)$value,lwd=2,col="red")

例 3.4 正規分布のモンテカルロ積分における値の精度と効率

「裾」と呼ばれる部分は正規分布において釣鐘型の山の頂点から遠く離れた部分のことを指しており，正規分布の平均から外れるということは起こる確率がきわめて低い部分である．そのため，正規乱数を基にしたモンテカルロ積分においては，山の頂点に近い部分の値は乱数でよく引きあてるので数は十分なのだが，裾の部分の「値が大きい/小さい」値は乱数でなかなか引き当てることができない．精度に大きな影響を与える部分でありながらもシミュレーションを重ねないと十分な数を得ることができないため，古典的なモンテカルロ積分の方法で精度を上げるためには，試行回数を多くしなければならない．

以下の実験では，実際に有効数字4ケタで正確な数値を出すために，108 回ものシミュレーションを行なっている．最終的な表3.1は，行が試行回数，列がx \sim \mathcal{N}(0,1) のときのxに対応している．一番下の行の値がboundという変数のそれぞれ値と同じになっていることから，108 もの試行回数が必要であることがわかる．

1
2
3
4
5
6
7

x <- rnorm(10^8)
bound <- qnorm(c(0.5, 0.75, 0.8, 0.9, 0.95, 0.99, 0.999, 0.9999))
res <- matrix(0, ncol=8, nrow=7)
for(i in 2:8)
  for(j in 1:8)
    res[i-1, j] <- mean(x[1:10^i] < bound[j])
matrix(as.numeric(format(res, digi=4)), ncol=8)

1
2
3
4
5
6
7
8
9

# > matrix(as.numeric(format(res, digi=4)), ncol=8)
# [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]
# [1,] 0.5200 0.7600 0.8000 0.9300 0.9800 1.0000 1.0000 1.0000
# [2,] 0.4760 0.7400 0.7880 0.8980 0.9550 0.9950 1.0000 1.0000
# [3,] 0.4980 0.7477 0.7967 0.8981 0.9480 0.9889 0.9993 0.9999
# [4,] 0.4994 0.7503 0.7996 0.8992 0.9504 0.9899 0.9989 0.9999
# [5,] 0.4999 0.7503 0.8002 0.9001 0.9501 0.9899 0.9990 0.9999
# [6,] 0.4999 0.7500 0.8000 0.9001 0.9500 0.9900 0.9990 0.9999
# [7,] 0.5000 0.7500 0.8000 0.9000 0.9500 0.9900 0.9990 0.9999