当記事は「統計学実践ワークブック(学術図書出版社)」の読解サポートにあたってChapter.32の「シミュレーション」に関して演習問題を中心に解説を行います。モンテカルロ法を用いた乱数生成や、生成した乱数を用いたモンテカルロ積分など、重要なトピックが多いので抑えておくと良いと思います。
Contents
本章のまとめ
下記で取りまとめを行なった。
https://www.hello-statisticians.com/explain-terms-cat/monte_carlo1.html
演習問題解説
例32.1
生成した乱数を$(U_1,V_1),…,(U_N,V_N)$のようにおくと、$M$は二項分布$\displaystyle Bin \left( \frac{\pi}{4},N \right)$に従うので、$M$の分散$V[M]$は下記のように表すことができる。
$$
\large
\begin{align}
V[M] = N \times \frac{\pi}{4} \times \left( 1-\frac{\pi}{4} \right)
\end{align}
$$
よって、$\displaystyle \hat{\pi} = \frac{4M}{N}$の分散$V[\hat{\pi}]$は下記のように計算できる。
$$
\large
\begin{align}
V[\hat{\pi}] &= V \left[ \frac{4M}{N} \right] \\
&= \frac{4^2}{N^2} V[M] \\
&= \frac{4^2}{N^2} \times N \times \frac{\pi}{4} \times \left( 1-\frac{\pi}{4} \right) \\
&= \frac{4\pi}{N} \left( 1-\frac{\pi}{4} \right) \\
&= \frac{\pi(4-\pi)}{N}
\end{align}
$$
ここで$SD[\hat{\pi}]=\sqrt{V[\hat{\pi}]}=0.01$を$N$に関して解く。
$$
\large
\begin{align}
SD[\hat{\pi}] &= 0.01 \\
SD[\hat{\pi}]^2 &= 10^{-4} \\
V[\hat{\pi}] &= 10^{-4} \\
\frac{\pi(4-\pi)}{N} &= 10^{-4} \\
N &= \pi(4-\pi) \times 10^{4} \\
&= 26967.6…
\end{align}
$$
上記に基づいて用いるサンプルの組の数を考えれば良い。
例32.2
確率変数$X$の累積分布関数を$F(x)$とおくと、$F(x)$は下記のように計算できる。
$$
\large
\begin{align}
F(x) &= \int_{0}^{x} 2 e^{-2x} dx \\
&= \left[ -e^{-2x} \right]_{0}^{x} \\
&= 1-e^{-2x}
\end{align}
$$
ここで$u=F(x)$とおき、逆関数$x = F^{-1}(u)$を考える。
$$
\large
\begin{align}
u &= F(x) \\
u &= 1-e^{-2x} \\
e^{-2x} &= 1-u \\
-2x &= \log{(1-u)} \\
x &= -\frac{1}{2} \log{(1-u)}
\end{align}
$$
よって、確率変数$X$に従う乱数を生成するには$U(0,1)$に従う乱数を生成し、$\displaystyle x = -\frac{1}{2} \log{(1-u)}$によって$x$を計算すれば良いことがわかる。
例32.3
採択・棄却法は下記の手順で行うことができる。
1) 一様分布$U(0,1)$に基づいて$x$と$u$を発生させる。
2) $\displaystyle r = \frac{1}{M_i}x(1-x)$を計算し、$u \leq r$ならば$x$を生成し、そうでなければ棄却する。
3) 1)と2)を必要数のサンプルが得られるまで繰り返す。
また、$l(x)=x(1-x)$の最大値が$\displaystyle x = \frac{1}{2}$のときに$\displaystyle \frac{1}{4}$となるので、$M_2$の方が効率良くサンプリングを行うことができる。採択確率は$M_1$のとき$\displaystyle \frac{1}{6}$、$M_2$のとき$\displaystyle \frac{2}{3}$であることも確認できる。
例32.4
下記のような計算を実行することで推定を行うことができる。
import numpy as np
np.random.seed(0)
x = 3.*np.random.rand(100000)+2.
g_hat = np.mean(np.e**(-x))
print("Estimated I: {:.3f}".format(3*g_hat))・実行結果
> print("Estimated I: {:.3f}".format(3*g_hat))
Estimated I: 0.129実際の結果が$\displaystyle \int_{2}^{5} e^{-x} dx = \left[-e^{-x}\right]_{2}^{5}=e^{-2}-e^{-5}=0.128…$であることより、この結果が妥当であることが確認できる。
例32.5
$$
\large
\begin{align}
g_1(x) &= e^x \\
g_2(x) &= g_1(x) – h(x) + \int_{0}^{1} h(x) dx \\
&= e^x – x + \frac{1}{2}
\end{align}
$$
$g_1(x), g_2(x)$は上記のように表すことができる。また、$I_1=I_2=e-1$が成立する。よって、$\sigma_1^2, \sigma_2^2$は下記のようになる。
$$
\large
\begin{align}
\sigma_1^2 &= \int_{0}^{1} (g_1(x)-I_1)^2 dx \\
&= \int_{0}^{1} (e^x-(e-1))^2 dx \\
&= \int_{0}^{1} (e^x-e+1)^2 dx = 0.2420 \\
\sigma_2^2 &= \int_{0}^{1} (g_2(x)-I_2)^2 dx \\
&= \int_{0}^{1} \left( e^x – x + \frac{1}{2} – (e-1) \right)^2 dx \\
&= \int_{0}^{1} \left( e^x – e – x + \frac{3}{2} \right)^2 dx = 0.0437
\end{align}
$$
これより$\hat{I}_2$の推定の方が$\hat{I}_1$の推定よりも精度が良いことがわかる。
例32.6
確率変数が$X \sim U(0,1)$であるので、$g_{1}(X), g_{2}(X)$の期待値は$I$であり、分散$\sigma_1^2, \sigma_2^2$は下記のように表すことができる。
$$
\large
\begin{align}
\sigma_1^2 &= \int_{0}^{1} (g_{1}(x)-I)^2 dx \\
\sigma_2^2 &= \int_{0}^{1} (g_{2}(x)-I)^2 dx
\end{align}
$$
上記に対し、$g_{1}(X), g_{1}(1-X)$の相関係数を$\rho=\rho[g_{1}(X),g_{1}(1-X)]$のようにおき、$\sigma_2^2$を$\sigma_1^2$を用いて表せるように式変形を行う。
$$
\large
\begin{align}
\sigma_1^2 &= \int_{0}^{1} (g_{2}(x)-I) dx \\
&= \int_{0}^{1} \left( \frac{g_{1}(x)+g_{1}(1-x)}{2} – I \right)^2 dx \\
&= \int_{0}^{1} \left( \frac{g_{1}(x)-I}{2} + \frac{g_{1}(1-x)-I}{2} \right)^2 dx \\
&= \int_{0}^{1} \left( \frac{g_{1}(x)-I}{2} \right)^2 + \left( \frac{g_{1}(1-x)-I}{2} \right)^2 + \frac{(g_{1}(x)-I)(g_{1}(1-x)-I)}{4} dx \\
&= \frac{\sigma_1^2}{2}(1+\rho)
\end{align}
$$
ここで$\sigma_1^2=0.242$であり、$\rho$は下記のように計算することができる。
$$
\large
\begin{align}
\rho &= \frac{1}{\sqrt{\sigma_1^2 \times \sigma_1^2}} \int_{0}^{1} (g_{1}(x)-I)(g_{1}(1-x)-I) dx \\
&= \frac{1}{\sigma_1^2} \int_{0}^{1} (e^{x}-I)(e^{1-x}-I) dx \\
&= \frac{-0.234}{0.242} = -0.9669…
\end{align}
$$
よって、推定値$\hat{I}_1, \hat{I}_2$の標準偏差$se(\hat{I}_1), se(\hat{I}_2)$は下記のように表すことができる。
$$
\large
\begin{align}
se(\hat{I}_1) &= \sqrt{V \left[ \frac{1}{2m} \sum_{i=1}^{2m} g_1(X_i) \right]} \\
&= \sqrt{\frac{1}{(2m)^2} \times 2m V \left[ g_1(X) \right]} \\
&= \sqrt{\frac{\sigma_1^2}{2m}} = \frac{\sigma_1}{\sqrt{2m}} \\
se(\hat{I}_2) &= \sqrt{V \left[ \frac{1}{m} \sum_{i=1}^{m} g_2(X_i) \right]} \\
&= \sqrt{\frac{1}{m^2} \times m V \left[ g_2(X) \right]} \\
&= \sqrt{\frac{1}{m} \frac{\sigma_1^2}{2}(1+\rho)} \\
&= \frac{\sigma_1}{\sqrt{2m}}\sqrt{1+\rho} \\
&= \frac{\sigma_1}{\sqrt{2m}}\sqrt{1+\rho} \\
&= \frac{\sqrt{0.033}\sigma_1}{\sqrt{2m}} < \frac{\sigma_1}{\sqrt{2m}} = se(\hat{I}_1)
\end{align}
$$
上記より$\hat{I}_2$の近似の方が精度が良いと考えることができる。
・追記
ワークブックの解答では$se(\hat{I}_2)$が$\displaystyle \frac{\sigma_1}{\sqrt{2m}}(1+\rho)$とされるが、正確には$\displaystyle \frac{\sigma_1}{\sqrt{2m}}\sqrt{1+\rho}$の間違いであると思われる。
例32.7
下記を実行することで計算を行うことができる。
import numpy as np
sample_mean = 11.083
sample_variance = 13.538
bias = 0.021
se_B = 1.010
m = 12.
print("[1] : {:.3f}".format(sample_mean-bias))
print("[2] : {:.3f}".format(se_B/np.sqrt(sample_variance/m)))・実行結果
[1] : 11.062
[2] : 0.951例32.8
$[1]$
下記を実行することで計算することができる。
import numpy as np
x = np.array([11., 9., 9., 12., 5., 8., 7., 12., 13., 14., 18., 15.])
x_sub = np.zeros(x.shape)
for i in range(x.shape[0]):
a = 0
for j in range(x.shape[0]):
if i != j:
a += x[j]
x_sub[i] = a/(x.shape[0]-1)
mean_x_sub = np.mean(x_sub)
print("mean-theta_j: {:.3f}".format(mean_x_sub))・実行結果
> print("mean of theta_j: {:.3f}".format(mean_x_sub))
mean-theta_j: 11.083$[2]$
下記を実行することで計算することができる。
import numpy as np
s2 = 0.10256
se_jack = np.sqrt(11*s2)
print("Standard error of jackknife estimator: {:.3f}".format(se_jack))・実行結果
> print("Standard error of jackknife estimator: {:.3f}".format(se_jack))
Standard error of jackknife estimator: 1.062問32.1
$[1]$
$\hat{\pi}$の標準偏差を$SD[\hat{\pi}]$と定義すると、標準偏差の定義により$SD[\hat{\pi}]$は下記のように考えることができる。
$$
\large
\begin{align}
SD[\hat{\pi}] &= \sqrt{V[\hat{\pi}]} = \sqrt{V \left[ 4 \times \frac{1}{n} \sum_{i=1}^{n} \sqrt{1-U_i^2} \right]} \\
&= \sqrt{\frac{4^2}{n^2} \times n V \left[ \sqrt{1-U^2} \right]} \\
&= \frac{4}{\sqrt{n}} \sqrt{ V \left[ \sqrt{1-U^2} \right] }
\end{align}
$$
ここで$\displaystyle V \left[ \sqrt{1-U^2} \right]$は下記のように計算することができる。
$$
\large
\begin{align}
V \left[ \sqrt{1-U^2} \right] &= E \left[ \sqrt{1-U^2}^2 \right] – E \left[ \sqrt{1-U^2} \right]^2 \\
&= \int_{0}^{1} (1-u^2) du – \left( \frac{\pi}{4} \right)^2 \\
&= \frac{2}{3} – \left( \frac{\pi}{4} \right)^2 = 0.04981…
\end{align}
$$
よって、$SD[\hat{\pi}]=0.01$となるのに必要な$n$は下記のように計算できる。
$$
\large
\begin{align}
SD[\hat{\pi}] &= 0.01 \\
\frac{4}{\sqrt{n}} \sqrt{ V \left[ \sqrt{1-U^2} \right] } &= 0.01 \\
\frac{4}{\sqrt{n}} \sqrt{ \frac{2}{3} – \left( \frac{\pi}{4} \right)^2 } &= 0.01 \\
n &= \frac{4^2}{0.01^2} \left( \frac{2}{3} – \left( \frac{\pi}{4} \right)^2 \right) \\
&= 7970.6…
\end{align}
$$
$[2]$
$\hat{\pi}$の標準偏差を$SD[\hat{\pi}]$とおくと、$[1]$と同様に下記のような計算を行うことができる。
$$
\large
\begin{align}
SD[\hat{\pi}] &= \sqrt{V[\hat{\pi}]} = \sqrt{V \left[ 4 \times \frac{1}{n} \sum_{i=1}^{n} \sqrt{1-U_i^2} + U_i – \frac{1}{2} \right]} \\
&= \frac{4}{\sqrt{n}} \sqrt{ V \left[ \sqrt{1-U^2} + U – \frac{1}{2} \right] } \\
&= \frac{4}{\sqrt{n}} \sqrt{ E \left[ \left( \sqrt{1-U^2} + U – \frac{1}{2} – \frac{\pi}{4} \right)^2 \right]} \\
&= \frac{4}{\sqrt{n}} \sqrt{ \int_{0}^{1} \left( \sqrt{1-u^2} + u – \frac{1}{2} – \frac{\pi}{4} \right)^2 du } \\
&= \frac{4}{\sqrt{n}} \sqrt{0.0144} \\
\end{align}
$$
よって、$SD[\hat{\pi}]=0.01$となるのに必要な$n$は下記のように計算できる。
$$
\large
\begin{align}
SD[\hat{\pi}] &= 0.01 \\
\frac{4}{\sqrt{n}} \sqrt{0.0144} &= 0.01 \\
n &= \frac{4^2}{0.01^2} \times 0.0144 = 2304
\end{align}
$$
参考
・準$1$級関連まとめ
https://www.hello-statisticians.com/toukeikentei-semi1
[…] ・自然科学の統計学 解答・統計学実践ワークブック 解答 […]