対照変量法（アンチセティック）

🎓 レベル：標準　|　重要度：A（必須）

📎 前提：収束率と誤差（√n則）　|　関連：制御変量法

要点（BLUF）

対照変量法：一様乱数 $U$ とその鏡像 $1-U$ のように負の相関を持つペアで関数を評価し、ペアの平均を推定に使います。
二つの誤差が逆方向に振れて相殺し、相関 $\rho$ が負の分だけ分散が下がる（理論： $1/(1+\rho)$ 倍）。
$\int_0^1 e^x dx$ で、同じ計算量のまま分散が約32倍縮小（ $e^U$ と $e^{1-U}$ の相関 $-0.97$ ）。

1. アイデア：誤差を逆向きに振らせる

モンテカルロ推定 $\frac{1}{n}\sum g(U_i)$ の分散は、各項の分散 $\sigma^2$ で決まります。もし負に相関する2つの推定を平均できれば、片方が上振れたとき他方が下振れて相殺し、平均の分散が下がります。

一様乱数 $U \sim \mathcal{U}(0,1)$ に対し $1-U$ も $\mathcal{U}(0,1)$ で、しかも $U$ が大きいとき $1-U$ は小さい（完全な負の依存）。 $g$ が単調なら $g(U)$ と $g(1-U)$ も負に相関します。そこでペアの平均

\hat{I}_{\text{anti}} = \frac{1}{n}\sum_{i=1}^{n} \frac{g(U_i) + g(1-U_i)}{2}

を推定量にします。

2. なぜ分散が下がるか

ペア1組の分散は

\text{Var}\!\left(\frac{g(U)+g(1-U)}{2}\right) = \frac{1}{4}\Big(\sigma^2 + \sigma^2 + 2\,\text{Cov}\big(g(U),g(1-U)\big)\Big) = \frac{\sigma^2}{2}\,(1+\rho)

ここで $\rho = \text{Corr}(g(U), g(1-U))$ 。同じ関数評価回数（プレーンは $2n$ 回、対照は $n$ ペア＝ $2n$ 回）で比べると、分散比は

\frac{\text{Var(対照)}}{\text{Var(プレーン)}} = 1 + \rho

$\rho < 0$ なら $1+\rho < 1$ で分散減少。 $\rho = -0.97$ なら $1+\rho = 0.03$ 、すなわち約33倍の改善です。 $g$ が単調なほど $\rho$ は強く負になり、効果が大きくなります。

3. 実測

import numpy as np

# 乱数シードを固定
rng = np.random.default_rng(30)
n = 100_000
reps = 2000

def plain(rng, m):                       # プレーンMC（2n標本）
    u = rng.random(m)
    return np.exp(u)

def antithetic(rng, m):                  # 対照変量（mペア=2m評価）
    u = rng.random(m)
    return (np.exp(u) + np.exp(1 - u)) / 2

plain_ests = np.array([plain(rng, 2*n).mean() for _ in range(reps)])
anti_ests  = np.array([antithetic(rng, n).mean() for _ in range(reps)])

print(f"プレーン   平均={plain_ests.mean():.5f}  分散={plain_ests.var():.2e}")
print(f"対照変量   平均={anti_ests.mean():.5f}  分散={anti_ests.var():.2e}")
print(f"分散減少率 = {plain_ests.var()/anti_ests.var():.2f} 倍")

u = rng.random(500_000)
print(f"corr(e^U, e^(1-U)) = {np.corrcoef(np.exp(u), np.exp(1-u))[0,1]:.4f}")

出力：

プレーン   平均=1.71829  分散=1.29e-06
対照変量   平均=1.71829  分散=4.02e-08
分散減少率 = 31.99 倍

corr(e^U, e^(1-U)) = -0.9677

出力の意味：両者とも真値 $e-1 = 1.71828$ を正しく推定（不偏性は保たれる）。しかし分散はプレーン $1.29\times10^{-6}$ に対し対照変量 $4.02\times10^{-8}$ で約32倍の改善。理論の $1/(1+\rho) = 1/(1-0.9677) = 31.0$ とよく一致します。標準誤差にすると $\sqrt{32} \approx 5.7$ 倍精度が上がる＝同じ精度をプレーンの約1/32のサンプルで達成できる、ということです。

数式の直観的意味

対照変量は「サンプルをわざと裏表ペアで取り、運の偏りを打ち消す」操作です。たまたま $U$ が小さい値ばかり出た（推定が下振れする）回でも、 $1-U$ は大きい値ばかりになり上振れ、平均すれば偏りが消える。 $1+\rho$ という分散比は「ペアの誤差がどれだけ逆向きか」の指標で、完全逆相関 $\rho=-1$ なら分散ゼロ（線形関数なら実際に厳密化）。鍵は** $g$ の単調性**——単調なら $U \uparrow$ で $g(U)\uparrow$ 、 $g(1-U)\downarrow$ となり負相関が保証されます。非単調な $g$ では $\rho$ が正になり逆効果になり得ます。

⚠️ よくある誤解・落とし穴

「いつでも分散が下がる」ではない： $g$ が非単調だと $g(U)$ と $g(1-U)$ が正相関になり、分散が増えることすらあります。単調性が効果の前提。
「サンプル数を2倍にしている」ではない： $n$ ペア＝ $2n$ 評価で、プレーン $2n$ と計算量は同じ。フェアな比較で約32倍です。
「 $1-U$ 以外は使えない」ではない：多次元なら各成分を反転、正規乱数なら $Z$ と $-Z$ など、負相関を作れる対なら何でも可。
「不偏性が崩れる」ではない： $1-U$ も一様なので $E[\hat{I}_{\text{anti}}] = I$ 。バイアスは入りません。
「制御変量と同じ」ではない：対照は入力に負相関ペアを仕込む。制御変量（制御変量法）は既知期待値の補助変数で補正する別手法です。

対応シミュレーション参照

本文の $\int_0^1 e^x dx$ 対照変量（default_rng(30)、32倍減）。