Day 021 — 確率の基礎②：サイコロとシミュレーション

📚 背景知識（読んでから問題へ）

確率を「実験」で確かめる

確率の理論値（例: サイコロで1が出る確率 = 1/6）は正しいはずだが、実際に10回投げても1が出ない日もある。

これは大数の法則と呼ばれる現象で、「試行回数を増やせば増やすほど、実験結果の確率は理論値に近づく」という法則。

モンテカルロ法

「乱数を大量に発生させて、確率や期待値を数値的に求める」手法。

シミュレーションで難しい確率問題を「力ずく」で解く方法としてKaggleでも使われる。

Python での乱数

📋 コピーimport numpy as np
np.random.seed(42)  # 再現性のために固定

# 0〜1の一様乱数を10個生成
np.random.random(10)

# 1〜6のサイコロ（整数）を1000回
np.random.randint(1, 7, 1000)  # 1以上7未満 = 1〜6

📝 問題

以下の3つをPythonで実装してください。

Q1. サイコロシミュレーション

サイコロを N=100, 1000, 10000, 100000回 ふったとき、各目（1〜6）が出た割合を計算し、理論値（1/6 ≈ 0.1667）との差が N と共にどう変化するかを確認せよ。

期待する出力（例）:

N=100:    [0.13, 0.17, 0.16, 0.18, 0.20, 0.16]  最大誤差: 0.033
N=1000:   [0.168, 0.165, 0.171, ...]              最大誤差: 0.005
N=10000:  ...                                      最大誤差: 0.002
N=100000: ...                                      最大誤差: 0.001

Q2. 誕生日のパラドックス（モンテカルロ法）

N人のグループの中に誕生日が同じ人が2人以上いる確率は、N=23人で50%を超えることが知られている。

モンテカルロシミュレーション（10000回試行）で、N=1〜60 の各Nについてこの確率を求め、確率が0.5を超える最小のNを答えよ。

Q3. 期待値の確認

サイコロを1回ふったとき、出た目の期待値（平均）の理論値は (1+2+3+4+5+6)/6 = 3.5 である。

N=10000回のシミュレーションで、実際の平均値を計算して理論値と比較せよ。

🔍 ヒント（段階的開示）

📋 コピーbirthdays = np.random.randint(1, 366, 23)  # 1〜365日
len(birthdays) != len(set(birthdays))       # 重複あり = True

📋 コピーimport numpy as np
np.random.seed(42)

# Q1
for N in [100, 1000, 10000, 100000]:
    rolls = np.random.randint(1, 7, N)
    counts = np.bincount(rolls, minlength=7)[1:]  # index 0を除く
    probs = counts / N
    max_error = np.max(np.abs(probs - 1/6))
    print(f"N={N:7d}: 最大誤差 {max_error:.4f}")

# Q2
def birthday_simulation(n_people, n_trials=10000):
    count = 0
    for _ in range(n_trials):
        birthdays = np.random.randint(1, 366, n_people)
        if len(set(birthdays)) < n_people:
            count += 1
    return count / n_trials

for n in range(1, 61):
    prob = birthday_simulation(n)
    if prob >= 0.5:
        print(f"N={n}人で確率={prob:.3f}（0.5超え）")
        break

# Q3
rolls = np.random.randint(1, 7, 10000)
print(f"シミュレーション平均: {np.mean(rolls):.4f}（理論値: 3.5000）")

✅ 模範解答

📋 コピーimport numpy as np
np.random.seed(42)

print("=" * 50)
print("Q1. サイコロシミュレーション（大数の法則）")
print("=" * 50)

for N in [100, 1000, 10000, 100000]:
    rolls = np.random.randint(1, 7, N)
    counts = np.bincount(rolls, minlength=7)[1:]  # インデックス0は不要
    probs = counts / N
    max_error = np.max(np.abs(probs - 1/6))
    print(f"N={N:7d}: 各目の確率 = {[f'{p:.3f}' for p in probs]}  最大誤差: {max_error:.4f}")

print()
print("=" * 50)
print("Q2. 誕生日のパラドックス")
print("=" * 50)

def birthday_prob(n_people, n_trials=10000):
    count = sum(
        1 for _ in range(n_trials)
        if len(set(np.random.randint(1, 366, n_people))) < n_people
    )
    return count / n_trials

threshold_n = None
for n in range(1, 61):
    prob = birthday_prob(n)
    if prob >= 0.5 and threshold_n is None:
        threshold_n = n
        print(f"→ N={n}人で確率が50%超え: {prob:.3f}")
    if n in [10, 20, 23, 30, 40, 50, 60]:
        print(f"  N={n:2d}人: 確率 = {birthday_prob(n):.3f}")

print()
print("=" * 50)
print("Q3. 期待値の確認")
print("=" * 50)

rolls = np.random.randint(1, 7, 10000)
simulated_mean = np.mean(rolls)
theoretical_mean = 3.5
print(f"シミュレーション平均: {simulated_mean:.4f}")
print(f"理論値:               {theoretical_mean:.4f}")
print(f"誤差:                 {abs(simulated_mean - theoretical_mean):.4f}")

🪜 Step-by-Step 解説

📐 数学・統計の補足（文系向け）

大数の法則を身近な例で:

コインを10回投げて全部表だったとしても、コインが偏っているわけではない。

100回、1000回投げると表と裏の比率は50:50に近づいていく。

モンテカルロ法の名前の由来:

カジノで有名なモナコの「モンテカルロ」から。ランダムなサンプリングで問題を解くイメージ。

🏆 Kaggleでの実践的な使い方

1. クロスバリデーションの安定性: fold数が多いほど精度評価が安定する（大数の法則）
2. Adversarial Validation: テストデータと学習データの分布確認にも確率的手法を使う
3. アンサンブルの乱数seed: 複数のseedで学習してバギング的に精度を安定させる

⚠️ よくある誤解・ミス

🚀 次のステップ

• 発展: モンテカルロ法で円周率πを近似する（点が円内に入る確率を使う）
• 次回予告: 条件付き確率（「Aが起きた後、Bが起きる確率」）

🎯 自己評価

自分の回答

気づき・メモ