読者です 読者をやめる 読者になる 読者になる

単純モンテカルロ積分 (1)

R 統計

単純モンテカルロ積分なんていう基礎的な事を復習してみました。

自由度 $\nu=1$、分散共分散行列が
\[
\mathbf{R}=\begin{pmatrix}1.0&0.7\\0.7&1.0\end{pmatrix}
\] の 2 変量 $t$ 分布について、
\[
\int_{a_1}^{b_1}\int_{a_2}^{b_2}\mathrm{MVT}_2(t_1,\,t_2 \,\mid\, \nu,\, \mathbf{R})\mathrm{d}t_1\mathrm{d}t_2
\] を求めてみます ($a_1\leq t_1 \leq b_1, a_2 \leq t_2 \leq b_2$ で囲まれた領域の確率を求めます)。
図示すると以下の様になります。

library(mvtnorm)

a <- c(-1,-2)
b <- c( 1, 2)

R  <- matrix(c(1,0.7,0.7,1),ncol=2)
N  <- 1000000
h  <- dmvt(c(0,0), sigma=R, log=FALSE)
r1 <- a[1] + (b[1]-a[1])*(runif(N))
r2 <- a[2] + (b[2]-a[2])*(runif(N))
r3 <- h*(runif(N))
r  <- cbind(r1,r2)
ri <- dmvt(r, sigma=R, log=FALSE)
pr <- (b[1]-a[1])*(b[2]-a[2])*h * sum(r3<ri)/N

単純モンテカルロ積分では、まず各軸の範囲や取り得る値に合うように (例えば、$t_1$ 軸方向ならば $[-1, 1]$) の一様乱数を生成します。
r1 は $t_1$、r2 は $t_2$、r3 は "密度" を表わす乱数です。
h は密度が一番高くなる所です、非心分布ではないので、中心の密度が一番高くなります (h は別に何にしても良いですが、取り得る最小の値を与えた方が収束が早いです)。

そして、各点において、乱数 r3 が密度 ri より小さくなる回数をカウントします sum(ri>r3)。
sum(ri>r3) をモンテカルロサンプルの数 N で割ってあげると、生成した領域の中で、乱数が非積分関数の内側に入った割合 (まだ体積じゃない) が求まります。
最後に生成した領域の体積をかけてあげると、その領域の確率 (積分値) が求まります。
\[
V=(b_1-a_1)(b_2-a_2)h\times\frac{1}{N}\sum_{i=1}^N I(r_{1i}, r_{2i}, r_{3i}),
\] \[
I(r_{1i}, r_{2i}, r_{3i})=
\begin{cases}
0, & \qquad r_3 \geq \mathrm{MVT}_2(r_{1i}, r_{2i} \mid\, \nu, \mathbf{R}), \\
1, & \qquad r_{3i} < \mathrm{MVT}_2(r_{1i}, r_{2i} \mid\, \nu, \mathbf{R}),
\end{cases}
\] \[
h=\max_{t_1, t_2} \{\mathrm{MVT}_2(t_1, t_2 \mid \nu, \mathbf{R})\}\quad(a_1\leq t_1 \leq b_1, a_2 \leq t_2 \leq b_2).
\]

答え合わせすると、

pr
pmvt(lower=a,upper=b,df=1,sigma=R)
[1] 0.4609525
[1] 0.4605710

大体同じになります。

単純モンテカルロ積分は効率が悪すぎるので、あまり利用されていません、もっと良いアルゴリズムがあります。
しかし、考え方、特に乱数を生成してからどう処理すると体積が求められるか、は重要なので復習しておいて良かったと思いました。

グラフ用スクリプト

path <- "C:/"
x <- y <- seq(-3, 3, len = 70)

rprod <- function(x, y) {
  r <- matrix(rep(0, length(x)*length(y)), ncol=length(x))
  for (i in 1:length(x)) {
  for (j in 1:length(y)) {
    r[i,j] <- dmvt(x=c(x[i],y[j]), sigma=matrix(c(1,0.7,0.7,1), ncol=2), df=1, log=FALSE)
  }}
  return(r)
}

r<-rprod(x,y)

postscript(
  file=paste(path, "plot.eps", sep=""), width = 9, height = 9,
  pointsize = 20, horizontal = FALSE, onefile = FALSE,
  paper = "special", family="Palatino"
)
par(mar=c(3.2, 3.2, 1.5, 0.5), mgp=c(2.0, 0.7, 0))
par(bty="l")

image(
  x, y, r, col = terrain.colors(240), axes = FALSE,
  breaks=seq(0, 0.24, 0.001),
  xlim=c(-3,3), ylim=c(-3,3),
  xlab=expression(bold(t)[1]),
  ylab=expression(bold(t)[2])
)
contour(
  x, y, r, levels = seq(0, 0.24, by=0.02), add = TRUE, col = "peru",
  xlim=c(-3,3), ylim=c(-3,3)
)
axis(1, at = seq(-3,3,1))
axis(2, at = seq(-3,3,1))
segments(-1, -2,  1, -2, xlim=c(-3,3), ylim=c(-3,3), lty=2)
segments(-1,  2,  1,  2, xlim=c(-3,3), ylim=c(-3,3), lty=2)
segments(-1, -2, -1,  2, xlim=c(-3,3), ylim=c(-3,3), lty=2)
segments( 1, -2,  1,  2, xlim=c(-3,3), ylim=c(-3,3), lty=2)

dev.off()

追記

念のためですが、棄却法です。