モンテカルロ法で積分値を推定する問題では、いくつかの分散減少法がしられています。
そのうちの一つである、antithetic variates method (負相関変量法・対称変量法) についてまとめてみようと思います。
では、積分
\[
\int_a^b f(x) \mathrm{d}x=
\int_0^1 f(a+(b-a)r) (b-a)
\mathrm{d}r=\int_0^1 g(r) \mathrm{d}r=
\theta
\] をモンテカルロ法によって推定することにします。
前回見た単純な推定量
\[
\hat{\theta}_{M}=
\frac{1}{N}\left\{\sum_{i=1}^N g(r_i)\right\}
\] は不偏推定量でした。
ただし、$r_i$ は $[0, 1)$ の一様分布に従う乱数とします。
antithetic variates method は、antithetic estimator
\[
\hat{\theta}_{A}=
\frac{1}{N/2}\left\{\sum_{i=1}^{N/2} \frac{g(r_i)+g(1-r_i)}{2}\right\}
\] を用いる方法です。
この推定量は、2 つの不偏推定量の平均になっています。
では、$\hat{\theta}_{M}$ と $\hat{\theta}_{A}$ の分散を比較して、どちらが良いのかを考えてみます。
ここで、$N/2 \in \mathbb{Z}^+$、$\mbox{Var}(g(r_i))=\s^2$ とします。
すると、
\[
\mbox{Var}(\hat{\theta}_{M})=
\frac{1}{N^2}\left\{\sum_{i=1}^N \mbox{Var}(g(r_i)) \right\}=
\frac{\s^2}{N},
\] \[
\begin{eqnarray}
\mbox{Var}(\hat{\theta}_{A}) &=&
\frac{1}{(N/2)^2}\left\{\sum_{i=1}^{N/2} \frac{\mbox{Var}(g(r_i))+
\mbox{Var}(1-g(r_i))+
2\cvar(g(r_i), g(1-r_i))}{2^2}\right\} \\ &=&
\frac{\s^2}{N}+\frac{\sum_{i=1}^{N/2} 2\cvar(g(r_i), g(1-r_i))}{N^2}
\end{eqnarray}
\] となります。
よって、$\cvar(g(r_i), g(1-r_i))$ が負ならば $\hat{\theta}_{A}$ の方が望ましく、正ならば $\hat{\theta}_{M}$ の方が良いという事が分かります。
どういう場合に $\cvar(g(r_i), g(1-r_i))$ が負になるかは、関数 $g(r)$ の形状によって決まります。
$g(r)$ が $r=0.5$ で対称な関数なら $\cvar(g(r_i), g(1-r_i))=0$ であり、$g(r)$ が有界で単調な関数なら $\cvar(g(r_i), g(1-r_i))<0$ となる事が知られています。
もちろん正になってしまう場合もあります。
ものすごく簡単に言うと、$r=0.5$ で対称でなく、凹凸ではない関数の方が $\hat{\theta}_{A}$ の分散が小さいという事になると思います。
では、$g(r)=a r + b$ の場合を考えてみます。
$r$ は一様乱数に従う確率変数なので、
\[
\mbox{Var}(g(r))=\mbox{Var}(a r + b)=\frac{a^2}{12}
\] また、
\[
\begin{eqnarray}
\cvar(g(r), g(1-r))&=&
\cvar(a r, -a r)= a^2 \{\mbox{E}(r)\}^2 -a^2 \mbox{E}(r^2) \\&=&
\frac{a^2}{4}-\frac{a^2}{3} =-\frac{a^2}{12}
\end{eqnarray}
\] であり、
\[
\mbox{Var}(\hat{\theta}_{M})=
\frac{a^2}{12},
\mbox{Var}(\hat{\theta}_{A})=0
\] となります。
つまり、$g(r)$ が線型写像なら、$\hat{\theta}_{A}$ は完璧な推定量になります。
計算コストについては、$\hat{\theta}_{A}$ は乱数の生成コストが $N/2$ になるので、近似精度の設定やコーディングのしかたによっては高速に計算する事ができます。
非積分関数が複雑な場合はどっちが良くなるか分かりませんが、プログラムは簡単なのでちょっと試しにやってみるのはいいかなと思います。
