単純モンテカルロ積分なんていう基礎的な事を復習してみました。
自由度 $\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()
追記
念のためですが、棄却法です。
