今回も
\[
\int_a^b f(x) \mathrm{d}x=
\theta
\] の様な積分値を推定する問題を考えます。
このままでは $[0,\, 1]$ の一様乱数を利用できないので、
\[
x=a+(b-a)r
\] \[
\frac{\mathrm{d}x}{\mathrm{d}r}=b-a
\] \[
0\leq r \leq 1
\] の置換積分を考えます。
これを用いると、
\[
\int_a^b f(x) \mathrm{d}x=
\int_0^1 f(a+(b-a)r) (b-a) \mathrm{d}r
\] となり、
\[
\hat{\theta}=
\frac{1}{N}\left\{\sum_{i=1}^N f(a+(b-a)r_i) (b-a)\right\}
\] の $r_i$ に一様乱数を代入する事で、積分値を推定することができます。
ここで、
\[
\mathrm{E}(\hat{\theta})=
\frac{1}{N}\sum_{i=1}^N \mathrm{E}\left\{ f(a+(b-a)r_i) (b-a)\right\}=\theta
\] \[
\begin{eqnarray}
\mathrm{E}\left\{ f(a+(b-a)r_i) (b-a)\right\}&=&
\int_0^1 f(a+(b-a)r_i) (b-a) \times \mathrm{Uniform}(0,\,1) \mathrm{d}r_i \\&=&
\int_0^1 f(a+(b-a)r_i) (b-a) \times 1 \mathrm{d}r_i \\&=&
\int_a^b f(x) \mathrm{d}x = \theta
\end{eqnarray}
\] となり、多次元でもヤコビアンの絶対値を求めるだけなので、あまり違いはありません。
では、前回 使った例
\[
\int_{-2}^{2}\int_{-1}^{1}\mathrm{MVT}_2(t_1, t_2 \mid \nu, \mathbf{R}) \mathrm{d}t_1\mathrm{d}t_2
\] \[
\nu=1,
\mathbf{R}=
\begin{pmatrix}1.0&0.7\\0.7&1.0\end{pmatrix}
\] を試してみます。
この場合は、変数変換すると以下のようになります。
\[
\int_0^1\int_0^1\mathrm{MVT}_2(a_1+(b_1-a_1)r_{1},
a_2+(b_2-a_2)r_{2}\mid \nu, \mathbf{R})(b_1-a_1)(b_2-a_2) \mathrm{d}r_{1} \mathrm{d}r_{2}
\]
library(mvtnorm) a <- c(-1,-2) b <- c( 1, 2) R <- matrix(c(1,0.7,0.7,1),ncol=2) N <- 1000000 r1 <- a[1] + (b[1]-a[1])*(runif(N)) r2 <- a[2] + (b[2]-a[2])*(runif(N)) r <- cbind(r1,r2) ri <- dmvt(r, sigma=R, log=FALSE) pr <- (b[1]-a[1])*(b[2]-a[2])*sum(ri)/N pr
実行すると、
> pr [1] 0.4604684
という具合に推定できます。
