# Monte Carlo 積分と分散低減：重点サンプリング・制御変量・層別 > 標本平均による積分の誤差が次元によらず σ/√N であることと、標準誤差・信頼区間の根拠を証明する。積型求積則の誤差 N^(-r/d) と比べたうえで、対称変量・制御変量・層別・条件づけ・重点サンプリングの分散を定理として評価する。 > https://rikai.mugen-giken.com/computer-science/stochastic-computation/monte-carlo-integration ## 0. この記事の要点 たくさんの変数が入った積分は、式のとおりに計算するとあっという間に手に負えなくなります。そこで、乱数を何回も引いて、その結果の平均から積分の値を見積もる、というやり方があります。この記事では、その見積もりがどれだけ当たるか（誤差をどう言い表すか）と、同じ回数の乱数でもっと当てる工夫（分散低減）を、証明つきで学びます。読み終えると、推定値に誤差の幅を添えて報告でき、目的に合わせて工夫を選べるようになります。 - 積分 $I = \mathbb{E}[f(X)]$ を標本平均で推定すると、推定量は不偏で、平均二乗誤差の平方根は $\sigma/\sqrt N$ です。この速さは**次元によりません**。標本分散から標準誤差を推定でき、中心極限定理と Slutsky の定理（確率収束する量で割っても分布収束が保たれるという定理）から漸近信頼区間が得られます（）。 - $d$ 次元で $r$ 階の滑らかさを使う積型求積則の誤差は $N^{-r/d}$ の程度で、$d > 2r$ では Monte Carlo 法の $N^{-1/2}$ が上回ります（）。 - Monte Carlo 法の改良は、指数 $-1/2$ ではなく**定数 $\sigma$ を小さくする**ことです。対称変量は単調な被積分関数で分散を減らし（）、制御変量は相関 $\rho$ のある既知の量を使って分散を $(1 - \rho^2)$ 倍にします（）。 - 層別サンプリングの比例配分は、分散から層間の成分を取り除きます（）。[条件付き期待値](/mathematics/probability/conditional-expectation#def-cond-exp)で置き換えた推定量の分散は元より大きくなりません（）。 - 重点サンプリングの分散は $\int f^2p^2/q - I^2$ で、最適な提案密度は $|f|p$ に比例します。提案の裾が目標より軽いと分散が無限大になります（、）。 - 正規化定数が未知のときは自己正規化推定量を使います。一致性をもち、偏りは $O(1/N)$ です（）。 ## 1. 動機 $d$ 次元の積分を数値的に求めるとき、各軸に $n$ 個の点を置く求積則は $N = n^d$ 回の関数評価を要します。$d = 10$ で各軸にわずか 10 点でも $10^{10}$ 回です。統計力学の分配関数、Bayes 統計の事後期待値、金融派生商品の価格、粒子輸送の計算では、$d$ は数十から数百万に及びます。 Monte Carlo 法は、積分を[期待値](/mathematics/probability/random-variables#def-expectation) $I = \mathbb{E}[f(X)]$ とみて、独立な標本 $X_1, \ldots, X_N$ の平均 $\frac1N\sum f(X_i)$ で推定します。大数の法則が収束を、中心極限定理が誤差の大きさ $\sigma/\sqrt N$ を保証します（[極限定理](/mathematics/probability/limit-theorems)の で $\pi$ の推定を例に見ました）。$N^{-1/2}$ は遅い収束です。精度を 1 桁上げるのに 100 倍の標本が要ります。しかしこの指数は次元に依存しません。高次元では、遅くても次元によらない方法が、速くても次元とともに劣化する方法に勝ちます。 この章の主題は 2 つあります。第一に、**推定値に誤差を付ける**ことです。Monte Carlo 法の出力は確率変数なので、「$0.7861$」という数値だけでは情報がなく、「$0.7861 \pm 0.0004$」と標準誤差を添えてはじめて結果になります。その根拠が §2 の定理です。第二に、**同じ計算量で誤差を小さくする**ことです。指数 $-1/2$ は独立な標本を使うかぎり変えられませんが、定数 $\sigma$ は推定量の作り方で何桁も変わります。§4〜§6 で、代表的な 5 つの方法の分散を定理として評価します。 以下、確率・期待値・分散を $\mathbb{P}$、$\mathbb{E}$、$\mathrm{Var}$ と書き、定義と基本性質は[確率変数と期待値](/mathematics/probability/random-variables)（、）に従います。 ## 2. Monte Carlo 推定量と誤差 $X$ を確率変数（$\mathbb{R}^d$ 値でよい）、$f$ を $f(X)$ が可積分になる可測関数とし、$I = \mathbb{E}[f(X)]$ を求めたい量とする。$X$ と同じ分布に従う独立な確率変数 $X_1, \ldots, X_N$ に対して $$ \hat I_N = \frac1N\sum_{i=1}^N f(X_i) $$ を $I$ の **Monte Carlo 推定量**という。$\sigma^2 = \mathrm{Var}(f(X))$ が有限のとき、$\sigma/\sqrt N$ を $\hat I_N$ の**標準誤差**という。標本分散 $$ s_N^2 = \frac{1}{N-1}\sum_{i=1}^N\bigl(f(X_i) - \hat I_N\bigr)^2 \qquad (N \ge 2) $$ を使った $s_N/\sqrt N$ を**標準誤差の推定値**という。 の設定で次が成り立つ。 1. $\mathbb{E}[\hat I_N] = I$ であり、$N \to \infty$ で $\hat I_N \to I$（概収束）。 2. $\sigma^2 < \infty$ ならば $\mathrm{Var}(\hat I_N) = \sigma^2/N$ で、平均二乗誤差の平方根は $\sqrt{\mathbb{E}[(\hat I_N - I)^2]} = \sigma/\sqrt N$ である。 3. $\sigma^2 < \infty$ ならば $\mathbb{E}[s_N^2] = \sigma^2$ であり、$s_N^2 \to \sigma^2$（概収束）。 4. $0 < \sigma^2 < \infty$ ならば、$T_N = \sqrt N(\hat I_N - I)/s_N$（$s_N = 0$ のときは $T_N = 0$ と定める）は標準正規分布に分布収束する。したがって任意の $z > 0$ について $$ \lim_{N\to\infty}\mathbb{P}\Bigl(\hat I_N - z\frac{s_N}{\sqrt N} \le I \le \hat I_N + z\frac{s_N}{\sqrt N}\Bigr) = 2\Phi(z) - 1 $$ である（$\Phi$ は標準正規分布の分布関数）。 $Y_i = f(X_i)$ と置きます。$Y_1, Y_2, \ldots$ は独立同分布で、$\mathbb{E}[Y_i] = I$ です。 **1.** 期待値の線形性から $\mathbb{E}[\hat I_N] = \frac1N\sum_i\mathbb{E}[Y_i] = I$ です。概収束は大数の強法則（）です。 **2.** 独立な $Y_i, Y_j$（$i \ne j$）の共分散は $0$ なので（）、分散の和の公式（ の 2 と 4）から $\mathrm{Var}(\hat I_N) = \frac{1}{N^2}\sum_i\mathrm{Var}(Y_i) = \sigma^2/N$ です。$\mathbb{E}[\hat I_N] = I$ なので $\mathbb{E}[(\hat I_N - I)^2] = \mathrm{Var}(\hat I_N)$ です。 **3.** $\sum_i(Y_i - \hat I_N)^2 = \sum_iY_i^2 - N\hat I_N^2$ なので $$ s_N^2 = \frac{N}{N-1}\Bigl(\frac1N\sum_{i=1}^NY_i^2 - \hat I_N^2\Bigr) $$ です。期待値を取ると、$\mathbb{E}[Y_i^2] = \sigma^2 + I^2$ と $\mathbb{E}[\hat I_N^2] = \mathrm{Var}(\hat I_N) + I^2 = \sigma^2/N + I^2$ から $$ \mathbb{E}[s_N^2] = \frac{N}{N-1}\Bigl(\sigma^2 + I^2 - \frac{\sigma^2}{N} - I^2\Bigr) = \sigma^2 $$ です。概収束については、$Y_i^2$ は可積分な独立同分布の列なので、大数の強法則から $\frac1N\sum Y_i^2 \to \sigma^2 + I^2$、また 1 から $\hat I_N^2 \to I^2$、$N/(N-1) \to 1$ です。よって $s_N^2 \to \sigma^2$ です。 **4.** 中心極限定理（）から $Z_N = \sqrt N(\hat I_N - I)/\sigma$ は標準正規分布に分布収束します。$R_N$ を、$s_N > 0$ のとき $s_N/\sigma$、$s_N = 0$ のとき $1$ と定めます。3 から $s_N \to \sigma > 0$（概収束）なので $R_N \to 1$（概収束）で、確率収束もします（ の 1）。Slutsky の定理（）から $Z_N/R_N$ は標準正規分布に分布収束します。$T_N$ と $Z_N/R_N$ が異なるのは $s_N = 0$ のときだけで、その確率は $0$ に収束するので、$T_N - Z_N/R_N$ は $0$ に確率収束し、再び Slutsky の定理（和の形）から $T_N$ も標準正規分布に分布収束します。 区間の式を示します。$s_N > 0$ のとき、$I$ が区間に入ることは $|T_N| \le z$ と同値です。$\Phi$ は連続なので、分布収束（）からすべての $x$ で $\mathbb{P}(T_N \le x) \to \Phi(x)$ です。任意の $\varepsilon > 0$ について $\mathbb{P}(T_N \le -z - \varepsilon) \le \mathbb{P}(T_N < -z) \le \mathbb{P}(T_N \le -z)$ で、両端は $\Phi(-z - \varepsilon)$ と $\Phi(-z)$ に収束し、$\varepsilon \to 0$ で一致するので、$\mathbb{P}(T_N < -z) \to \Phi(-z)$ です。よって $$ \mathbb{P}(|T_N| \le z) = \mathbb{P}(T_N \le z) - \mathbb{P}(T_N < -z) \to \Phi(z) - \Phi(-z) = 2\Phi(z) - 1 $$ です。$\mathbb{P}(s_N = 0) \to 0$ なので、区間の被覆確率も同じ極限をもちます。 2 が Monte Carlo 法の基本の式です。誤差 $\sigma/\sqrt N$ に次元 $d$ は現れません（$\sigma$ は $f$ と分布で決まり、次元とともに増えることも減ることもあります）。4 は、未知の $\sigma$ を標本から推定した $s_N$ で置き換えても、区間 $\hat I_N \pm z\,s_N/\sqrt N$ の被覆確率が漸近的に正しいことを保証します。$z = 1$ で約 $68\%$、$z = 1.96$ で $95\%$、$z = 3$ で約 $99.7\%$ です。この区間を、信頼係数 $2\Phi(z) - 1$ の**漸近信頼区間**といいます（信頼集合の定義は 。「漸近」は、被覆確率が極限でだけ保証されることを指します）。この章と以降の章のコードでは、推定値に必ず $s_N/\sqrt N$ を添え、理論値との差を標準誤差の倍数で示します。 仮定 $\sigma^2 < \infty$ は本質的です。$f(X)$ が可積分でも分散が無限大なら、$\hat I_N$ は収束しますが $s_N^2$ は発散し、区間は意味をもちません。重点サンプリングで提案分布の裾が軽いときに実際に起きます（§6）。 有限の $N$ で何が言えるかも確かめておきます。 の設定で $0 < \sigma^2 < \infty$ とする。 1. すべての $N \ge 1$ と $k > 0$ について $\mathbb{P}\bigl(|\hat I_N - I| \ge k\sigma/\sqrt N\bigr) \le 1/k^2$。 2. $\rho_3 = \mathbb{E}[|f(X) - I|^3] < \infty$ ならば、絶対定数 $C$ があって、すべての $N \ge 1$ と $z > 0$ について $$ \Bigl|\mathbb{P}\Bigl(|\hat I_N - I| \le \frac{z\sigma}{\sqrt N}\Bigr) - \bigl(2\Phi(z) - 1\bigr)\Bigr| \le \frac{2C\rho_3}{\sigma^3\sqrt N}. $$ **1.** の 1・2 から $\hat I_N$ の平均は $I$、分散は $\sigma^2/N$ です。Chebyshev の不等式（）を $a = k\sigma/\sqrt N$ で使うと、確率は $(\sigma^2/N)/a^2 = 1/k^2$ 以下です。 **2.** $Z_N = \sqrt N(\hat I_N - I)/\sigma$ の分布関数を $F_N$ とし、$\delta = C\rho_3/(\sigma^3\sqrt N)$ と置きます。Berry–Esseen の定理（）から、すべての $x$ で $|F_N(x) - \Phi(x)| \le \delta$ です。$\mathbb{P}(|Z_N| \le z) = F_N(z) - \mathbb{P}(Z_N < -z)$ で、$\mathbb{P}(Z_N < -z) = \lim_{x \uparrow -z}F_N(x)$ です。各 $x$ で $|F_N(x) - \Phi(x)| \le \delta$ なので、$x \uparrow -z$ の極限を取ると、$\Phi$ の連続性から $|\mathbb{P}(Z_N < -z) - \Phi(-z)| \le \delta$ です。三角不等式から $$ \bigl|\mathbb{P}(|Z_N| \le z) - (\Phi(z) - \Phi(-z))\bigr| \le 2\delta $$ です。 1 は分布の形によらない保証で、$95\%$ の区間には $k = \sqrt{20} \approx 4.47$ が要ります。正規近似の $1.96$ の 2.3 倍の幅です。2 は、正規近似による被覆確率の誤差が $N^{-1/2}$ の速さで消えること、ただしその定数が歪み $\rho_3/\sigma^3$ に比例することを示します。稀な事象の指示関数（$\mathbb{P}(f = 1) = p$ が小さい）では $\rho_3/\sigma^3 \approx 1/\sqrt p$ なので、$Np$ が大きくなるまで正規近似は信用できません。 1 標本あたりの分散が $\sigma^2$、計算時間が $t$ の推定量で、全体の計算時間 $T$ を使うと、標本数は $N = T/t$、分散は $\sigma^2t/T$ である。$1/(\sigma^2t)$ をこの推定量の**効率**という。2 つの推定量の比較は、分散ではなく $\sigma^2t$ で行う。 分散を半分にする工夫でも、1 標本の計算時間が 3 倍になるなら損です。以下の定理は分散だけを述べますが、実際の選択ではこの積を見ます。 ## 3. 求積則との比較 Monte Carlo 法が有利になるのはどこからかを、定理で確かめます。求積則とは、積分を「関数の値の重みつき和」で近似する公式のことです（台形則や Simpson 則が代表例です）。1 次元の公式を各軸に使って $d$ 次元に広げたものが積型求積則で、その誤差を見積もるのが次の定理の目的です。証明は、変数を 1 つずつ「積分」から「求積則」に置き換えていき、1 回の置き換えの誤差を 1 次元の誤差評価で抑える、という方針です。$[0,1]$ 上の求積則 $Q g = \sum_{j} w_jg(x_j)$ から、$[0,1]^d$ 上の**積型求積則** $$ Q^{\otimes d}f = \sum_{j_1, \ldots, j_d} w_{j_1}\cdots w_{j_d}\,f(x_{j_1}, \ldots, x_{j_d}) $$ を作ります。節点が $J$ 個なら、評価回数は $N = J^d$ です。 $r \ge 1$、$n \ge 1$ とする。$[0,1]$ 上の求積則 $Q_n$ が、重み $w_j \ge 0$、$\sum_jw_j = 1$ をもち、ある定数 $K$ について、すべての $g \in C^r[0,1]$ で $$ \Bigl|\int_0^1 g(x)\,dx - Q_ng\Bigr| \le K\,n^{-r}\max_{[0,1]}|g^{(r)}| $$ を満たすとする。 1. $f : [0,1]^d \to \mathbb{R}$ が連続で、各 $k$ について $\partial^rf/\partial x_k^r$ が存在して連続ならば $$ \Bigl|\int_{[0,1]^d} f(x)\,dx - Q_n^{\otimes d}f\Bigr| \le K\,n^{-r}\sum_{k=1}^d\max_{[0,1]^d}\Bigl|\frac{\partial^rf}{\partial x_k^r}\Bigr|. $$ 2. 複合台形則（$n$ 個の小区間、$r = 2$、$K = \frac1{12}$）と複合 Simpson 則（$n$ は偶数、$r = 4$、$K = \frac1{180}$）はこの仮定を満たし、評価回数 $N = (n+1)^d$ について誤差は $2^rK N^{-r/d}\sum_k\max|\partial_k^rf|$ 以下である。 3. 複合台形則の積型則を $f(x) = x_1^2$ に使ったときの誤差は、ちょうど $\frac{1}{6n^2}$ で、$\frac16N^{-2/d}$ より大きい。 **1.** $k = 0, 1, \ldots, d$ について、「最初の $k$ 個の変数には $Q_n$ を、残りの $d - k$ 個の変数には積分を」使った量 $$ A_k = \sum_{j_1, \ldots, j_k} w_{j_1}\cdots w_{j_k}\int_{[0,1]^{d-k}} f(x_{j_1}, \ldots, x_{j_k}, t_{k+1}, \ldots, t_d)\,dt_{k+1}\cdots dt_d $$ を考えます。$A_0 = \int f$、$A_d = Q_n^{\otimes d}f$ なので、$\int f - Q_n^{\otimes d}f = \sum_{k=1}^d(A_{k-1} - A_k)$ です。$f$ は連続なので、Fubini の定理（）から $A_{k-1}$ の中の積分を「$t_k$ について先に積分する」形に書けて、 $$ A_{k-1} - A_k = \sum_{j_1, \ldots, j_{k-1}} w_{j_1}\cdots w_{j_{k-1}}\int_{[0,1]^{d-k}}\Bigl(\int_0^1 g(t_k)\,dt_k - Q_ng\Bigr)dt_{k+1}\cdots dt_d $$ です。ここで $g(t) = f(x_{j_1}, \ldots, x_{j_{k-1}}, t, t_{k+1}, \ldots, t_d)$ は、他の変数を固定した $t$ の 1 変数関数です。$g \in C^r[0,1]$ で $|g^{(r)}| \le \max|\partial_k^rf|$ なので、仮定から括弧の中の絶対値は $Kn^{-r}\max|\partial_k^rf|$ 以下です。重みは非負で和が $1$、積分領域の体積は $1$ なので、$|A_{k-1} - A_k| \le Kn^{-r}\max|\partial_k^rf|$ です。$k$ について足して結論を得ます。 **2.** 複合則の誤差（）で $[a,b] = [0,1]$、$h = 1/n$ とすると、台形則の誤差の絶対値は $\frac{h^2}{12}|g''(\xi)| \le \frac1{12}n^{-2}\max|g''|$、Simpson 則は $\frac{h^4}{180}|g^{(4)}(\xi'')| \le \frac1{180}n^{-4}\max|g^{(4)}|$ です。重みは、台形則が $\frac h2, h, \ldots, h, \frac h2$、Simpson 則が $\frac h3, \frac{4h}3, \frac{2h}3, \ldots, \frac{4h}3, \frac h3$ で、どちらも非負で和が $1$ です（定数関数 $1$ を正確に積分するからです）。節点は $n + 1$ 個なので $N = (n+1)^d \le (2n)^d$、すなわち $n \ge \frac12N^{1/d}$ で、$n^{-r} \le 2^rN^{-r/d}$ です。 **3.** $f$ は $x_1$ だけの関数で、他の変数についての和は重みの和 $1$ を掛けるだけなので、$Q_n^{\otimes d}f$ は 1 次元の複合台形則 $T_h$ を $g(x) = x^2$ に使った値です。 の 1 で $g'' = 2$ が定数なので、$\int_0^1g - T_h(g) = -\frac{h^2}{12}\cdot2 = -\frac{1}{6n^2}$ です。$N = (n+1)^d > n^d$ から $n^{-2} > N^{-2/d}$ です。 1 と 2 は、$r$ 階の滑らかさを使う積型則の誤差が $N^{-r/d}$ の程度**以下**であること、3 は台形則でこの次数が実際に達成されてしまうことを示します。Monte Carlo 法の $\sigma N^{-1/2}$ と指数を比べると、 $$ \frac rd < \frac12 \iff d > 2r $$ のとき Monte Carlo 法の収束が速くなります。台形則（$r = 2$）なら 5 次元以上、Simpson 則（$r = 4$）なら 9 次元以上です。逆に低次元の滑らかな積分では、求積則が圧倒的に有利です。 の 3 は積型の台形則についての下界です。より一般に、$r$ 階の導関数が有界な関数の全体に対しては、$N$ 個の点を使うどんな決定的な求積則でも、最悪の誤差は $N^{-r/d}$ の定数倍より小さくできないことが知られています（Bakhvalov）。証明は、$N$ 個の点をすべて避ける小さなこぶ状の関数を作る議論で、Novak の本にあります。関数の族を狭める（混合導関数が有界、など）と事情が変わり、それが準 Monte Carlo 法の主題です。 次のコードは、被積分関数 $f(x) = \prod_{k=1}^d\frac\pi2\sin(\pi x_k)$（積分の値は $1$）について、積型 Simpson 則と Monte Carlo 法を同じ評価回数で比べます。 ```python # Python 3.9.6, NumPy 2.0.2 import numpy as np f1 = lambda x: 0.5 * np.pi * np.sin(np.pi * x) # 1 次元の因子。積分は 1 def simpson_1d(n): # n 個の小区間（n は偶数）の複合 Simpson 則 x = np.linspace(0, 1, n + 1) w = np.ones(n + 1); w[1:-1:2] = 4; w[2:-1:2] = 2 return (w * f1(x)).sum() / (3 * n) rng = np.random.default_rng(20261008) print(" d n N=(n+1)^d 積型 Simpson の誤差 Monte Carlo の誤差 (標準誤差)") for d, n in ((2, 30), (6, 4), (10, 2)): N = (n + 1) ** d q = simpson_1d(n) ** d # 被積分関数が積なので、積型則は 1 次元の値の d 乗 y = np.prod(f1(rng.random((N, d))), axis=1) est, se = y.mean(), y.std(ddof=1) / np.sqrt(N) print(f"{d:2d} {n:2d} {N:10d} {abs(q - 1):.3e} {abs(est - 1):.3e} ({se:.3e})") ``` ```text d n N=(n+1)^d 積型 Simpson の誤差 Monte Carlo の誤差 (標準誤差) 2 30 961 1.338e-06 6.344e-03 (2.308e-02) 6 4 15625 1.376e-02 8.516e-03 (1.260e-02) 10 2 59049 5.859e-01 1.618e-02 (1.124e-02) ``` 2 次元では Simpson 則が 4 桁よく、6 次元でほぼ並び、10 次元では Simpson 則の誤差が $59\%$ に達するのに対して Monte Carlo 法は約 $1.6\%$ です。Monte Carlo 法の実際の誤差は、3 つの場合とも標準誤差の 1.5 倍以内にあります。定理が予言するとおり、交代は $d = 2r = 8$ の前後で起きています。 ## 4. 対称変量と制御変量 ここからは、同じ $I$ を、より分散の小さい確率変数の期待値として書き直す方法を扱います。 ### 4.1. 対称変量 $U \sim \mathrm{U}(0,1)$ なら $1 - U \sim \mathrm{U}(0,1)$ です。$f(U)$ と $f(1 - U)$ は同じ分布に従い、$U$ が大きいとき $1 - U$ は小さいので、$f$ が単調なら 2 つの値は逆向きに動きます。その平均を取ると揺らぎが打ち消し合います。 $U_1, \ldots, U_{N/2}$（$N$ は偶数）を $[0,1]^d$ 上の一様分布に従う独立な確率変数とする。$\boldsymbol{1} = (1, \ldots, 1)$ として $$ \hat I_N^{\mathrm{anti}} = \frac1N\sum_{i=1}^{N/2}\bigl(f(U_i) + f(\boldsymbol{1} - U_i)\bigr) $$ を**対称変量**（antithetic variates）による推定量という。関数の評価回数は $N$ である。 効果の根拠は、単調な関数どうしの[共分散](/mathematics/probability/random-variables#def-variance-covariance)の符号です。次の補題は、同じ向きに増える 2 つの量は一緒に動く（共分散が非負）という当たり前の感覚を、多次元でも確かめるものです。証明は、1 次元で示したあと、変数を 1 つずつ積分して次元を下げる帰納法で進みます。 $g, h : [0,1]^d \to \mathbb{R}$ をボレル可測で、各変数について単調非減少（他の変数を固定したとき、その変数の関数として非減少）とし、$U$ を $[0,1]^d$ 上の一様分布に従う確率変数、$g(U), h(U) \in L^2$ とする。このとき $$ \mathrm{Cov}\bigl(g(U), h(U)\bigr) \ge 0. $$ $d$ についての帰納法で示します。 **$d = 1$。** $U'$ を $U$ と独立で同じ分布に従う確率変数とします。$g, h$ が非減少なので、$g(U) - g(U')$ と $h(U) - h(U')$ は（$U \ge U'$ ならともに非負、$U \le U'$ ならともに非正で）積が常に非負です。期待値を取って展開すると、独立性（）と同分布性から $$ 0 \le \mathbb{E}\bigl[(g(U) - g(U'))(h(U) - h(U'))\bigr] = 2\,\mathbb{E}[g(U)h(U)] - 2\,\mathbb{E}[g(U)]\,\mathbb{E}[h(U)] = 2\,\mathrm{Cov}(g(U), h(U)) $$ です。 **$d \ge 2$。** $u = (u', t)$（$u' \in [0,1]^{d-1}$、$t \in [0,1]$）と書きます。Fubini の定理（）を $g^2, h^2$ に使うと、ほとんどすべての $t$ について $g(\cdot, t), h(\cdot, t) \in L^2([0,1]^{d-1})$ です。そのような $t$ について $$ G(t) = \int g(u', t)\,du', \qquad H(t) = \int h(u', t)\,du' $$ と置きます。$g(\cdot, t)$、$h(\cdot, t)$ は $d - 1$ 変数の各変数について非減少なので、帰納法の仮定から $$ \int g(u', t)h(u', t)\,du' \ge G(t)H(t) $$ です。$G, H$ は $t$ について非減少で（$g, h$ が $t$ について非減少だからです）、$G(t)^2 \le \int g(u', t)^2du'$ から $L^2(0,1)$ に属します。$d = 1$ の場合を $G, H$ に使うと $\int_0^1G(t)H(t)\,dt \ge \int_0^1G\,\int_0^1H$ です。Fubini の定理でつなぐと $$ \mathbb{E}[g(U)h(U)] = \int_0^1\Bigl(\int g(u',t)h(u',t)\,du'\Bigr)dt \ge \int_0^1G(t)H(t)\,dt \ge \int_0^1G\,\int_0^1H = \mathbb{E}[g(U)]\,\mathbb{E}[h(U)] $$ です。 この補題は、有限和についての Chebyshev の和の不等式（同じ向きに並んだ 2 つの数列の積の平均は、平均の積以上）の積分版です。 $f : [0,1]^d \to \mathbb{R}$ をボレル可測で、各変数について単調（変数ごとに、非減少か非増加かのどちらかに決まっている）とし、$U$ を $[0,1]^d$ 上の一様分布に従う確率変数、$f(U) \in L^2$、$\sigma^2 = \mathrm{Var}(f(U))$ とする。 1. $\mathrm{Cov}\bigl(f(U), f(\boldsymbol{1} - U)\bigr) \le 0$。 2. $\hat I_N^{\mathrm{anti}}$ は不偏で、 $$ \mathrm{Var}\bigl(\hat I_N^{\mathrm{anti}}\bigr) = \frac{\sigma^2 + \mathrm{Cov}(f(U), f(\boldsymbol{1} - U))}{N} \le \frac{\sigma^2}{N} $$ である。すなわち、同じ評価回数 $N$ の単純な推定量 $\hat I_N$ の分散以下である。 **1.** まず $f$ がすべての変数について非減少の場合を考えます。$g(u) = f(u)$、$h(u) = -f(\boldsymbol{1} - u)$ と置くと、どちらも各変数について非減少です（$u_k$ が増えると $1 - u_k$ が減り、$f$ が減るので $-f$ が増えます）。$\boldsymbol{1} - U$ は $U$ と同じ分布に従うので $h(U) \in L^2$ です。 から $\mathrm{Cov}(f(U), -f(\boldsymbol{1} - U)) \ge 0$、すなわち主張が従います。 一般の場合は、$f$ が非増加である変数の番号の集合を $D$ とし、写像 $\tau : [0,1]^d \to [0,1]^d$ を「$k \in D$ の成分を $1 - u_k$ に替え、他はそのまま」と定めます。$\tilde f = f \circ \tau$ はすべての変数について非減少です。$V = \tau(U)$ は $[0,1]^d$ 上の一様分布に従い、$\tau \circ \tau$ は恒等写像で、$\tau(\boldsymbol{1} - u) = \boldsymbol{1} - \tau(u)$ なので、 $$ f(U) = \tilde f(V), \qquad f(\boldsymbol{1} - U) = \tilde f\bigl(\tau(\boldsymbol{1} - U)\bigr) = \tilde f(\boldsymbol{1} - V) $$ です。前半を $\tilde f$ と $V$ に使えば $\mathrm{Cov}(f(U), f(\boldsymbol{1} - U)) = \mathrm{Cov}(\tilde f(V), \tilde f(\boldsymbol{1} - V)) \le 0$ です。 **2.** $W_i = \frac12(f(U_i) + f(\boldsymbol{1} - U_i))$ と置くと、$\hat I_N^{\mathrm{anti}} = \frac{1}{N/2}\sum_{i=1}^{N/2}W_i$ で、$W_i$ は独立同分布です。$\mathbb{E}[W_i] = \frac12(I + I) = I$ なので不偏です。$\mathrm{Var}(f(\boldsymbol{1} - U)) = \sigma^2$ なので、分散の和の公式（ の 2 と 4）から $$ \mathrm{Var}(W_i) = \frac14\bigl(\sigma^2 + \sigma^2 + 2\,\mathrm{Cov}(f(U), f(\boldsymbol{1} - U))\bigr) = \frac{\sigma^2 + \mathrm{Cov}(f(U), f(\boldsymbol{1} - U))}{2} $$ です。 の 2 を $W_i$ の $N/2$ 個の平均に使うと $\mathrm{Var}(\hat I_N^{\mathrm{anti}}) = \mathrm{Var}(W_1)/(N/2)$ で、主張の式になります。不等号は 1 です。 単調性の仮定を外すと、対称変量は逆効果になりえます。$f(u) = (u - \frac12)^2$ では $f(1 - u) = f(u)$ なので共分散は $\sigma^2$ で、分散は $2\sigma^2/N$ と 2 倍になります（同じ点を 2 回評価しているのと同じです）。逆関数法で作った乱数（$X = F^{-}(U)$ は $U$ について非減少）を単調な関数に入れる計算では、仮定が自然に満たされます。標準誤差は $W_i$ の標本分散から推定します。$f(U_i)$ と $f(\boldsymbol{1} - U_i)$ を独立な $N$ 個の標本として扱ってはいけません。 ### 4.2. 制御変量 $Y = f(X)$ と同じ標本から計算でき、期待値 $\nu = \mathbb{E}[C]$ が**既知**の確率変数 $C = c(X)$ を**制御変量**という。係数 $\beta \in \mathbb{R}$ に対して $$ \hat I_N^{\mathrm{cv}}(\beta) = \frac1N\sum_{i=1}^N\bigl(Y_i - \beta(C_i - \nu)\bigr) $$ を制御変量による推定量という。 $Y, C \in L^2$、$\sigma^2 = \mathrm{Var}(Y)$、$\sigma_C^2 = \mathrm{Var}(C) > 0$、$\rho = \mathrm{Cov}(Y, C)/(\sigma\sigma_C)$（$\sigma = 0$ のときは $\rho = 0$）とする。 1. 任意の定数 $\beta$ について $\hat I_N^{\mathrm{cv}}(\beta)$ は不偏で、$\mathrm{Var}\bigl(\hat I_N^{\mathrm{cv}}(\beta)\bigr) = \frac1N\bigl(\sigma^2 - 2\beta\,\mathrm{Cov}(Y, C) + \beta^2\sigma_C^2\bigr)$。 2. この分散は $\beta^* = \mathrm{Cov}(Y, C)/\sigma_C^2$ で最小になり、最小値は $(1 - \rho^2)\sigma^2/N$ である。 3. $\beta$ を標本から決めた確率変数 $\hat\beta \in L^2$ で置き換えると、偏りは $\mathbb{E}[\hat I_N^{\mathrm{cv}}(\hat\beta)] - I = -\mathrm{Cov}(\hat\beta, \bar C_N)$（$\bar C_N = \frac1N\sum C_i$）で、 $$ \bigl|\mathbb{E}[\hat I_N^{\mathrm{cv}}(\hat\beta)] - I\bigr| \le \frac{\sigma_C}{\sqrt N}\sqrt{\mathrm{Var}(\hat\beta)} $$ を満たす。$\hat\beta$ が $(C_1, \ldots, C_N)$ と独立（別の標本から推定した場合）なら、偏りは $0$ である。 **1.** $W_i = Y_i - \beta(C_i - \nu)$ は独立同分布で、$\mathbb{E}[W_i] = I - \beta(\nu - \nu) = I$ です。$\mathrm{Var}(W_i) = \mathrm{Var}(Y - \beta C) = \sigma^2 - 2\beta\,\mathrm{Cov}(Y, C) + \beta^2\sigma_C^2$（ の 2〜4）で、 の 2 から $N$ で割ったものが推定量の分散です。 **2.** $\beta$ の 2 次式を平方完成すると $$ \sigma^2 - 2\beta\,\mathrm{Cov}(Y, C) + \beta^2\sigma_C^2 = \sigma_C^2\Bigl(\beta - \frac{\mathrm{Cov}(Y, C)}{\sigma_C^2}\Bigr)^2 + \sigma^2 - \frac{\mathrm{Cov}(Y, C)^2}{\sigma_C^2} $$ なので、$\beta = \beta^*$ で最小値 $\sigma^2 - \mathrm{Cov}(Y,C)^2/\sigma_C^2 = (1 - \rho^2)\sigma^2$ を取ります。Cauchy–Schwarz の不等式（）から $\rho^2 \le 1$ で、これは非負です。 **3.** $\hat I_N^{\mathrm{cv}}(\hat\beta) = \bar Y_N - \hat\beta(\bar C_N - \nu)$ で、$\mathbb{E}[\bar Y_N] = I$、$\mathbb{E}[\bar C_N] = \nu$ なので、偏りは $-\mathbb{E}[\hat\beta(\bar C_N - \nu)] = -\mathrm{Cov}(\hat\beta, \bar C_N)$ です。Cauchy–Schwarz の不等式と $\mathrm{Var}(\bar C_N) = \sigma_C^2/N$ から評価が従います。$\hat\beta$ が $\bar C_N$ と独立なら共分散は $0$ です（）。 $\rho = \pm0.9$ で分散は $0.19$ 倍、$\rho = \pm0.99$ で約 $0.02$ 倍です。制御変量は、被積分関数の「解析的に積分できる近似」を引いて、残りだけを Monte Carlo 法で積分する方法と読めます。実際には $\mathrm{Cov}(Y, C)$ は未知なので、$\beta^*$ を標本共分散と標本分散の比 $\hat\beta$ で推定します。3 から、$\hat\beta$ の標準偏差が $O(N^{-1/2})$ であれば（4 次モーメントが有限な通常の状況です）偏りは $O(1/N)$ で、標準誤差 $O(N^{-1/2})$ に比べて無視できます。偏りを完全に避けたいときは、小さな予備の標本で $\hat\beta$ を決め、本番の標本ではそれを定数として使います。 $I = \int_0^1e^x\,dx = e - 1$ を $Y = e^U$（$U \sim \mathrm{U}(0,1)$）の期待値として推定します。必要な量はすべて閉じた形で求まります。 - **単純な推定量。** $\sigma^2 = \mathbb{E}[e^{2U}] - I^2 = \frac{e^2 - 1}{2} - (e - 1)^2 \approx 0.24204$。 - **対称変量。** $\mathbb{E}[e^Ue^{1-U}] = e$ なので $\mathrm{Cov}(e^U, e^{1-U}) = e - (e-1)^2 \approx -0.23421$。 から分散は $(\sigma^2 - 0.23421)/N \approx 0.00782/N$ で、単純な推定量の $0.0323$ 倍（約 $1/31$）。 - **制御変量 $C = U$。** $\nu = \frac12$、$\sigma_C^2 = \frac1{12}$、$\mathbb{E}[Ue^U] = \int_0^1xe^x\,dx = 1$ から $\mathrm{Cov}(e^U, U) = 1 - \frac{e-1}{2} \approx 0.14086$。$\beta^* = 12 \times 0.14086 \approx 1.6903$、$\rho^2 = 0.14086^2/(0.24204/12) \approx 0.98372$。 から分散は $0.0163$ 倍（約 $1/61$）。 - **層別（10 層、比例配分）。** §5 の と から $0.0110$ 倍（約 $1/91$）。 どの方法も、関数の評価回数は同じ $N$ です。 の「当たり外れ」の推定量 $\mathbf{1}\{U_1^2 + U_2^2 \le 1\}$ の分散は $\frac\pi4(1 - \frac\pi4) \approx 0.1685$ でした。$U_2$ について先に積分すると $\mathbb{P}(U_2 \le \sqrt{1 - U_1^2} \mid U_1) = \sqrt{1 - U_1^2}$ なので、$Y = \sqrt{1 - U^2}$ も期待値 $\pi/4$ をもちます。その分散は $\mathbb{E}[1 - U^2] - \frac{\pi^2}{16} = \frac23 - \frac{\pi^2}{16} \approx 0.0498$ で、$0.296$ 倍です（§5 の の例になっています）。 この $Y$ にさらに制御変量 $C = U^2$ を使います。$\nu = \frac13$、$\sigma_C^2 = \frac15 - \frac19 = \frac4{45}$、$\mathbb{E}[U^2\sqrt{1 - U^2}] = \int_0^{\pi/2}\sin^2\theta\cos^2\theta\,d\theta = \frac\pi{16}$ から $$ \mathrm{Cov}(Y, C) = \frac\pi{16} - \frac\pi4\cdot\frac13 = -\frac{\pi}{48}, \qquad \rho^2 = \frac{(\pi/48)^2}{\frac4{45}\bigl(\frac23 - \frac{\pi^2}{16}\bigr)} \approx 0.9674 $$ です。分散はさらに $1 - \rho^2 \approx 0.0326$ 倍になり、最初の推定量と比べると $0.296 \times 0.0326 \approx 0.0096$ 倍、約 $1/104$ です。$Y$ は $U$ について単調減少なので対称変量も使えて、その場合の分散の比は数値積分で $0.275$ 倍です（$\mathrm{Cov}(Y(U), Y(1-U)) \approx -0.0361$）。この例では、被積分関数の形（$1 - u^2/2$ に近い）に合った制御変量のほうがはるかに効きます。 ## 5. 条件づけと層別 ### 5.1. 条件付き期待値による置き換え 以下の $\mathcal{G}$ は、「すでに値が分かっているものとする情報」を表します。たとえば $\mathcal{G}$ が $U_1$ の値の情報なら、$U_1$ は固定して残りの乱数だけを平均する、ということです。具体例は で見られます。 $Y \in L^2$、$\mathcal{G}$ を部分 $\sigma$-加法族とすると、 $$ \mathrm{Var}(Y) = \mathbb{E}\bigl[\mathrm{Var}(Y \mid \mathcal{G})\bigr] + \mathrm{Var}\bigl(\mathbb{E}[Y \mid \mathcal{G}]\bigr). $$ $\hat Y = \mathbb{E}[Y \mid \mathcal{G}]$ と置きます（条件付き期待値の定義は ）。$\hat Y \in L^2$ です（ の 1）。条件付き分散の定義（）とタワー・プロパティ（）から $$ \mathbb{E}\bigl[\mathrm{Var}(Y \mid \mathcal{G})\bigr] = \mathbb{E}\bigl[\mathbb{E}[Y^2 \mid \mathcal{G}]\bigr] - \mathbb{E}[\hat Y^2] = \mathbb{E}[Y^2] - \mathbb{E}[\hat Y^2] $$ です。また $\mathbb{E}[\hat Y] = \mathbb{E}[Y]$ から $\mathrm{Var}(\hat Y) = \mathbb{E}[\hat Y^2] - (\mathbb{E}[Y])^2$ です。2 つを足すと $\mathbb{E}[Y^2] - (\mathbb{E}[Y])^2 = \mathrm{Var}(Y)$ です。 $Y \in L^2$、$I = \mathbb{E}[Y]$、$\mathcal{G}$ を部分 $\sigma$-加法族とし、$\hat Y = \mathbb{E}[Y \mid \mathcal{G}]$ と置く。このとき $\mathbb{E}[\hat Y] = I$ で、 $$ \mathrm{Var}(\hat Y) = \mathrm{Var}(Y) - \mathbb{E}\bigl[\mathrm{Var}(Y \mid \mathcal{G})\bigr] \le \mathrm{Var}(Y) $$ である。したがって、$\hat Y$ の独立な複製の標本平均は $I$ の不偏推定量で、その分散は $Y$ の標本平均の分散以下である。 $\mathbb{E}[\hat Y] = I$ はタワー・プロパティ（）です。分散の式は で、条件付き分散はほとんど確実に非負なので（条件付き Jensen の不等式。）、その期待値は非負です。標本平均についての主張は の 1・2 です。 使い方は、「乱数の一部について解析的に積分できるなら、その分は乱数を使わずに積分してしまう」です。 で $\mathbf{1}\{U_1^2 + U_2^2 \le 1\}$ を $\sqrt{1 - U_1^2}$ に替えたのは、$\mathcal{G} = \sigma(U_1)$ についての条件付き期待値を取ったことにあたります。統計学の Rao–Blackwell の定理と同じ計算なので、この操作を Rao–Blackwell 化とも呼びます。効率の観点では、$\hat Y$ の計算が $Y$ より重くなる分を で差し引いて評価します。 ### 5.2. 層別サンプリング 層別サンプリングは、値域をいくつかの区画（層）に分け、各層から決まった個数ずつ標本を取る方法です。どの層に何個の標本を割り当てるか（配分）が分散を左右します。 $X$ の値域を互いに交わらない可測集合 $A_1, \ldots, A_K$（**層**）に分け、$p_k = \mathbb{P}(X \in A_k) > 0$、$\sum_kp_k = 1$ とする。層 $k$ の中での条件付き分布 $\mathbb{P}(X \in B \mid X \in A_k) = \mathbb{P}(X \in B \cap A_k)/p_k$ に従う確率変数を $X^{(k)}$ とし、 $$ \mu_k = \mathbb{E}[f(X^{(k)})], \qquad \sigma_k^2 = \mathrm{Var}(f(X^{(k)})) $$ と置く。層 $k$ から $n_k \ge 1$ 個の独立な標本 $X^{(k)}_1, \ldots, X^{(k)}_{n_k}$ を取り（異なる層の標本も独立）、 $$ \hat I^{\mathrm{str}} = \sum_{k=1}^K p_k\cdot\frac{1}{n_k}\sum_{i=1}^{n_k} f(X^{(k)}_i) $$ を**層別推定量**という。$n_1 + \cdots + n_K = N$ を満たす $(n_k)$ を**配分**といい、$n_k = Np_k$ を**比例配分**という。 の設定で $f(X) \in L^2$、$I = \mathbb{E}[f(X)]$、$\sigma^2 = \mathrm{Var}(f(X))$ とする。 1. $I = \sum_kp_k\mu_k$ であり、$\sigma^2 = \sum_kp_k\sigma_k^2 + \sum_kp_k(\mu_k - I)^2$（層内分散と層間分散への分解）。 2. $\hat I^{\mathrm{str}}$ は不偏で、$\mathrm{Var}(\hat I^{\mathrm{str}}) = \sum_k p_k^2\sigma_k^2/n_k$。 3. 比例配分（$Np_k$ がすべて整数のとき）では $$ \mathrm{Var}(\hat I^{\mathrm{str}}) = \frac1N\sum_kp_k\sigma_k^2 = \frac{\sigma^2}{N} - \frac1N\sum_kp_k(\mu_k - I)^2 \le \frac{\sigma^2}{N}. $$ 4. 正の実数 $n_1, \ldots, n_K$ で $\sum_kn_k = N$ を満たすもの全体の上で、$\sum_kp_k^2\sigma_k^2/n_k$ は $\frac1N\bigl(\sum_kp_k\sigma_k\bigr)^2$ 以上である。$\sigma_k$ がすべて正のとき、等号は $n_k = Np_k\sigma_k/\sum_jp_j\sigma_j$（**Neyman 配分**）で成り立つ。また $\bigl(\sum_kp_k\sigma_k\bigr)^2 \le \sum_kp_k\sigma_k^2$ である。 **1.** 可積分な $g(X)$ について、条件付き分布の定義から $\mathbb{E}[g(X)\mathbf{1}_{A_k}(X)] = p_k\,\mathbb{E}[g(X^{(k)})]$ です（$g$ が指示関数のときは定義そのもので、一般の $g$ へは期待値の定義の標準的な手順で延びます。）。層は値域の分割なので $\sum_k\mathbf{1}_{A_k}(X) = 1$ で、 $$ \mathbb{E}[g(X)] = \sum_kp_k\,\mathbb{E}[g(X^{(k)})] $$ です。$g = f$ として $I = \sum_kp_k\mu_k$ です。$g = (f - I)^2$ とすると、$\mathbb{E}[(f(X^{(k)}) - I)^2] = \sigma_k^2 + (\mu_k - I)^2$（平均 $\mu_k$ のまわりの分散と、平均のずれの二乗）から分解の式を得ます。 **2.** 層 $k$ の標本平均 $\bar Y_k$ は平均 $\mu_k$、分散 $\sigma_k^2/n_k$ をもち（ の 1・2）、異なる層の標本平均は独立です。よって $\mathbb{E}[\hat I^{\mathrm{str}}] = \sum_kp_k\mu_k = I$、$\mathrm{Var}(\hat I^{\mathrm{str}}) = \sum_kp_k^2\sigma_k^2/n_k$ です（ の 2・4 と ）。 **3.** 2 に $n_k = Np_k$ を代入すると $\frac1N\sum_kp_k\sigma_k^2$ で、1 の分解から残りが従います。 **4.** 数列についての Cauchy–Schwarz の不等式から $$ \Bigl(\sum_kp_k\sigma_k\Bigr)^2 = \Bigl(\sum_k\frac{p_k\sigma_k}{\sqrt{n_k}}\cdot\sqrt{n_k}\Bigr)^2 \le \Bigl(\sum_k\frac{p_k^2\sigma_k^2}{n_k}\Bigr)\Bigl(\sum_kn_k\Bigr) = N\sum_k\frac{p_k^2\sigma_k^2}{n_k} $$ です。等号は 2 つのベクトルが比例するとき、すなわち $p_k\sigma_k/\sqrt{n_k} = \lambda\sqrt{n_k}$、$n_k = p_k\sigma_k/\lambda$ のときで、$\sum n_k = N$ から $\lambda = \sum_jp_j\sigma_j/N$ です。最後の不等式も Cauchy–Schwarz で、$\bigl(\sum_k\sqrt{p_k}\cdot\sqrt{p_k}\sigma_k\bigr)^2 \le \sum_kp_k\cdot\sum_kp_k\sigma_k^2 = \sum_kp_k\sigma_k^2$ です。 3 は、比例配分の層別サンプリングが**決して損をしない**ことを述べています。取り除かれるのは層間分散 $\sum_kp_k(\mu_k - I)^2$ で、層の平均 $\mu_k$ が層ごとに大きく違うほど、つまり $f$ の変動が層の分け方でよく説明されるほど効きます。4 の Neyman 配分は、ばらつきの大きい層に多くの標本を割り当てます。ただし $\sigma_k$ は未知なので予備の標本で推定する必要があり、推定がまずいと比例配分より悪くなりえます。標準誤差は、層ごとの標本分散 $s_k^2$ から $\sqrt{\sum_kp_k^2s_k^2/n_k}$ で推定します（各層に 2 個以上の標本が要ります）。 の積分で、$[0,1)$ を $K$ 個の等しい区間 $A_k = [\frac{k-1}K, \frac kK)$ に分けます。$p_k = 1/K$ で、層 $k$ の標本は $X^{(k)} = (k - 1 + V)/K$（$V \sim \mathrm{U}(0,1)$）です。$e^{X^{(k)}} = e^{(k-1)/K}e^{V/K}$ なので $$ \sigma_k^2 = e^{2(k-1)/K}v_K, \qquad v_K = \mathrm{Var}(e^{V/K}) = \frac{K(e^{2/K} - 1)}{2} - K^2(e^{1/K} - 1)^2 $$ です。等比級数の和から $$ \sum_{k=1}^Kp_k\sigma_k^2 = \frac{v_K}{K}\cdot\frac{e^2 - 1}{e^{2/K} - 1} $$ で、$K = 10$ では $v_{10} \approx 9.216\times10^{-4}$、$\sum_kp_k\sigma_k^2 \approx 0.002659$ です。$\sigma^2 \approx 0.24204$ の $0.0110$ 倍で、分散は約 $1/91$ になります。層の幅を $1/K$ にすると層内分散は $1/K^2$ に比例して減るので（滑らかな $f$ では $\sigma_k^2 \approx f'(x)^2/(12K^2)$）、1 次元では層を細かくするほど効きます。$d$ 次元で各軸を $K$ 等分すると層の数が $K^d$ になり、高次元ではすべての変数を層別することはできません。重要な少数の変数だけを層別します。 次のコードは、 の 4 つの方法を $N = 10^4$ で $200$ 回ずつ独立に繰り返し、推定量の分散の比を閉じた形の理論値と比べます。 ```python # Python 3.9.6, NumPy 2.0.2 import numpy as np e = np.e I = e - 1 s2 = (e**2 - 1) / 2 - (e - 1) ** 2 # Var(e^U) cov_anti = e - (e - 1) ** 2 # Cov(e^U, e^{1-U}) cov_cv = 1 - (e - 1) / 2 # Cov(e^U, U) K = 10 vK = K * (np.exp(2 / K) - 1) / 2 - (K * (np.exp(1 / K) - 1)) ** 2 # Var(e^{V/K}), V ~ U(0,1) within = vK * (e**2 - 1) / (K * (np.exp(2 / K) - 1)) # sum_k p_k sigma_k^2 theory = {"単純": 1.0, "対称変量": (s2 + cov_anti) / s2, "制御変量": 1 - cov_cv**2 / (s2 / 12), "層別(10層)": within / s2} rng = np.random.default_rng(20261008) N, R = 10**4, 200 est = {k: np.empty(R) for k in theory} beta = 12 * cov_cv # 最適な係数 beta*（既知として使う） for r in range(R): u = rng.random(N) est["単純"][r] = np.exp(u).mean() h = u[: N // 2] est["対称変量"][r] = 0.5 * (np.exp(h) + np.exp(1 - h)).mean() est["制御変量"][r] = (np.exp(u) - beta * (u - 0.5)).mean() x = (np.arange(N) % K + u) / K # 各層に N/K 個（比例配分） est["層別(10層)"][r] = np.exp(x).mean() v0 = est["単純"].var(ddof=1) rel = np.sqrt(4 / (R - 1)) # 分散の比の相対標準誤差（正規近似） print(f"sigma^2 = {s2:.6f}, beta* = {beta:.6f}") for k in theory: m, v = est[k].mean(), est[k].var(ddof=1) ratio = v / v0 print(f"{k:10s} 平均の偏り = {(m - I)/np.sqrt(v/R):+.2f} SE 分散比 = {ratio:.5f} ± {ratio*rel:.5f} 理論 = {theory[k]:.5f}") ``` ```text sigma^2 = 0.242036, beta* = 1.690309 単純 平均の偏り = +0.69 SE 分散比 = 1.00000 ± 0.14178 理論 = 1.00000 対称変量 平均の偏り = +1.10 SE 分散比 = 0.02794 ± 0.00396 理論 = 0.03233 制御変量 平均の偏り = -0.16 SE 分散比 = 0.01512 ± 0.00214 理論 = 0.01628 層別(10層) 平均の偏り = +0.27 SE 分散比 = 0.01074 ± 0.00152 理論 = 0.01099 ``` 4 つの推定量とも、200 回の平均は真値から標準誤差の 1.1 倍以内にあり（不偏性と矛盾しません）、分散の比は理論値と標準誤差の 1.2 倍以内で一致しています。分散の比の標準誤差は、200 個の標本から求めた 2 つの分散の比のばらつきを正規近似 $\sqrt{4/(R-1)}$ で見積もったものです。 ## 6. 重点サンプリング ### 6.1. 提案分布と重み これまでの方法は、目標の分布 $p$ から標本を取っていました。重点サンプリングは、**別の分布 $q$ から標本を取り、[密度](/computer-science/stochastic-computation/sampling-from-distributions#def-density-vector)の比で補正します**。被積分関数が大きい領域に標本を集めるのが目的です。 $p$ を $\mathbb{R}^d$ 上の確率密度、$f$ を $\int|f|p < \infty$ を満たす可測関数、$I = \int f(x)p(x)\,dx$ とする。確率密度 $q$ が、$f(x)p(x) \ne 0$ となるすべての $x$ で $q(x) > 0$ を満たすとき、$q$ を**提案密度**という。 $$ w(x) = \frac{p(x)}{q(x)} \quad (q(x) > 0), \qquad w(x) = 0 \quad (q(x) = 0) $$ を**重み**という。密度 $q$ に従う独立な $X_1, \ldots, X_N$ に対して $$ \hat I_N^{\mathrm{is}} = \frac1N\sum_{i=1}^N f(X_i)\,w(X_i) $$ を**重点サンプリング推定量**という。以下 $\mathbb{E}_q$、$\mathrm{Var}_q$ は $X \sim q$ についての期待値と分散を表す。 の設定で次が成り立つ。 1. $\mathbb{E}_q[f(X)w(X)] = I$。したがって $\hat I_N^{\mathrm{is}}$ は不偏である。 2. $\mathrm{Var}\bigl(\hat I_N^{\mathrm{is}}\bigr) = \dfrac1N\Bigl(\displaystyle\int_{\{q > 0\}}\frac{f(x)^2p(x)^2}{q(x)}\,dx - I^2\Bigr)$（右辺の積分が $+\infty$ のときは分散も $+\infty$）。 3. $c = \int|f|p > 0$ とする。すべての提案密度 $q$ について $\int_{\{q>0\}}f^2p^2/q \ge c^2$ で、$q^* = |f|p/c$ のとき等号が成り立つ。すなわち分散の最小値は $(c^2 - I^2)/N$ である。特に $f \ge 0$ ならば $c = I$ で、$q^* = fp/I$ による推定量の分散は $0$ である。 **1.** $\{q > 0\}$ の上で $fwq = fp$ です。$\{q = 0\}$ の上では提案密度の条件から $fp = 0$ です。よって $$ \mathbb{E}_q[f(X)w(X)] = \int_{\{q>0\}} f(x)\frac{p(x)}{q(x)}q(x)\,dx = \int_{\{q>0\}}f(x)p(x)\,dx = \int f(x)p(x)\,dx = I $$ です。不偏性は の 1 です。 **2.** 同じ計算で $\mathbb{E}_q[(fw)^2] = \int_{\{q>0\}}f^2p^2/q$ です。 の 2 と $\mathrm{Var} = \mathbb{E}[(\cdot)^2] - I^2$ から従います。 **3.** $\{q = 0\}$ の上で $|f|p = 0$ なので、積分についての Cauchy–Schwarz の不等式から $$ c = \int_{\{q>0\}}\frac{|f|p}{\sqrt q}\cdot\sqrt q\,dx \le \Bigl(\int_{\{q>0\}}\frac{f^2p^2}{q}\,dx\Bigr)^{1/2}\Bigl(\int q\,dx\Bigr)^{1/2} = \Bigl(\int_{\{q>0\}}\frac{f^2p^2}{q}\,dx\Bigr)^{1/2} $$ です。$q^* = |f|p/c$ は確率密度で、$fp \ne 0$ のところで正なので提案密度です。このとき $\{q^* > 0\}$ の上で $f^2p^2/q^* = c|f|p$ なので、積分は $c\cdot c = c^2$ です。$f \ge 0$ なら $c = \int fp = I$ です。 3 は目標を示しますが、そのままでは使えません。$q^* = fp/I$ を書くには、求めたい $I$ 自身が要るからです。この定理の読み方は、「**$q$ を $|f|p$ の形に近づけるほど分散が減る**」です。$f$ が狭い領域でだけ大きい（稀な事象の確率、鋭いピークをもつ尤度）とき、$p$ からの標本はほとんどその領域に入らず、単純な推定量の相対誤差は非常に大きくなります。その領域に標本を集める $q$ を使うと、分散が何桁も減ります。 逆に、$q$ の選び方がまずいと分散は増え、無限大にもなります。2 の被積分関数 $f^2p^2/q$ は、$q$ が $p$ より速く $0$ に近づく領域で大きくなります。 目標を $p = N(0, \sigma_p^2)$、提案を $q = N(0, \sigma_q^2)$ とします。重みは $w(x) = \frac{\sigma_q}{\sigma_p}\exp\bigl(-\frac{x^2}{2}(\frac1{\sigma_p^2} - \frac1{\sigma_q^2})\bigr)$ で、 $$ \mathbb{E}_q[w^2] = \int\frac{p^2}{q}\,dx = \frac{\sigma_q}{\sqrt{2\pi}\,\sigma_p^2}\int_{-\infty}^\infty\exp\Bigl(-x^2\Bigl(\frac{1}{\sigma_p^2} - \frac{1}{2\sigma_q^2}\Bigr)\Bigr)dx $$ です。指数の係数 $a = \frac1{\sigma_p^2} - \frac1{2\sigma_q^2}$ が正のとき、すなわち $$ \sigma_q^2 > \frac{\sigma_p^2}{2} $$ のときに限り積分は収束し、Gauss 積分 $\int e^{-ax^2}dx = \sqrt{\pi/a}$（ で $x \mapsto x/\sqrt a$）から $$ \mathbb{E}_q[w^2] = \frac{\sigma_q}{\sqrt{2\pi}\,\sigma_p^2}\sqrt{\frac{2\pi\sigma_p^2\sigma_q^2}{2\sigma_q^2 - \sigma_p^2}} = \frac{\sigma_q^2}{\sigma_p\sqrt{2\sigma_q^2 - \sigma_p^2}} $$ です。$\sigma_q^2 \le \sigma_p^2/2$ では $\mathbb{E}_q[w^2] = \infty$ で、$f \equiv 1$（あるいは裾で $0$ に近づかない任意の $f$）の推定量の分散は無限大です。このとき $\hat I_N^{\mathrm{is}}$ は大数の法則で収束はしますが、ときどき現れる巨大な重みが標本平均を跳ね上げ、$s_N$ による誤差の見積もりは当てになりません。提案の分散を目標より大きめに取る（$\sigma_q = 1.5\sigma_p$ なら $\mathbb{E}_q[w^2] = 2.25/\sqrt{3.5} \approx 1.203$）のが安全な側です。 ```python # Python 3.9.6, NumPy 2.0.2 import numpy as np def weights(sq, n, rng): # 目標 N(0,1)、提案 N(0, sq^2) x = sq * rng.standard_normal(n) return sq * np.exp(-0.5 * x**2 * (1 - 1 / sq**2)) for sq in (1.5, 0.6): rng = np.random.default_rng(20261008) cond = 2 * sq**2 > 1 ew2 = sq**2 / np.sqrt(2 * sq**2 - 1) if cond else np.inf print(f"sigma_q = {sq}: E_q[w^2] の理論値 = {ew2:.4f}") for n in (10**3, 10**4, 10**5, 10**6, 10**7): w = weights(sq, n, rng) print(f" N = {n:>8d} 平均 = {w.mean():.4f} 二乗平均 = {(w**2).mean():8.4f} 最大の重み = {w.max():9.2f}") ``` ```text sigma_q = 1.5: E_q[w^2] の理論値 = 1.2027 N = 1000 平均 = 0.9962 二乗平均 = 1.1928 最大の重み = 1.50 N = 10000 平均 = 1.0026 二乗平均 = 1.2066 最大の重み = 1.50 N = 100000 平均 = 0.9978 二乗平均 = 1.1979 最大の重み = 1.50 N = 1000000 平均 = 1.0002 二乗平均 = 1.2028 最大の重み = 1.50 N = 10000000 平均 = 1.0000 二乗平均 = 1.2027 最大の重み = 1.50 sigma_q = 0.6: E_q[w^2] の理論値 = inf N = 1000 平均 = 0.9827 二乗平均 = 2.4232 最大の重み = 26.76 N = 10000 平均 = 0.9779 二乗平均 = 2.8601 最大の重み = 52.37 N = 100000 平均 = 0.9920 二乗平均 = 7.0902 最大の重み = 510.61 N = 1000000 平均 = 0.9973 二乗平均 = 9.3320 最大の重み = 1163.03 N = 10000000 平均 = 0.9999 二乗平均 = 14.9321 最大の重み = 4688.85 ``` $\sigma_q = 1.5$ では重みは $1.5$ 以下に有界で、二乗平均は理論値 $1.2027$ に落ち着きます。$\sigma_q = 0.6$（$\sigma_q^2 = 0.36 < 0.5$）では、重みの平均は $1$ に近づいていきますが、二乗平均は $N$ を増やすたびに大きくなり、最大の重みも増え続けます。この標本から計算した標準誤差は、$N$ をいくら増やしても信用できません。 $Z \sim N(0,1)$ について $I = \mathbb{P}(Z > 4) \approx 3.167\times10^{-5}$ を推定します。単純な推定量 $\mathbf{1}\{Z > 4\}$ の分散は $I(1 - I)$ で、相対標準誤差は $\sqrt{(1 - I)/(NI)} \approx 178/\sqrt N$ です。$N = 10^6$ でも $18\%$ です。 提案を平行移動した正規分布 $q = N(\mu, 1)$ とします。重みは $w(x) = \varphi(x)/\varphi(x - \mu) = e^{-\mu x + \mu^2/2}$ で、 の 2 の積分は $$ \int_4^\infty\frac{\varphi(x)^2}{\varphi(x - \mu)}\,dx = \int_4^\infty\varphi(x)\,e^{-\mu x + \mu^2/2}\,dx = e^{\mu^2}\int_4^\infty\varphi(x + \mu)\,dx = e^{\mu^2}\bigl(1 - \Phi(4 + \mu)\bigr) $$ です（$\varphi(x)e^{-\mu x} = \varphi(x + \mu)e^{\mu^2/2}$ を使いました）。$\mu = 4$ では $e^{16}(1 - \Phi(8)) \approx 5.528\times10^{-9}$ で、1 標本あたりの分散は $5.528\times10^{-9} - I^2 \approx 4.525\times10^{-9}$ です。単純な推定量の分散 $3.167\times10^{-5}$ の約 $1/7000$ で、相対標準誤差は $2.12/\sqrt N$ になります。 ```python # Python 3.9.6, NumPy 2.0.2, SciPy 1.13.1, mpmath 1.3.0 import numpy as np from scipy.special import erfc from mpmath import mp, erfc as mp_erfc, sqrt as mp_sqrt, exp as mp_exp mp.dps = 30 p_ref = float(mp_erfc(4 / mp_sqrt(2)) / 2) print(f"参照値 P(Z > 4): mpmath = {p_ref:.12e}, scipy = {0.5 * erfc(4 / np.sqrt(2)):.12e}") rng = np.random.default_rng(20261008) N, mu = 10**6, 4.0 z = rng.standard_normal(N) y = (z > 4).astype(float) # 単純な推定 x = z + mu # 提案 N(4, 1) からの標本 w = np.exp(-mu * x + mu**2 / 2) # 重み p(x)/q(x) g = (x > 4) * w for name, v in (("単純", y), ("重点サンプリング", g)): est, se = v.mean(), v.std(ddof=1) / np.sqrt(N) print(f"{name:10s} 推定値 = {est:.4e} ± {se:.1e} 相対標準誤差 = {se/est:.2%} 真値との差 = {(est - p_ref)/se:+.2f} SE") var_is = float(mp_exp(16) * mp_erfc(8 / mp_sqrt(2)) / 2) - p_ref**2 # 1 標本あたりの分散（閉じた形） print(f"1 標本あたりの分散: 単純 = {p_ref*(1-p_ref):.3e}, 重点 = {var_is:.3e}, 比 = {p_ref*(1-p_ref)/var_is:.0f}") print(f"重みの有効標本サイズ / N = {w.sum()**2 / (w**2).sum() / N:.1e}") ``` ```text 参照値 P(Z > 4): mpmath = 3.167124183312e-05, scipy = 3.167124183312e-05 単純 推定値 = 3.2000e-05 ± 5.7e-06 相対標準誤差 = 17.68% 真値との差 = +0.06 SE 重点サンプリング 推定値 = 3.1664e-05 ± 6.7e-08 相対標準誤差 = 0.21% 真値との差 = -0.10 SE 1 標本あたりの分散: 単純 = 3.167e-05, 重点 = 4.525e-09, 比 = 6999 重みの有効標本サイズ / N = 1.2e-04 ``` 同じ $N = 10^6$ で、相対標準誤差は $17.7\%$ から $0.21\%$ になりました（比は $\sqrt{6999} \approx 84$）。どちらの推定値も参照値から標準誤差の $0.1$ 倍以内です。最後の行の「有効標本サイズ」は次の節で定義します。この値が極端に小さいのに推定がうまくいっている理由もそこで説明します。 ### 6.2. 自己正規化推定量と有効標本サイズ Bayes 統計の事後分布や Boltzmann 分布では、目標の密度が $p = \tilde p/Z$ の形で、正規化定数 $Z = \int\tilde p$ が未知です。このとき重み $p/q$ は計算できませんが、比の形にすれば $Z$ が打ち消されます。 $\tilde p \ge 0$ を $0 < Z = \int\tilde p < \infty$ を満たす可測関数、$p = \tilde p/Z$、$q$ を $\tilde p(x) > 0$ となるすべての $x$ で $q(x) > 0$ を満たす確率密度とし、$\tilde w = \tilde p/q$（$q = 0$ のところでは $0$）と置く。密度 $q$ に従う独立な $X_1, \ldots, X_N$ に対して $$ \tilde I_N = \frac{\sum_{i=1}^N\tilde w(X_i)f(X_i)}{\sum_{i=1}^N\tilde w(X_i)} $$ を**自己正規化重点サンプリング推定量**という（分母が $0$ のときは $0$ と定める）。また $$ N_{\mathrm{eff}} = \frac{\bigl(\sum_i\tilde w(X_i)\bigr)^2}{\sum_i\tilde w(X_i)^2} $$ を重みの**有効標本サイズ**という。 の設定で、$\int|f|p < \infty$、$I = \int fp$ とする。 1. $N \to \infty$ で $\tilde I_N \to I$（概収束）。 2. $\mathbb{E}_q[\tilde w^2(1 + f^2)] < \infty$ かつ $\tau^2 := \mathbb{E}_q[\tilde w^2(f - I)^2]/Z^2 > 0$ ならば、$\sqrt N(\tilde I_N - I)$ は $N(0, \tau^2)$ に分布収束する。 3. $f$ が有界（$|f| \le M_f$）で、$\mathbb{E}_q[\tilde w^2] < \infty$ かつ $\mathbb{P}_q(\tilde w(X) > 0) = 1$ ならば $$ \mathbb{E}[\tilde I_N] - I = -\frac{\mathbb{E}_q[\tilde w^2(f - I)]}{NZ^2} + R_N, \qquad |R_N| \le \frac{2M_f\,\mathrm{Var}_q(\tilde w)}{NZ^2}. $$ 特に偏りは $O(1/N)$ である。 4. $1 \le N_{\mathrm{eff}} \le N$（分母が正のとき）であり、$\mathbb{E}_q[\tilde w^2] < \infty$ ならば $N_{\mathrm{eff}}/N \to Z^2/\mathbb{E}_q[\tilde w^2] = 1/\mathbb{E}_q[w^2]$（概収束）。ここで $w = \tilde w/Z = p/q$。 次の証明のうち 3（偏りの評価）は、推定量が「分子の平均 ÷ 分母の平均」という比の形で、比の期待値は期待値の比にならないために少し手間がかかります。分母の逆数を、分母が $Z$ に近いことを使って $Z$ のまわりで展開し、各項の期待値を評価する、という方針です。 $a_i = \tilde w(X_i)(f(X_i) - I)$、$b_i = \tilde w(X_i)$ と置き、$A_N = \frac1N\sum a_i$、$B_N = \frac1N\sum b_i$ とします。 の 1 と同じ計算で $$ \mathbb{E}_q[b_1] = \int\tilde p = Z, \qquad \mathbb{E}_q[a_1] = \int\tilde p\,(f - I) = Z(I - I) = 0 $$ です。$B_N > 0$ のとき $\tilde I_N - I = \frac{\sum b_if(X_i)}{\sum b_i} - I = A_N/B_N$ です。 **1.** 大数の強法則（）から $A_N \to 0$、$B_N \to Z > 0$（概収束）です。よって確率 $1$ で、十分大きい $N$ では $B_N > 0$ で、$\tilde I_N - I = A_N/B_N \to 0$ です。 **2.** 仮定から $a_1 \in L^2$ で、$\mathrm{Var}_q(a_1) = \mathbb{E}_q[\tilde w^2(f - I)^2] = Z^2\tau^2 > 0$ です。中心極限定理（）から $\sqrt NA_N/(Z\tau)$ は標準正規分布に分布収束します。$B_N/Z \to 1$（概収束、したがって確率収束）なので、Slutsky の定理（）から $\sqrt NA_N/(\tau B_N)$ も標準正規分布に分布収束します。 の 4 の証明と同じく、$B_N = 0$ となる確率は $0$ に収束するので、$\sqrt N(\tilde I_N - I)/\tau$ も同じ極限をもちます。 **3.** 仮定から確率 $1$ で $B_N > 0$ です。$\tilde I_N$ は $f(X_i)$ の重みつき平均なので $|\tilde I_N| \le M_f$、また $|I| \le M_f$ から $|A_N/B_N| \le 2M_f$ です。恒等式 $$ \frac{1}{B} = \frac1Z - \frac{B - Z}{Z^2} + \frac{(B - Z)^2}{Z^2B} \qquad (B > 0) $$ （右辺を通分すると $\frac{ZB - B(B-Z) + (B-Z)^2}{Z^2B} = \frac{Z^2}{Z^2B}$）に $A_N$ を掛けて期待値を取ります。$\mathbb{E}[A_N] = 0$ です。$a_i$ と $b_j$（$i \ne j$）は独立で $\mathbb{E}[a_i] = 0$ なので、 $$ \mathbb{E}[A_N(B_N - Z)] = \frac{1}{N^2}\sum_i\mathbb{E}[a_i(b_i - Z)] = \frac1N\mathbb{E}_q[a_1b_1] = \frac1N\mathbb{E}_q[\tilde w^2(f - I)] $$ です。第 3 項は、$|A_N/B_N| \le 2M_f$ から絶対値が $2M_f(B_N - Z)^2/Z^2$ 以下で、その期待値は $2M_f\,\mathrm{Var}(B_N)/Z^2 = 2M_f\,\mathrm{Var}_q(\tilde w)/(NZ^2)$ 以下です。以上を合わせて主張の式を得ます。 **4.** $b_i \ge 0$ なので $(\sum b_i)^2 \ge \sum b_i^2$、また Cauchy–Schwarz の不等式から $(\sum b_i\cdot1)^2 \le N\sum b_i^2$ です。よって $1 \le N_{\mathrm{eff}} \le N$ です。$N_{\mathrm{eff}}/N = B_N^2/(\frac1N\sum b_i^2)$ で、大数の強法則から分子は $Z^2$、分母は $\mathbb{E}_q[\tilde w^2]$ に概収束します。 2 は、比 $A_N/B_N$ に対する delta 法（ の 2 変数版）の結論を、Slutsky の定理で直接導いたものです。漸近分散は、正規化した重み $w = p/q$ で書くと $\tau^2 = \mathbb{E}_q[w^2(f - I)^2]$ で、標準誤差は $$ \widehat{\mathrm{se}} = \frac{\sqrt{\sum_i\tilde w_i^2(f_i - \tilde I_N)^2}}{\sum_i\tilde w_i} $$ で推定します。3 から、自己正規化推定量は不偏ではありませんが、偏り $O(1/N)$ は標準誤差 $O(N^{-1/2})$ より高次で、$N$ が大きければ無視できます。 4 の有効標本サイズは、重みのばらつきの指標です。すべての重みが等しければ $N_{\mathrm{eff}} = N$、1 個の重みだけが大きければ $N_{\mathrm{eff}} \approx 1$ です。$(f - I)^2$ と $w^2$ が $q$ のもとでほぼ無相関なら $\tau^2 \approx \mathbb{E}_q[w^2]\,\mathbb{E}_q[(f-I)^2]$ となり、「重みのせいで標本が $N/\mathbb{E}_q[w^2] \approx N_{\mathrm{eff}}$ 個分の働きしかしない」と読めます。ただし $N_{\mathrm{eff}}$ は $f$ を見ない診断です。 では $\mathbb{E}_q[w^2] = e^{16}$ で、$N_{\mathrm{eff}}/N$ の極限は $e^{-16} \approx 1.1\times10^{-7}$ と絶望的に見えますが、巨大な重みが現れるのは $x < 4$ の側で、そこでは $f = \mathbf{1}\{x > 4\}$ が $0$ なので、推定には影響しません（出力の $1.2\times10^{-4}$ は、$10^6$ 個の標本ではこの極限にすら届いていないことも示しています）。特定の $f$ のための提案分布は、ほかの $f$ には使えないということです。事後分布のように多くの $f$ の期待値を同じ標本から求める場合には、$N_{\mathrm{eff}}$ が $N$ に比べて極端に小さくないことを確かめます。 - **多重レベル Monte Carlo 法**（Giles）は、離散化の細かさが違う近似の差を制御変量のように積み重ね、粗い近似に多くの標本を、細かい近似に少ない標本を割り当てます。確率微分方程式の期待値の計算量を大きく減らします。 - **乱数の共通化**は、2 つの設定の差 $\mathbb{E}[f_1] - \mathbb{E}[f_2]$ を推定するとき、両方に同じ乱数を使って正の相関を作り、差の分散を減らします。対称変量と逆向きの発想です。 - **稀な事象の分割法**は、事象に近づいた標本を複製して、稀な領域の標本を増やします。 これらは Glasserman の本と Asmussen–Glynn の本に詳しく書かれています。この章では定理としては扱いません。 ## 7. 演習 の推定量 $\hat\pi_N$（標準偏差 $1.642/\sqrt N$）で、$|\hat\pi_N - \pi| \le 10^{-4}$ を確率 $95\%$ で達成するのに必要な $N$ を、(1) 正規近似（ の 4）、(2) Chebyshev の不等式（ の 1）のそれぞれで求めよ。 の制御変量つきの推定量を使うと (1) はどう変わるか。 (1) $1.96 \times 1.642/\sqrt N \le 10^{-4}$ から $N \ge (1.96 \times 1.642 \times 10^4)^2 \approx 1.04\times10^9$。 (2) $1/k^2 = 0.05$ から $k = \sqrt{20} \approx 4.472$ で、$N \ge (4.472 \times 1.642 \times 10^4)^2 \approx 5.39\times10^9$。正規近似の $20/1.96^2 \approx 5.2$ 倍です。 制御変量つきの推定量は 1 標本あたりの分散が約 $1/104$ なので、(1) は $1.04\times10^9/104 \approx 1.0\times10^7$ になります。必要な標本数は分散に比例するので、分散低減の倍率がそのまま標本数の節約になります。 $U \sim \mathrm{U}(0,1)$、$f(u) = (u - \frac12)^2$ とする。対称変量による推定量の分散を単純な推定量の分散と比べよ。また、$f(u) = u^2$ ではどうなるか。 $f(u) = (u - \frac12)^2$ は $f(1 - u) = f(u)$ を満たすので、$\mathrm{Cov}(f(U), f(1-U)) = \mathrm{Var}(f(U)) = \sigma^2$ です。 の 2 の等式（単調性は不等号にだけ使われています）から分散は $2\sigma^2/N$ で、単純な推定量の 2 倍です。$f$ が単調でないので、定理の不等号は保証されません。 $f(u) = u^2$ は単調増加です。$\mathbb{E}[U^2] = \frac13$、$\sigma^2 = \frac15 - \frac19 = \frac4{45}$、$\mathbb{E}[U^2(1-U)^2] = \frac13 - \frac24 + \frac15 = \frac1{30}$ から $\mathrm{Cov} = \frac1{30} - \frac19 = -\frac7{90}$ です。分散は $(\frac4{45} - \frac7{90})/N = \frac{1}{90N}$ で、単純な推定量の $\frac1{90}/\frac4{45} = \frac18$ 倍です。 を一般化して、$\mathbb{P}(Z > a)$（$a > 0$）を提案 $N(\mu, 1)$ の重点サンプリングで推定する。(1) 推定量の 2 次モーメントが $m(\mu) = e^{\mu^2}(1 - \Phi(a + \mu))$ であることを示せ。(2) $m$ を最小にする $\mu^*$ は $2\mu^*R(a + \mu^*) = 1$ を満たすことを示せ。ここで $R(x) = (1 - \Phi(x))/\varphi(x)$ は Mills 比である。(3) $\mu^* > a$ を示せ。 (1) の計算で $4$ を $a$ に替えるだけです。$\int_a^\infty\varphi(x)^2/\varphi(x-\mu)\,dx = e^{\mu^2}\int_a^\infty\varphi(x + \mu)\,dx = e^{\mu^2}(1 - \Phi(a + \mu))$。 (2) $m'(\mu) = e^{\mu^2}\bigl(2\mu(1 - \Phi(a + \mu)) - \varphi(a + \mu)\bigr) = e^{\mu^2}\varphi(a+\mu)\bigl(2\mu R(a + \mu) - 1\bigr)$ です。$\mu \le 0$ では括弧の中が負で $m$ は減少するので、最小点は $\mu > 0$ にあり、そこで $m'(\mu^*) = 0$、すなわち $2\mu^*R(a + \mu^*) = 1$ です（$\mu \to \infty$ で $m(\mu) \to \infty$ となることは、Mills 比の下界 $R(x) > x/(1 + x^2)$ から $m(\mu) > e^{\mu^2}\varphi(a+\mu)\frac{a+\mu}{1 + (a+\mu)^2}$ で、指数部分が $\mu^2 - (a+\mu)^2/2 \to \infty$ となることから分かります。よって最小点は存在します）。 (3) Mills 比の上界 $R(x) < 1/x$（$x > 0$。 の 2）を $x = a + \mu^*$ に使うと、$1 = 2\mu^*R(a + \mu^*) < 2\mu^*/(a + \mu^*)$ から $a + \mu^* < 2\mu^*$、すなわち $\mu^* > a$ です。最適な平行移動は、境界 $a$ より少し奥です。$a = 4$ では数値的に $\mu^* \approx 4.12$ で、2 次モーメントは $\mu = 4$ の $5.528\times10^{-9}$ から $5.488\times10^{-9}$ へ $0.7\%$ 減るだけなので、実用上は $\mu = a$ で十分です。 $X \sim \mathrm{U}(0,1)$、$f(x) = \mathbf{1}\{x > 0.9\}\cdot 10(x - 0.9)$ とし、層を $A_1 = [0, 0.9]$、$A_2 = (0.9, 1)$ とする。$\sigma^2$、比例配分の分散、Neyman 配分とその分散を求めよ。 $A_1$ の上で $f = 0$ なので $\mu_1 = 0$、$\sigma_1 = 0$。$A_2$ の上では $f(X^{(2)}) = V$（$V \sim \mathrm{U}(0,1)$）なので $\mu_2 = \frac12$、$\sigma_2^2 = \frac1{12}$。$p_1 = 0.9$、$p_2 = 0.1$ です。 - $I = 0.1 \times \frac12 = 0.05$。$\mathbb{E}[f^2] = 0.1 \times \frac13$ から $\sigma^2 = \frac1{30} - 0.0025 \approx 0.03083$。 - 比例配分: の 3 から $\frac1N\sum p_k\sigma_k^2 = \frac{0.1}{12N} \approx 0.00833/N$。単純な推定量の $0.27$ 倍。確認として層間分散は $0.9(0 - 0.05)^2 + 0.1(0.5 - 0.05)^2 = 0.0225$ で、$0.00833 + 0.0225 = 0.03083$ です。 - Neyman 配分: $n_k \propto p_k\sigma_k$ から $n_1 = 0$、$n_2 = N$（$\sigma_1 = 0$ の層には、$\mu_1 = 0$ が既知なので標本が要りません）。分散は $\frac1N(\sum p_k\sigma_k)^2 = \frac{0.01}{12N} \approx 0.000833/N$ で、比例配分のさらに $1/10$、単純な推定量の $0.027$ 倍です。 Neyman 配分は、ここでは「$f$ が $0$ でない領域にだけ標本を置く」という重点サンプリングと同じことをしています。実際、$q$ を $A_2$ 上の一様分布とした重点サンプリング推定量は $0.1 \times V$ の標本平均で、分散は同じ $\frac{0.01}{12N}$ です。 ## 参考文献 1. A. B. Owen, *Monte Carlo theory, methods and examples*, 2013 — 著者が公開している本（[artowen.su.domains/mc](https://artowen.su.domains/mc/)）。第 8 章（分散低減）、第 9 章（重点サンプリング）。 2. P. Glasserman, *Monte Carlo Methods in Financial Engineering*, Springer, 2004 — 第 4 章（分散低減法）。 3. S. Asmussen, P. W. Glynn, *Stochastic Simulation: Algorithms and Analysis*, Springer, 2007 — 第 III 章（出力の解析）、第 V 章（分散低減）、第 VI 章（稀な事象）。 4. C. P. Robert, G. Casella, *Monte Carlo Statistical Methods*, 2nd ed., Springer, 2004 — 第 3 章（Monte Carlo 積分）、第 4 章（Monte Carlo 法の分散の制御）。 5. J. M. Hammersley, D. C. Handscomb, *Monte Carlo Methods*, Methuen, 1964 — 対称変量・制御変量の古典的な解説。 6. E. Novak, *Deterministic and Stochastic Error Bounds in Numerical Analysis*, Lecture Notes in Mathematics 1349, Springer, 1988 — 決定的な求積則の下界と Monte Carlo 法の比較。 7. M. B. Giles, "Multilevel Monte Carlo path simulation", *Operations Research* 56 (2008), 607–617.