この問題のテーマ・前提知識
- テーマ: モンテカルロ法による円周率 \(\pi\) の推定（二項分布を使う方法・積分近似を使う方法の2通り）
- 前提知識: 二項分布 \(B(n, p)\) の期待値・分散 / 確率変数の線形変換と分散 / 定積分の計算
問題の概要
2つの小問からなる。[1] \([0,1]\) 上の一様乱数ペア \((U, V)\) を \(N\) 組生成し、単位円内に落ちた個数 \(M\) を使って \(\hat{\pi} = 4M/N\) で \(\pi\) を推定する方法について、期待値・分散・標準偏差を求め、SD が 0.01 以下となる \(N\) の下限を求める。[2] \(U\) のみを使って \(\hat{\pi} = \frac{4}{n}\sum\sqrt{1-U_i^2}\) で推定する方法について、同様に分散と \(n\) の下限を求める。
解説
解法の方針：[1] 各試行が二項分布 \(B(N, \pi/4)\) に従うことを利用し、\(\hat{\pi} = 4M/N\) の期待値・分散を線形変換の公式で求める。[2] \(\sqrt{1-U^2}\) の期待値・分散を定積分で計算し、標本平均の分散から \(n\) の下限を導く。
[1]
\(N\) 組の乱数のうち、\(U^2+V^2<1\) となる組の個数 \(M\) は、二項分布 \(B(N, \pi/4)\) に従う。
\[E[M]=N\frac{\pi}{4}\]
\[V[M]=N\frac{\pi}{4}\left(1-\frac{\pi}{4}\right)\]
より、
\[E[\hat{\pi}]=E\left[\frac{4M}{N}\right]=\frac{4}{N}E[M]=\frac{4}{N}\cdot N\frac{\pi}{4}=\pi\]
\[V[\hat{\pi}]=V\left[\frac{4M}{N}\right]=\frac{4^2}{N^2}V[M]=\frac{4^2}{N^2}\cdot N\frac{\pi}{4}\left(1-\frac{\pi}{4}\right)=\frac{4}{N}\pi\left(1-\frac{\pi}{4}\right)\]
とわかる。
従って、\(\hat{\pi}\) の標準偏差 \(SD[\hat{\pi}]\) は、
\[SD[\hat{\pi}]=\sqrt{V[\hat{\pi}]}=\sqrt{\frac{\pi(4-\pi)}{N}}\leq 0.01\]
となり、これを満たす \(N\) は、
\[N\geq\frac{\pi(4-\pi)}{0.01^2}\approx \frac{2.7}{0.0001}=27000\]
[2]
\begin{eqnarray}
V[\hat{\pi}] &=& V\left[4\times\frac{1}{n}\sum_{i=1}^n\sqrt{1-U_i^2}\right]\\
&=& \frac{4^2}{n^2}\sum_{i=1}^nV\left[\sqrt{1-U_i^2}\right] \\
&=& \frac{4^2}{n^2}\sum_{i=1}^n\left\{E[1-U_i^2]-\left(E\left[\sqrt{1-U_i^2}\right]\right)^2\right\} \\
&=& \frac{4^2}{n^2}\sum_{i=1}^n\left\{\int_0^1(1-u^2)du-\left(\frac{\pi}{4}\right)^2\right\} \\
&=& \frac{4^2}{n^2}\sum_{i=1}^n\left\{\left[u-\frac{1}{3}u^3\right]_0^1-\frac{\pi^2}{4^2}\right\} \\
&=& \frac{4^2}{n^2}\sum_{i=1}^n\left\{\left(1-\frac{1}{3}\right)-\frac{\pi^2}{4^2}\right\} \\
&=& \frac{4^2}{n^2}\sum_{i=1}^n\left\{\frac{2}{3}-\frac{\pi^2}{4^2}\right\} \\
&=& \frac{4^2}{n^2}\left(\frac{2}{3}-\frac{\pi^2}{4^2}\right)n \\
&=& \frac{4^2}{n}\left(\frac{2}{3}-\frac{\pi^2}{4^2}\right) \\
\end{eqnarray}
従って、\(\hat{\pi}\) の標準偏差 \(SD[\hat{\pi}]\) は、
\[SD[\hat{\pi}]=\sqrt{\frac{4^2}{n}\left(\frac{2}{3}-\frac{\pi^2}{4^2}\right)}\leq 0.01\]
となり、これを満たす \(n\) は、
\[n\geq\frac{16\times0.05}{0.01^2}\approx 8000\]
よくある質問
Q. [1] で \(M\) が二項分布に従う理由がわかりません。
各試行 \((U_i, V_i)\) が単位円内に入る確率は正方形の面積 1 に対する四分円の面積 \(\pi/4\) であり、各試行は独立。「\(N\) 回の独立なベルヌーイ試行で成功確率が一定」という二項分布の定義をそのまま満たすため、\(M \sim B(N, \pi/4)\) となる。
Q. [2] で \(E[\sqrt{1-U^2}]=\pi/4\) となる理由は？
\(U\) が \([0,1]\) 上の一様分布に従うとき、\(E[\sqrt{1-U^2}]=\int_0^1\sqrt{1-u^2}\,du\) となる。この積分は単位円の四分円の面積に等しいため \(\pi/4\) になる。これが [2] の推定量が不偏推定量である根拠でもある。
Q. [1] と [2] でどちらの方法が精度が高いですか？
[2] の積分近似法の方が同じサンプル数で分散が小さく（\(n \approx 8000\) vs \(N \approx 27000\)）、精度が高い。これは \(V\) を積分で条件付き消去する分散低減法（条件付け法）の一種であり、乱数のばらつきの一部を解析的に取り除けるため、モンテカルロ法の改良手法として重要な考え方。
コメント