共通乱数法（CRN）｜シミュレーション

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

📎 前提：対照変量法（アンチセティック）　|　関連：シミュレーション最適化

要点（BLUF）

共通乱数法（CRN, Common Random Numbers）：複数のシステム構成を比較するとき、各構成のシミュレーションに同じ乱数列を使うことで、構成間の差の分散を下げます。
「比較したいのは差なのに、各構成のノイズが独立だと差が見えにくい」——同じ乱数を共有すれば、共通由来のばらつきが相殺されます。
2つの M/M/1 待ち行列（ $\mu=1.0$ と $1.05$ ）の平均待ち時間差の分散が、独立乱数より約10倍小さくなることを実証します。

1. 問題：差を見たいのにノイズが邪魔

シミュレーションでよくあるのは「構成Aと構成B、どちらが良いか」という比較です（窓口1個 vs 2個、発注点 s=20 vs 30…）。知りたいのは差 $\Delta = E[Y_A] - E[Y_B]$ 。これを推定 $\hat{\Delta} = \bar{Y}_A - \bar{Y}_B$ するとき、分散は

\text{Var}(\hat{\Delta}) = \text{Var}(\bar{Y}_A) + \text{Var}(\bar{Y}_B) - 2\,\text{Cov}(\bar{Y}_A, \bar{Y}_B)

各構成を独立な乱数でシミュすると共分散はゼロ、分散は単純な和。これだと、差 $\Delta$ が小さいとき、各構成のノイズに埋もれて差が判別できません。

2. アイデア：同じ乱数で正の相関を作る

CRN は、構成AとBを同じ乱数列で走らせます。すると両者が共通の「運」を共有するので、 $\hat{Y}_A$ と $\hat{Y}_B$ が正に相関（ $\text{Cov} > 0$ ）し、上式の $-2\text{Cov}$ が効いて分散が下がる。

直観は明快——「同じ客の流れ・同じサービス時間で2つの窓口数を比べれば、条件が揃うので差が純粋に構成の違いだけを反映する」。乱数が違うと、Aが空いていたのは構成が良いのか、たまたま客が少なかったのか区別できません。対照変量（対照変量法（アンチセティック））が負の相関を作って和の分散を下げたのに対し、CRN は正の相関を作って差の分散を下げます。

3. 実装：2つの M/M/1 を比較する

$\lambda=0.8$ で、サービス率 $\mu=1.0$ と $\mu=1.05$ の M/M/1 の平均待ち時間差を推定します。CRN では、到着間隔と「基準サービス時間（レート1の指数）」を両系で共有し、 $\mu$ でスケールするだけにします。

import numpy as np

def mean_wait(inter, base_serv, mu):
    """共有した到着・基準サービスから平均待ち時間を計算"""
    serv = base_serv / mu                       # 基準Exp(1)をExp(mu)にスケール
    n = len(inter)
    arr = np.cumsum(inter)
    start = np.empty(n); depart = np.empty(n)
    start[0] = arr[0]; depart[0] = start[0] + serv[0]
    for i in range(1, n):
        start[i] = max(arr[i], depart[i-1]); depart[i] = start[i] + serv[i]
    return (start - arr).mean()

lam = 0.8; n = 2000; reps = 2000

# 共通乱数法：1複製内で両系が同じ乱数を共有
rng = np.random.default_rng(70)
diffs_crn = np.empty(reps)
for r in range(reps):
    inter = rng.exponential(1/lam, n)
    base_serv = rng.exponential(1.0, n)
    diffs_crn[r] = mean_wait(inter, base_serv, 1.0) - mean_wait(inter, base_serv, 1.05)

# 独立乱数：各系で別々の乱数
rng = np.random.default_rng(71)
diffs_ind = np.empty(reps)
for r in range(reps):
    iA = rng.exponential(1/lam, n); sA = rng.exponential(1.0, n)
    iB = rng.exponential(1/lam, n); sB = rng.exponential(1.0, n)
    diffs_ind[r] = mean_wait(iA, sA, 1.0) - mean_wait(iB, sB, 1.05)

print(f"平均待ち時間の差（理論 4.000 - 3.048 = 0.952）")
print(f"  CRN 推定差={diffs_crn.mean():.4f}  分散={diffs_crn.var(ddof=1):.5f}")
print(f"  独立 推定差={diffs_ind.mean():.4f}  分散={diffs_ind.var(ddof=1):.5f}")
print(f"  分散減少率 = {diffs_ind.var(ddof=1)/diffs_crn.var(ddof=1):.1f} 倍")

出力：

平均待ち時間の差（理論 4.000 - 3.048 = 0.952）
  CRN 推定差=0.9296  分散=0.12261
  独立 推定差=0.9457  分散=1.24445
  分散減少率 = 10.1 倍

出力の意味：理論的な待ち時間差は $W_q(\mu{=}1.0) - W_q(\mu{=}1.05) = 4.000 - 3.048 = 0.952$ 。CRN・独立どちらの推定差もこれに近い（不偏性は保たれる）が、差の分散は CRN 0.123 に対し独立 1.244 で約10倍の改善。つまり同じ精度で2構成の優劣を判定するのに、CRN は独立法の約1/10のサンプルで済む。サービス率をわずか0.05上げる効果のような小さな差を検出したいときほど、CRN は強力です。

4. 実装の作法と注意

乱数の同期（synchronization）が命：両構成で「同じ乱数が同じ目的に使われる」ようにそろえる。本コードは到着・基準サービスを共有し $\mu$ でスケールするだけにして同期を保ちました。
乱数ストリームを用途別に分ける：到着用・サービス用などストリームを分けておくと、構成変更で乱数の使われ方がずれにくい（同期崩れの防止）。
正の相関が前提：構成変更が出力に単調・滑らかに効くと正相関が強い。乱数が複雑に絡む系では相関が弱まり効果が減ります（まれに逆効果）。
比較・最適化の土台：CRN はシミュレーション最適化や複数案のランキング・選択で標準的に使われます。

数式の直観的意味

CRN は「条件をそろえて、純粋に構成の違いだけを見る」操作です。差の分散 $\text{Var}(\hat{\Delta}) = \text{Var}_A + \text{Var}_B - 2\text{Cov}$ で、 $\text{Cov}$ を大きくすれば差の分散が縮む。同じ乱数を使うと、両系が「同じ客の波・同じ運」を浴びるので出力が一緒に上下し（正相関）、その共通成分が引き算で消える。残るのは構成の違いに起因する差だけ——だからノイズが小さい。対照変量（対照変量法（アンチセティック））の鏡像で、あちらは「和の分散を負相関で消す」、こちらは「差の分散を正相関で消す」。どちらも相関を設計してばらつきを相殺するという分散減少の共通原理の現れです。「Aが良かったのは運か実力か」を切り分けるために条件をそろえる、という実験計画の発想とも一致します。

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

「同じ乱数を使うと結果が偏る」ではない：各構成の推定は不偏のまま。変わるのは差の分散だけ。
「同期は気にしなくていい」ではない：乱数の使われ方が両構成でずれると相関が崩れ、効果が消える（最悪逆効果）。同期が CRN の生命線。
「絶対性能の推定にも効く」ではない：CRN が下げるのは差の分散。各構成単体の性能推定の分散は下げません。
「いつでも正相関」ではない：構成変更が出力に非単調・不連続に効くと相関が弱い・負になることも。事前に相関を確認。
「対照変量と同じ」ではない：対照は1つの推定で負相関ペア、CRN は2つの系で正相関を共有。目的（和 vs 差）が逆です。

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

本文の2つの M/M/1 比較（CRN default_rng(70)／独立 default_rng(71)、差の分散10倍減）。