応用

確率論入門　第12章 §12.4　★★★ 期待値・積分を乱数の平均で近似する Monte Carlo 法と、標本から分布関数全体を一様に推定する Glivenko–Cantelli の定理

Monte Carlo 積分

$d$ 次元の立方体の上の可積分な関数 $g$ について、次の積分を数値的に求めたいとする。

$$I=\int_{[0,1]^d}g\,\mathrm{d}\mathrm{Leb}$$

$\bm U_1,\bm U_2,\dots$ を、立方体の上の一様分布に従う独立同分布の確率ベクトルとする。第5章 §5.3 の変数変換公式から、$E[g(\bm U_1)]=I$ である。そこで

$$\hat I_n=\frac1n\sum_{k=1}^{n}g(\bm U_k)\tag{12.12}$$

を積分の推定値とする。これが Monte Carlo 積分である。各項 $g(\bm U_k)$ は独立同分布で、$E[\left\vert g(\bm U_1)\right\vert ]=\int\left\vert g\right\vert <\infty$ である。よって Kolmogorov の強法則（定理 12.6）から、$\hat I_n\to I$ a.s. となる。

誤差の大きさは分散で測る。次の分散が有限であるとする。

$$\sigma_g^2=V[g(\bm U_1)]=\int g^2-I^2<\infty$$

このとき推定値の平均は $I$、分散は $\sigma_g^2/n$ であり、

$$\sqrt{E\bigl[(\hat I_n-I)^2\bigr]}=\frac{\sigma_g}{\sqrt n},\qquad P\Bigl(\left\vert \hat I_n-I\right\vert \ge\frac{c\,\sigma_g}{\sqrt n}\Bigr)\le\frac1{c^2}\tag{12.13}$$

となる（2つ目は Chebyshev の不等式）。誤差が $n^{-1/2}$ に比例し、その係数が次元によらないことが Monte Carlo 法の特徴である。

格子点で関数値をとる台形公式と比べる。格子点を $d$ 次元に並べると、1辺あたりの点の数は $n^{1/d}$ になる。よって台形公式の誤差は $n^{-2/d}$ の程度になり、5次元以上では Monte Carlo の $n^{-1/2}$ より遅くなる。

実際の計算では $\sigma_g$ も未知である。そこで同じ乱数から、次の量で推定する。

$$\hat\sigma_n^2=\frac1n\sum_k(g(\bm U_k)-\hat I_n)^2$$

確認テスト 12.1-4 と同じ議論（強法則版）から、$\hat\sigma_n^2\to\sigma_g^2$ a.s. である。(12.13) の確率の評価は、第13章の中心極限定理で精密になる。

$$P\Bigl(\left\vert \hat I_n-I\right\vert \ge\frac{1.96\,\sigma_g}{\sqrt n}\Bigr)\to0.05$$

一様分布でない分布 $f$ に関する期待値 $E[h(X)]$ も、同じ分布の乱数を作れれば同じ方法で求まる。乱数の作り方には、第3章 §3.5 の逆変換法 $F^{-1}(U)$ や、第7章 §7.2 の Box–Muller 法がある。

経験分布関数と Glivenko–Cantelli の定理

$X_1,X_2,\dots$ を分布関数 $F$ の独立同分布の列とする。最初の $n$ 個の観測値から作る

$$F_n(x)=\frac1n\sum_{k=1}^{n}\mathbf 1_{\{X_k\le x\}}=\frac{\#\{k\le n:X_k\le x\}}{n}\qquad(x\in\mathbb{R})\tag{12.14}$$

を経験分布関数と呼ぶ。これは観測値の各点で高さ $1/n$ ずつ跳ぶ階段関数である。それ自体が1つの分布関数（観測値に等しい重みを置いた離散分布の分布関数）でもある。

点 $x$ を固定すると、$\mathbf 1_{\{X_k\le x\}}$ は独立同分布で、期待値がその点の分布関数の値に等しい有界な確率変数である。よって強法則から、$F_n(x)\to F(x)$ a.s. である。同じく $\mathbf 1_{\{X_k<x\}}$ に強法則を使うと、左極限についても次が成り立つ。

$$F_n(x-)=\frac1n\#\{k\le n:X_k<x\}\to F(x-)=P(X_1<x)\quad\text{a.s.}$$

しかし確率 1 の事象は点ごとに違い、点は非可算個ある。よってこのままでは「すべての点で同時に」収束すると結論することはできない。分布関数の単調性を使うと、有限個の点での収束から一様収束が得られる。

証明

自然数 $m$ を固定する。第3章 §3.5 の一般化逆関数 $F^{-1}(u)=\inf\{x:F(x)\ge u\}$ を使い、分点を次のように置く。

$$t_j=F^{-1}(j/m)\qquad(j=1,\dots,m-1)$$

さらに $t_0=-\infty$、$t_m=+\infty$ とする。端点では次のように約束する。

$$F(-\infty)=F_n(-\infty)=0,\qquad F(+\infty-)=F_n(+\infty-)=1$$

$0<j/m<1$ だから、分点 $t_1\le\dots\le t_{m-1}$ は有限の実数である。第3章の補題 $F^{-1}(u)\le x\iff u\le F(x)$ から、次の2つが成り立つ。

• $x=t_j$ とすると、$j/m\le F(t_j)$ である。
• $x<t_j$ なら、$F^{-1}(j/m)>x$ である。よって $F(x)<j/m$ である。左から近づけて、$F(t_j-)\le j/m$ を得る。

2つを合わせると、$j=1,\dots,m$ で次が成り立つ（両端では約束した値で確かめられる）。

$$F(t_j-)-F(t_{j-1})\le\frac jm-\frac{j-1}{m}=\frac1m$$

有限個の点での最大のずれを

$$\Delta_n=\max_{0\le j\le m}\Bigl(\left\vert F_n(t_j)-F(t_j)\right\vert +\left\vert F_n(t_j-)-F(t_j-)\right\vert \Bigr)$$

と置く（両端の項は約束により 0）。各項は強法則によって 0 に近づく。よって $2(m-1)$ 個の確率 1 の事象の共通部分の上で、$\Delta_n\to0$ である。

任意の実数は、ある $j$ について $t_{j-1}\le x<t_j$ を満たす（空の区間は飛ばす）。2つの分布関数は単調非減少だから

$$F_n(x)\le F_n(t_j-)\le F(t_j-)+\Delta_n\le F(t_{j-1})+\frac1m+\Delta_n\le F(x)+\frac1m+\Delta_n,$$

$$F_n(x)\ge F_n(t_{j-1})\ge F(t_{j-1})-\Delta_n\ge F(t_j-)-\frac1m-\Delta_n\ge F(x)-\frac1m-\Delta_n$$

である。したがって $D_n\le\frac1m+\Delta_n$ で、$\limsup_nD_n\le\frac1m$ a.s. となる。$m=1,2,\dots$ に対する可算個の確率 1 の事象の共通部分の上で、上極限は 0 である。

分布関数が連続なら、$D_n$ は観測値を小さい順に並べた $X_{(1)}\le\dots\le X_{(n)}$ を使って、有限回の計算で求まる。経験分布関数は区間 $[X_{(i)},X_{(i+1)})$ の上で一定値 $i/n$ をとる。真の分布関数は単調なので、ずれの上限は跳びの直前か直後で達成される。

$$D_n=\max_{1\le i\le n}\max\Bigl(\frac in-F(X_{(i)}),\ F(X_{(i)})-\frac{i-1}{n}\Bigr)\tag{12.15}$$

確認テスト

オリジナル

1. 12.4-1 点 $x$ を固定すると、$nF_n(x)\sim\mathrm{Bin}(n,F(x))$ であることを説明せよ。また $E[F_n(x)]$ と $V[F_n(x)]$ を求めよ。
2. 12.4-2 積分 $I=\int_0^1\mathrm{e}^x\,\mathrm{d} x$ を、次の推定値で推定する（$U_k$ は独立な一様乱数）。

$$\hat I_n=\frac1n\sum_k\mathrm{e}^{U_k}$$

分散 $\sigma^2=V[\mathrm{e}^{U_1}]$ を求めよ。$\sqrt{E[(\hat I_n-I)^2]}\le0.001$ となる最小の $n$ を求めよ。3. 12.4-3 12.4-2 で、$U$ と $1-U$ の組を使う次の推定値を考える（対称変量法）。

$$\check I_n=\frac1n\sum_{k=1}^{n}\frac12\bigl(\mathrm{e}^{U_k}+\mathrm{e}^{1-U_k}\bigr)$$

分散 $V\bigl[\frac12(\mathrm{e}^{U}+\mathrm{e}^{1-U})\bigr]$ を求めよ。この推定値は関数の値を $2n$ 回計算する。同じ回数を使う $\hat I_{2n}$ と分散を比べよ。4. 12.4-4 密度 $f$ と $q$ をとる。$q$ は集合 $\{x:g(x)f(x)\neq0\}$ の上で正とする。また $\int\left\vert g\right\vert f\,\mathrm{d} x<\infty$ とする。密度 $q$ から独立に作った列による重点抽出の推定値を、$\tilde I_n$ とする。これが確率 1 で $\int gf\,\mathrm{d} x$ に収束することを示せ。