← 目次 ｜ ← 前: 15 ｜ 次: 17 →

16. 離散事象シミュレーション

テーマ: 窓口待ち行列 M/M/1 のシミュレーションと理論値の照合

学習点: random モジュール, イベント駆動ループ, 統計量の収集, 理論式との突き合わせ（検証の作法）

依存: 標準ライブラリのみ ／ 難易度: 中級

実行方法

uv run 16_queue_mm1.py

スクリプト冒頭の PEP 723 メタデータ(# /// script)により、必要なライブラリは uv が自動的に仮想環境へ導入します。事前の pip install は不要です。

解説

何をするプログラムか

銀行の窓口やコールセンターのように、客がランダムに到着し 1 つの窓口で順番に処理される状況を「待ち行列 M/M/1 モデル」として 20 万人分シミュレートし、平均待ち時間 Wq・平均系内時間 W・平均系内客数 L を理論式と突き合わせるプログラムです。

経営上の関心は「窓口の稼働率 ρ を上げると待ち時間はどうなるか」です。稼働率 95% は一見効率的ですが、待ち時間は稼働率 50% の約 19 倍に膨らみます。シミュレーションと理論の両面からこの非線形性を確認します。

コードの読みどころ

• simulate_mm1() は rng.expovariate(lam) で到着間隔を、rng.expovariate(mu) でサービス時間を指数分布から生成します。random.Random(seed) で乱数生成器を固定しているため、何度実行しても同じ結果が再現されます。
• イベント駆動ループの核心は start = max(t_arrive, t_free) の 1 行です。到着時刻と「窓口が空く時刻」の遅い方がサービス開始時刻になり、wait = start - t_arrive が待ち時間になります。時刻を 1 秒ずつ進めるのではなく、イベント（到着）単位で時計を飛ばすのが離散事象シミュレーションの発想です。
• 統計量は waits と systems のリストに収集し、statistics.mean() で平均を取ります。系内客数 L は直接数えず、リトルの法則 l_sim = lam * w_sim で換算しているのも読みどころです。
• 待ち時間の分位点は sorted(waits) の上で waits_sorted[int(n*q)] と添字計算で取り出しており、平均では見えない「上位 1% の客の待ち時間」を可視化します。

理論的背景

到着がポアソン過程（率 λ）、サービスが指数分布（率 μ）、窓口 1 つのモデルを M/M/1 と呼びます。利用率 ρ=λ/μ<1 のとき定常状態が存在し、平均系内客数 L = ρ/(1-ρ)、平均系内時間 W = 1/(μ-λ)、平均待ち時間 Wq = ρ/(μ-λ) が成り立ちます。また、任意の定常な待ち行列で成り立つリトルの法則 L = λW を使えば、時間の統計から客数の統計へ変換できます。分母の (1-ρ) が待ち時間の発散を生む点が要点です。

実行結果の見方

最初の表では、ρ=0.3〜0.9 まで実測と理論がほぼ一致していることを確認してください。ρ=0.95 では実測 21.0 分に対し理論 19.0 分とぶれが大きくなりますが、これは混雑時ほど自己相関が強く収束が遅いためで、シミュレーションを理論で検証する「検証の作法」の教材になっています。

「混雑の非線形性」の棒グラフでは、稼働率 95%→98% で待ち時間が 19 分→49 分と約 2.6 倍になる点、最後の分布では平均 3.97 分に対し 99% 分位点が 22.45 分・最大 58.6 分と、平均だけでは裾のリスクを見落とすことが読み取れます。

ソースコード

# /// script
# requires-python = ">=3.11"
# dependencies = []
# ///
"""16: 離散事象シミュレーション -------------------------------------------
テーマ: 窓口待ち行列 M/M/1 のシミュレーションと理論値の照合
学習点: random モジュール, イベント駆動ループ, 統計量の収集,
        理論式との突き合わせ（検証の作法）
根拠: 到着がポアソン過程(率λ)、サービスが指数分布(率μ)、窓口1つのとき
      利用率 ρ=λ/μ<1 で定常状態が存在し、
      平均系内客数 L = ρ/(1-ρ)、平均系内時間 W = 1/(μ-λ)、
      平均待ち時間 Wq = ρ/(μ-λ)。またリトルの法則 L = λW が成り立つ。
"""
import random
import statistics as st


def simulate_mm1(lam: float, mu: float, n: int = 200_000, seed: int = 1):
    """n人の客を処理する M/M/1 をシミュレートし待ち時間の統計を返す。"""
    rng = random.Random(seed)
    t_arrive = 0.0
    t_free = 0.0          # 窓口が空く時刻
    waits, systems = [], []
    for _ in range(n):
        t_arrive += rng.expovariate(lam)          # 到着間隔 ~ Exp(λ)
        start = max(t_arrive, t_free)
        wait = start - t_arrive
        service = rng.expovariate(mu)             # サービス時間 ~ Exp(μ)
        t_free = start + service
        waits.append(wait)
        systems.append(wait + service)
    busy_time = sum(s - w for s, w in zip(systems, waits))
    return waits, systems, busy_time / t_free


def main() -> None:
    mu = 1.0        # 平均サービス時間 1 分
    print(f"{'λ':>6}{'ρ':>7}{'Wq実測':>10}{'Wq理論':>10}"
          f"{'W実測':>9}{'W理論':>9}{'L実測':>9}{'L理論':>9}")
    print("-" * 70)
    for lam in [0.3, 0.5, 0.7, 0.8, 0.9, 0.95]:
        waits, systems, util = simulate_mm1(lam, mu)
        wq_sim, w_sim = st.mean(waits), st.mean(systems)
        wq_th, w_th = lam / (mu * (mu - lam)), 1 / (mu - lam)
        l_sim = lam * w_sim                       # リトルの法則で換算
        l_th = lam / (mu - lam)
        print(f"{lam:>6.2f}{lam/mu:>7.2f}{wq_sim:>10.3f}{wq_th:>10.3f}"
              f"{w_sim:>9.3f}{w_th:>9.3f}{l_sim:>9.3f}{l_th:>9.3f}")

    print("  ※ ρ が 1 に近いほど自己相関が強く収束が遅いため、"
          "同じ客数でも実測は理論値からぶれる。")

    print("\n[混雑の非線形性] ρ が 1 に近づくと待ち時間は発散する")
    for rho in [0.5, 0.8, 0.9, 0.95, 0.98, 0.99]:
        print(f"  ρ={rho:>5.2f}: Wq = {rho/(1-rho):>8.2f} 分 "
              f"{'█' * min(int(rho/(1-rho)*1.2), 60)}")
    print("  ※ 稼働率を95%→98%に上げると待ち時間は約2.6倍になる。")
    print("    「窓口を遊ばせない」運用が顧客待ち時間を爆発させる典型例。")

    print("\n[待ち時間分布] λ=0.8, μ=1.0")
    waits, _, util = simulate_mm1(0.8, 1.0, 200_000)
    waits_sorted = sorted(waits)
    n = len(waits_sorted)
    print(f"  窓口稼働率(実測) = {util:.4f} (理論 0.8000)")
    print(f"  待たずに済んだ客の割合 = {sum(w < 1e-12 for w in waits)/n:.3f} "
          "(理論 1-ρ = 0.200)")
    for q in [0.50, 0.75, 0.90, 0.95, 0.99]:
        print(f"  {q:>5.0%} 分位点: {waits_sorted[int(n*q)]:>7.3f} 分")
    print(f"  平均 {st.mean(waits):.3f} 分 / 最大 {max(waits):.2f} 分")
    print("  ※ 平均だけを見ると、上位1%の顧客が被る待ち時間を見落とす。")


if __name__ == "__main__":
    main()

実行結果

     λ      ρ      Wq実測      Wq理論      W実測      W理論      L実測      L理論
----------------------------------------------------------------------
  0.30   0.30     0.430     0.429    1.431    1.429    0.429    0.429
  0.50   0.50     0.999     1.000    2.000    2.000    1.000    1.000
  0.70   0.70     2.314     2.333    3.315    3.333    2.320    2.333
  0.80   0.80     3.973     4.000    4.974    5.000    3.979    4.000
  0.90   0.90     8.862     9.000    9.863   10.000    8.876    9.000
  0.95   0.95    21.035    19.000   22.036   20.000   20.934   19.000
  ※ ρ が 1 に近いほど自己相関が強く収束が遅いため、同じ客数でも実測は理論値からぶれる。

[混雑の非線形性] ρ が 1 に近づくと待ち時間は発散する
  ρ= 0.50: Wq =     1.00 分 █
  ρ= 0.80: Wq =     4.00 分 ████
  ρ= 0.90: Wq =     9.00 分 ██████████
  ρ= 0.95: Wq =    19.00 分 ██████████████████████
  ρ= 0.98: Wq =    49.00 分 ██████████████████████████████████████████████████████████
  ρ= 0.99: Wq =    99.00 分 ████████████████████████████████████████████████████████████
  ※ 稼働率を95%→98%に上げると待ち時間は約2.6倍になる。
    「窓口を遊ばせない」運用が顧客待ち時間を爆発させる典型例。

[待ち時間分布] λ=0.8, μ=1.0
  窓口稼働率(実測) = 0.7993 (理論 0.8000)
  待たずに済んだ客の割合 = 0.201 (理論 1-ρ = 0.200)
    50% 分位点:   2.303 分
    75% 分位点:   5.741 分
    90% 分位点:  10.249 分
    95% 分位点:  13.782 分
    99% 分位点:  22.451 分
  平均 3.973 分 / 最大 58.63 分
  ※ 平均だけを見ると、上位1%の顧客が被る待ち時間を見落とす。

---

← 目次 ｜ ← 前: 15 ｜ 次: 17 →