制御変量法｜シミュレーション

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

📎 前提：対照変量法（アンチセティック）　|　関連：層化サンプリング

要点（BLUF）

制御変量法：目的量 $f$ と相関し、しかも期待値が既知の補助変数 $c$ を使い、 $f$ の推定を「 $c$ の誤差」で補正します。
最適係数 $b^\* = \text{Cov}(f,c)/\text{Var}(c)$ で補正すると、分散は $(1-\rho^2)$ 倍。相関 $\rho$ が強いほど劇的に減ります。
$\int_0^1 e^x dx$ を制御変量 $c=x$ （期待値0.5、相関0.99）で補正し、分散を約55倍縮小。

1. アイデア：既知の答えで答え合わせ

推定したいのは $\theta = E[f(X)]$ 。ここで、 $f$ と相関し期待値 $E[c] = \mu_c$ が分かっている別の関数 $c(X)$ があるとします。新しい推定量を

\hat{\theta}_{\text{cv}} = \frac{1}{n}\sum_{i=1}^n \Big( f(X_i) - b\,(c(X_i) - \mu_c) \Big)

と作ります。 $E[c(X) - \mu_c] = 0$ なので、どんな係数 $b$ でも $E[\hat{\theta}_{\text{cv}}] = \theta$ （不偏）。狙いは「 $c$ の標本平均が $\mu_c$ からどれだけずれたか」を見て、 $f$ の推定の同じ向きのずれを差し引くことです。 $f$ と $c$ が強く相関していれば、 $c$ のずれは $f$ のずれの良い予測子になり、誤差が消えます。

2. 最適係数と分散

$\hat{\theta}_{\text{cv}}$ の分散は $b$ の2次式で、最小化すると

b^\* = \frac{\text{Cov}(f, c)}{\text{Var}(c)},\qquad \text{Var}(\hat{\theta}_{\text{cv}}) = \text{Var}(f)\,(1 - \rho^2)

ここで $\rho = \text{Corr}(f, c)$ 。分散は $(1-\rho^2)$ 倍になります。 $\rho = 0.99$ なら $1-\rho^2 = 0.0199$ 、すなわち約50倍の改善。 $b^\*$ は $f$ を $c$ に回帰したときの傾きそのもので、制御変量法は「 $f$ から、 $c$ で線形に説明できる部分を取り除く」回帰だと見られます。実務では $b^\*$ を標本から推定します（わずかなバイアスが入るが $n$ 大で無視可能）。

3. 実測

$\int_0^1 e^x dx$ （ $f = e^U$ ）の制御変量に $c = U$ （期待値 $\mu_c = 0.5$ ）を使います。 $e^U$ と $U$ はどちらも単調増加で強く相関します。

import numpy as np

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

def plain(rng):
    u = rng.random(n)
    return np.exp(u).mean()

def control_variate(rng):
    u = rng.random(n)
    f = np.exp(u)
    c = u                                 # 制御変量（期待値0.5）
    b = np.cov(f, c)[0,1] / c.var()       # 最適係数 b* を標本推定
    return (f - b*(c - 0.5)).mean()

p_ests = np.array([plain(rng) for _ in range(reps)])
c_ests = np.array([control_variate(rng) for _ in range(reps)])

print(f"プレーン   平均={p_ests.mean():.5f}  分散={p_ests.var():.2e}")
print(f"制御変量   平均={c_ests.mean():.5f}  分散={c_ests.var():.2e}")
print(f"分散減少率 = {p_ests.var()/c_ests.var():.2f} 倍")

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

出力：

プレーン   平均=1.71830  分散=2.38e-06
制御変量   平均=1.71828  分散=4.34e-08
分散減少率 = 54.82 倍

corr(e^U, U) = 0.9918

出力の意味：両者とも真値 1.71828 を不偏に推定しつつ、分散はプレーン $2.38\times10^{-6}$ に対し制御変量 $4.34\times10^{-8}$ で約55倍減。理論の $1/(1-\rho^2) = 1/(1-0.9918^2) = 61$ にほぼ届きます（ $b^\*$ を標本推定したぶん少し届かない）。相関 0.99 という強い線形関係が効いています。

4. 制御変量の選び方

期待値が解析的に分かること（これが制約）。多項式・線形項・正規分布の積率・粗いモデルの厳密解などが定番。
目的量と強く相関すること。相関がゼロなら効果ゼロ、無理に使うと推定の手間だけ増えます。
複数の制御変量を同時に使う多変量版もあり、回帰の重相関係数で分散が決まります。
金融では「解析解のある単純オプションを制御変量に、複雑オプションを推定」が典型。

数式の直観的意味

制御変量は「カンニングできる似た問題で答え合わせをする」発想です。 $c$ は答え（ $\mu_c$ ）が分かっている問題。同じ乱数でそれを解いてみて、 $\mu_c$ から $\Delta$ ずれていたら、相関する本命 $f$ もおおよそ $b^\*\Delta$ ずれているはず——だから差し引く。 $(1-\rho^2)$ は回帰の「説明できなかった残差の割合」そのもので、相関が強いほど残差が小さく、分散が消えます。対照変量（対照変量法（アンチセティック））が「入力に負相関ペアを仕込む」のに対し、制御変量は「既知の答えを持つ相棒で補正する」——どちらも相関を使う点は共通です。

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

「 $c$ の期待値を知らなくても使える」ではない： $\mu_c$ が既知であることが大前提。未知なら補正項の期待値がゼロにならず、バイアスが入ります。
「 $b$ は1でいい」ではない：最適は $b^\* = \text{Cov}/\text{Var}$ 。 $b=1$ は $f$ と $c$ のスケールが一致する特殊ケースだけ。
「相関が弱くても効く」ではない：効果は $\rho^2$ 依存。相関 0.3 なら分散は 0.91 倍（ほぼ無意味）。強い相関の相棒を探すことが肝心。
「 $b^\*$ を同じ標本で推定すると無効」ではない：わずかなバイアスは入りますが $n$ 大で無視でき、実用上問題ありません。気になれば別標本で $b^\*$ を推定。
「対照変量と併用できない」ではない：併用可能。独立な分散減少効果を重ねられます。

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

本文の制御変量（default_rng(31)、55倍減、 $c=U$ ）。