新年初投稿、ちょっと調べ物をしたのでメモ。
確率変数 $X$ の期待値 $\theta=\mathrm{E}(X)$ について、モンテカルロ法による近似を考える。
何らかの方法で、$X$ と同じ確率分布に従う $n$ 個の独立な乱数列 $X_i$ ($i=1, 2, \ldots, i, \ldots, n$) が得られるとき、
\[
\hat{\theta}_n=
\frac{1}{n}\sum_{i=1}^{n} X_i
\] により、期待値を近似する事ができる。
$n \rightarrow \infty$ のとき、$\hat{\theta_n} \rightarrow \theta$。
この近似の誤差について考えてみる。
ここで、中心極限定理を用いれば、
\[
Y_n = \hat{\theta_n} - \theta \overset{d}{\rightarrow} N(0, \s^2/n)
\] であるから、$n$ が充分大きいとき、$\s_n^2 = \mathrm{Var}(\hat{\theta}_n) = \mathrm{Var}(Y_n)$ より、
\[
\s_n^2=\s^2/n
\] となる。
通常用いるときは、$\s^2$ が未知なので、これを推定量で置き換えた、
\[
\hat{\s}_n^2=
\hat{\s}^2/n=
\left(\frac{1}{n-1}\sum_{i=1}^{n}(X_i-\hat{\theta}_n)^2\right) / n
\] を用いて分散を評価する。
色々と文献を見ると、信頼区間 $[\hat{\theta}_n - k\sqrt{\hat{\s}_n^2}, \hat{\theta}_n + k\sqrt{\hat{\s}_n^2}]$ ($k$ は定数) の幅が、ある一定の値 $\epsilon$ より小さくなるようにアルゴリズムが組まれている事が多い。
中心極限定理を用いると、$k=2$ の時は約95%信頼区間、$k=3$ の時は約99%信頼区間になる。
$k=3$ になっているものが多分多いと思う、ソフトウェアのマニュアルに書いていない事もあるので、ソースコードをちゃんと読めないと困る事も・・・
仮説検定の Type I error rate や Power に関するシミュレーションを行う場合、確率変数 $X$ は 0 (非棄却) と 1 (棄却) をとり、期待値は Type I error rate (Power) の値になる。
したがって、$X$ は、発生確率 $p=$Type I error rate (Power) の二項分布に従う確率変数と見なして良いだろう。
二項分布の仮定を置く事ができるなら、分散は、
\[
\hat{\s}_n^{2(b)}=
\frac{\hat{\theta_n}(1-\hat{\theta}_n)}{n}
\] である。
二項分布の仮定を置かなくても求める事ができる、
\[
\begin{eqnarray}
\hat{\s}_n^2 &=&
\left(\frac{1}{n-1}\sum_{i=1}^{n}(X_i-\hat{\theta}_n)^2\right) / n \\&=&
\frac{1}{n-1}\left(\sum_{i=1}^{n}X_i^2-n\hat{\theta}_n^2\right) / n \\&=&
\frac{1}{n-1}\left(\frac{\sum_{i=1}^{n}X_i^2}{n}-\hat{\theta}_n^2\right)
\end{eqnarray}
\] 確率変数 $X$ は $\{0, 1\}$ を取るから、$X_i^2=X_i$ と置き換えられるので、
\[
\begin{eqnarray}
\hat{\s}_n^2 &=&
\frac{1}{n-1}\left(\hat{\theta}_n-\hat{\theta}_n^2\right) \\&=&
\frac{\hat{\theta}_n(1-\hat{\theta}_n)}{n-1}
\end{eqnarray}
\] $n$ の影響が出る、もしくは中心極限定理が使えないようなオーダでシミュレーションを行う事はまれだと思うので、どっちを使ってもいいんだろう。
でも統一して考える事ができるし、過小になることはないので、今後は $n-1$ で割って評価しようと思った。
モンテカルロ積分の場合もほとんど同じ。
randomized quasi-Monte Carlo methods についても調べたけど面白かった。
