# 2018 Q2(1)-2  2024.10.5

import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
N_values = [50, 70, 100]  # 総球数 N のバリエーション
M_steps = 10  # 赤球のステップ数を指定
n_trials = 10000  # シミュレーションの試行回数

# グラフ描画の準備
plt.figure(figsize=(10, 6))

# 各 N に対してシミュレーションと理論値を計算
for N in N_values:
    M_values = np.arange(2, N + 1, M_steps)  # 赤球数をステップで変化 (M >= 2 に設定)
    
    # 理論値の計算
    P_Xi_Xj_theory = (M_values / N) * ((M_values - 1) / (N - 1))
    
    # シミュレーションによる確率計算
    P_Xi_Xj_simulation = []
    
    for M in M_values:
        count_Xi_Xj_1 = 0
        for _ in range(n_trials):
            # 箱の中の球を準備（赤球 M 個、青球 N-M 個）
            balls = [1] * M + [0] * (N - M)
            
            # 非復元抽出で2つの球を無作為に抽出
            drawn_balls = np.random.choice(balls, size=2, replace=False)
            
            # 2回連続で赤球が引かれたか確認
            if drawn_balls[0] == 1 and drawn_balls[1] == 1:
                count_Xi_Xj_1 += 1
        
        # シミュレーション結果の確率
        P_Xi_Xj_simulation.append(count_Xi_Xj_1 / n_trials)
    
    # グラフに追加描画
    plt.plot(M_values, P_Xi_Xj_theory, label=f'理論値 (N={N})', linestyle='--')
    plt.plot(M_values, P_Xi_Xj_simulation, label=f'シミュレーション (N={N})', marker='o', linestyle='None')

# グラフのラベルとタイトル設定
plt.xlabel('赤球の数 (M)')
plt.ylabel('P(X_i = 1, X_j = 1) の確率')
plt.title('異なる N における P(X_i = 1, X_j = 1) の理論値とシミュレーション結果の比較')
plt.legend()
plt.grid(True)
plt.show()