出力解析（過渡・定常・複製）

🎓 レベル：発展　|　重要度：A（必須）

📎 前提：待ち行列のシミュレーション　|　関連：推定量の信頼区間

要点（BLUF）

DES の出力は確率的なので、点推定だけでなく過渡の扱いと信頼区間まで含めて統計的に解析する必要があります。
過渡（warm-up）：空状態から始めると、定常状態に達するまで系が「軽い」ので、定常値を過小評価します。
対策は (1) ウォームアップ区間の除去、(2) 独立な複製（replication）による信頼区間。M/M/1 で立ち上がりと収束を実証します。

1. 問題：空状態スタートの過渡バイアス

待ち行列シミュレーションは普通「行列ゼロ・サーバ空き」から始めます。すると最初の数客はほとんど待たず、系が定常の混雑に達するまで時間がかかる。この過渡期（warm-up period）のデータを混ぜると、定常状態の平均待ち時間を過小評価します。

$i$ 番目に到着した客の待ち時間を、多数の独立な複製で平均すると、立ち上がりの様子が見えます（ $\rho=0.9$ 、定常値 $W_q=9$ ）。

import numpy as np

def mm1_waits(lam, mu, n_customers, seed):
    """単一窓口の各客の待ち時間（イベント駆動）"""
    rng = np.random.default_rng(seed)
    inter = rng.exponential(1/lam, n_customers)
    serv = rng.exponential(1/mu, n_customers)
    arrival = np.cumsum(inter)
    start = np.empty(n_customers); depart = np.empty(n_customers)
    start[0] = arrival[0]; depart[0] = start[0] + serv[0]
    for i in range(1, n_customers):
        start[i] = max(arrival[i], depart[i-1]); depart[i] = start[i] + serv[i]
    return start - arrival

lam, mu = 0.9, 1.0
R, N = 3000, 200
W = np.array([mm1_waits(lam, mu, N, seed=1000+r) for r in range(R)])
avg_by_index = W.mean(axis=0)            # i番目の客の平均待ち（R複製平均）
Wq_th = (lam/mu)/(mu-lam)

print(f"定常値 Wq = {Wq_th:.3f}")
for i in [0, 4, 19, 49, 99, 199]:
    print(f"客 #{i+1:>4}: {R}複製平均の待ち = {avg_by_index[i]:.3f}")

出力：

定常値 Wq = 9.000
客 #   1: 3000複製平均の待ち = 0.000
客 #   5: 3000複製平均の待ち = 1.363
客 #  20: 3000複製平均の待ち = 3.315
客 #  50: 3000複製平均の待ち = 4.977
客 # 100: 3000複製平均の待ち = 6.352
客 # 200: 3000複製平均の待ち = 7.368

出力の意味：1番目の客は必ず待ち0、そこから $1.36 \to 3.32 \to 4.98 \to 6.35 \to 7.37$ と、定常値9へゆっくり立ち上がっていきます（ $\rho=0.9$ は混雑が強く収束が遅い）。200客目でもまだ7.4で定常に達していません。この立ち上がり区間（warm-up）を平均に含めると、定常 $W_q$ を過小評価する——これが過渡バイアスです。

2. 対策1：ウォームアップ除去と長時間化

短い run の平均が定常値をどれだけ過小評価するか、run 長を変えて確かめます。

import numpy as np

def mm1_waits(lam, mu, n_customers, seed):
    rng = np.random.default_rng(seed)
    inter = rng.exponential(1/lam, n_customers); serv = rng.exponential(1/mu, n_customers)
    arrival = np.cumsum(inter)
    start = np.empty(n_customers); depart = np.empty(n_customers)
    start[0]=arrival[0]; depart[0]=start[0]+serv[0]
    for i in range(1, n_customers):
        start[i]=max(arrival[i], depart[i-1]); depart[i]=start[i]+serv[i]
    return start - arrival

lam, mu = 0.9, 1.0
Wq_th = (lam/mu)/(mu-lam)
for n in [50, 200, 1000, 10000]:
    means = np.array([mm1_waits(lam, mu, n, seed=2000+r).mean() for r in range(500)])
    print(f"n={n:>6}: 平均待ち(500複製平均) = {means.mean():.3f}  (定常 {Wq_th:.3f})")

出力：

n=    50: 平均待ち(500複製平均) = 3.526  (定常 9.000)
n=   200: 平均待ち(500複製平均) = 5.823  (定常 9.000)
n=  1000: 平均待ち(500複製平均) = 8.021  (定常 9.000)
n= 10000: 平均待ち(500複製平均) = 8.668  (定常 9.000)

出力の意味：run が短いほど過渡の比重が大きく、 $n=50$ では 3.5（定常の4割）と激しく過小評価。 $n$ を増やすと過渡が薄まり 8.67 まで改善しますが、まだ9に届かない。対策は**(a) 先頭の warm-up 区間を捨ててから集計する**、(b) run を十分長くして過渡を相対的に無視できるようにするの2本立てです。

3. 対策2：複製による信頼区間

run を1本だけ走らせても、その平均がどれだけ信頼できるか分かりません。独立シードで複数回（複製）走らせ、複製間の平均と標準誤差から信頼区間を作ります。

import numpy as np

def mm1_mean_wait(lam, mu, n, seed):
    rng = np.random.default_rng(seed)
    inter = rng.exponential(1/lam, n); serv = rng.exponential(1/mu, n)
    arrival = np.cumsum(inter)
    start = np.empty(n); depart = np.empty(n)
    start[0]=arrival[0]; depart[0]=start[0]+serv[0]
    for i in range(1, n):
        start[i]=max(arrival[i], depart[i-1]); depart[i]=start[i]+serv[i]
    return (start - arrival).mean()

lam, mu = 0.9, 1.0
Wq_th = (lam/mu)/(mu-lam)
means = np.array([mm1_mean_wait(lam, mu, 5000, seed=3000+r) for r in range(40)])
se = means.std(ddof=1) / np.sqrt(len(means))
print(f"40複製(各n=5000): {means.mean():.3f} +/- {1.96*se:.3f}  (定常 {Wq_th:.3f})")

出力：

40複製(各n=5000): 9.042 +/- 0.949  (定常 9.000)

出力の意味：40複製それぞれが独立なので、複製平均の標準偏差から標準誤差が計算でき、95%信頼区間 $9.042 \pm 0.949 = [8.09, 9.99]$ が定常値9.0を含みます。各 run を長く（ $n=5000$ ）取って過渡を薄め、複製で誤差を評価する——これが DES 出力解析の王道です。複製法のほかに、1本の長い run を区間に分けるバッチ平均法もあります。

数式の直観的意味

過渡バイアスは「系がまだ温まっていない」状態の混入です。M/M/1 の待ち時間は前の客の状態を引き継ぐ（自己相関が強い）ので、空スタートの影響が長く尾を引く—— $\rho$ が1に近いほど緩和時間が長く、収束が遅い。だから単純な long-run 平均は過小評価する。複製法が機能するのは、各 run を独立シードで回すと run 間が独立になり、中心極限定理が複製平均に効いて、 $\sigma_{\text{run}}/\sqrt{R}$ で標準誤差が測れるから。1本の run 内のデータは自己相関しているので素朴な標準誤差が使えない——ここが「複製 or バッチ」が必要な理由です。

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

「1回走らせた平均が答え」ではない：確率的なので必ず誤差。複製で信頼区間を付けます。
「run 内の各客で標準誤差を計算」してはいけない：客の待ち時間は強く自己相関するので、独立を仮定した標準誤差は過小。複製かバッチ平均を使う。
「warm-up はいつもゼロでよい」ではない：混雑系（ $\rho$ 大）ほど過渡が長い。除去しないと過小評価。
「長く回せば warm-up は無視できる」は半分正しい：相対的に薄まりますが、 $\rho \to 1$ では収束が極端に遅く、十分長くする必要があります。
「定常分布が常に存在する」ではない： $\rho \ge 1$ では定常がなく、過渡という概念自体が成立しません（行列は無限に伸びる）。

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

本文の過渡曲線（default_rng(1000+r)）・short-run バイアス（default_rng(2000+r)）・複製信頼区間（default_rng(3000+r)）。