# 2018 Q2(5)-2  2024.10.10

import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
N_true = 1000  # 真の青球の数 N (これを推定する)
K_values = [100, 300, 500]  # 赤球の数 K のバリエーション
n_values = [100, 200, 300]  # 抽出する回数 n のバリエーション
n_trials = 10000  # シミュレーションの試行回数

# グラフレイアウト設定
fig, axes = plt.subplots(len(K_values), len(n_values), figsize=(15, 12))
axes = axes.flatten()

plot_idx = 0

# 各 K, n の組み合わせに対してシミュレーションを実行
for K in K_values:
    for n in n_values:
        N_hat_values = []

        for _ in range(n_trials):
            # 青球 N 個 + 赤球 K 個のボールを作成
            balls = [0] * N_true + [1] * K  # 0が青球、1が赤球
            drawn_balls = np.random.choice(balls, size=n, replace=False)  # 非復元抽出
            X = np.sum(drawn_balls)  # 抽出された赤球の個数

            # X を使って N の推定量 N_hat を計算
            if X > 0:  # Xが0でない場合
                N_hat = K * (n - X) / X
                N_hat_values.append(N_hat)

        # シミュレーション結果から N_hat の分散を計算
        N_hat_var = np.var(N_hat_values)
        
        # 理論値 ε の計算
        epsilon_theory = np.sqrt((N_true + K) * (N_true + K - n) / (n * N_true * K))
        
        # シミュレーションから求めた ε の計算
        epsilon_simulation = np.sqrt(N_hat_var) / N_true
       
        # 推定量 N_hat の分布をヒストグラムで表示
        ax = axes[plot_idx]
        ax.hist(N_hat_values, bins=100, density=True, alpha=0.7, color='b', edgecolor='black')
        ax.axvline(x=N_true, color='r', linestyle='--')  # 真の N を破線で表示
        ax.set_title(f'K={K}, n={n}')
        ax.set_xlabel('推定量 N_hat')
        ax.set_ylabel('確率密度')

        # εの値と真のNをまとめて表示
        ax.text(0.95, 0.95, f'真の N = {N_true}\nε (理論) = {epsilon_theory:.4f}\nε (シミュ) = {epsilon_simulation:.4f}',
                transform=ax.transAxes, verticalalignment='top', horizontalalignment='right',
                bbox=dict(facecolor='white', alpha=0.7), fontsize=10)
        
        ax.grid(True)

        plot_idx += 1

# グラフ全体の調整
plt.tight_layout()
plt.show()

# 2018 Q2(5)-1  2024.10.9

import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
N_true = 100  # 真の青球の数 N (これを推定する)
K = 30  # 赤球の数 K (既知)
n = 20  # 抽出する回数
n_trials = 10000  # シミュレーションの試行回数

# 推定量 Nハットのシミュレーション
N_hat_values = []

for _ in range(n_trials):
    # 青球 N 個 + 赤球 K 個のボールを作成
    balls = [0] * N_true + [1] * K  # 0が青球、1が赤球
    drawn_balls = np.random.choice(balls, size=n, replace=False)  # 非復元抽出
    X = np.sum(drawn_balls)  # 抽出された赤球の個数
    
    # X を使って N の推定量 N_hat を計算
    if X > 0:  # Xが0でない場合
        N_hat = K * (n - X) / X
        N_hat_values.append(N_hat)

# シミュレーション結果の分析
N_hat_mean = np.mean(N_hat_values)
N_hat_var = np.var(N_hat_values)

# 結果表示
print(f"推定量 N_hat の平均: {N_hat_mean}")
print(f"推定量 N_hat の分散: {N_hat_var}")

# ヒストグラム表示
plt.figure(figsize=(10, 6))
plt.hist(N_hat_values, bins=50, density=True, alpha=0.7, color='b', edgecolor='black')
plt.axvline(x=N_true, color='r', linestyle='--', label=f'真の N = {N_true}')
plt.title(f'推定量 N_hat の分布 (N={N_true}, K={K}, n={n})')
plt.xlabel('推定量 N_hat')
plt.ylabel('確率密度')
plt.legend()
plt.grid(True)
plt.show()

推定量 N_hat の平均: 126.09541441704363
推定量 N_hat の分散: 8070.195087794488

# 2018 Q2(4)  2024.10.8

import numpy as np
import matplotlib.pyplot as plt

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

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

for N in N_values:
    M_values = np.arange(M_steps, N + 1, M_steps)  # 赤球数をステップで変化
    
    # 理論値の計算
    E_X_theory = n * M_values / N
    V_X_theory = n * M_values / N * (N - M_values) / N * (N - n) / (N - 1)
    
    # シミュレーションによる期待値と分散の計算
    E_X_simulation = []
    V_X_simulation = []
    
    for M in M_values:
        X_values = []
        for _ in range(n_trials):
            balls = [1] * M + [0] * (N - M)
            drawn_balls = np.random.choice(balls, size=n, replace=False)
            X_values.append(np.sum(drawn_balls))  # 赤球の数
        
        E_X_simulation.append(np.mean(X_values))
        V_X_simulation.append(np.var(X_values))
    
    # 期待値のグラフに追加描画
    plt.subplot(2, 1, 1)  # 期待値のグラフ (上側)
    plt.plot(M_values, E_X_theory, label=f'理論値 (N={N})', linestyle='--')
    plt.plot(M_values, E_X_simulation, label=f'シミュレーション (N={N})', marker='o', linestyle='None')
    
    # 分散のグラフに追加描画
    plt.subplot(2, 1, 2)  # 分散のグラフ (下側)
    plt.plot(M_values, V_X_theory, label=f'理論値 (N={N})', linestyle='--')
    plt.plot(M_values, V_X_simulation, label=f'シミュレーション (N={N})', marker='o', linestyle='None')

# グラフの設定
plt.subplot(2, 1, 1)
plt.xlabel('赤球の数 (M)')
plt.ylabel('期待値 E[X]')
plt.title('異なる N における期待値 E[X] の理論値とシミュレーション結果の比較')
plt.legend()
plt.grid(True)

plt.subplot(2, 1, 2)
plt.xlabel('赤球の数 (M)')
plt.ylabel('分散 V[X]')
plt.title('異なる N における分散 V[X] の理論値とシミュレーション結果の比較')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

# 2018 Q2(3)  2024.10.7

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import hypergeom

# パラメータの設定
N_values = [50, 70, 100]  # 総球数 N のバリエーション (列ごとに固定)
M_values = [10, 20, 30, 40, 50]  # 赤球の数 M のバリエーション (行ごとに変化)
n = 10  # 抽出する回数
n_trials = 10000  # シミュレーションの試行回数

# グラフのレイアウト設定 (5行x3列)
fig, axes = plt.subplots(5, 3, figsize=(15, 20))  # 5行3列のサブプロット

# 各 M に対してシミュレーションと理論値を計算
for row, M in enumerate(M_values):  # 行が M を表す
    for col, N in enumerate(N_values):  # 列が N を表す
        # 理論値の計算
        x_values = np.arange(0, n + 1)  # 可能な赤球の数 x の範囲
        rv = hypergeom(N, M, n)  # 超幾何分布のインスタンス
        P_X_theory = rv.pmf(x_values)  # 理論値（超幾何分布）
        
        # シミュレーションによる確率計算
        drawn_red_balls = []
        for _ in range(n_trials):
            balls = [1] * M + [0] * (N - M)  # 赤球 M 個、青球 N-M 個のリスト
            drawn_balls = np.random.choice(balls, size=n, replace=False)  # n 回非復元抽出
            x = np.sum(drawn_balls)  # 赤球の個数を数える
            drawn_red_balls.append(x)
        
        # 各サブプロットにグラフを描画
        ax = axes[row, col]
        ax.plot(x_values, P_X_theory, 'bo-', label=f'理論値 (N={N}, M={M})')
        ax.hist(drawn_red_balls, bins=np.arange(-0.5, n + 1.5, 1), density=True, alpha=0.6, color='r', label='シミュレーション結果')
        ax.set_xlabel('赤球の個数 (x)')
        ax.set_ylabel('確率 P(X = x)')
        ax.set_title(f'N={N}, M={M}, n={n}')
        ax.grid(True)

# 全体のレイアウト調整
plt.tight_layout()
plt.show()

# 2018 Q2(2)  2024.10.6

import numpy as np
import matplotlib.pyplot as plt

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

# 期待値 E[X_i] のグラフ描画
plt.figure(figsize=(10, 6))
for N in N_values:
    M_values = np.arange(2, N + 1, M_steps)  # 赤球数をステップで変化
    
    # 理論値の計算
    E_Xi_theory = M_values / N
    
    # シミュレーションによる期待値の計算
    E_Xi_simulation = []
    
    for M in M_values:
        count_Xi_1 = 0
        for _ in range(n_trials):
            balls = [1] * M + [0] * (N - M)
            drawn_ball = np.random.choice(balls, size=1, replace=False)
            count_Xi_1 += drawn_ball[0]
        
        E_Xi_simulation.append(count_Xi_1 / n_trials)
    
    # グラフに追加描画
    plt.plot(M_values, E_Xi_theory, label=f'理論値 (N={N})', linestyle='--')
    plt.plot(M_values, E_Xi_simulation, label=f'シミュレーション (N={N})', marker='o', linestyle='None')

plt.xlabel('赤球の数 (M)')
plt.ylabel('期待値 E[X_i]')
plt.title('異なる N における 期待値 E[X_i] の理論値とシミュレーション結果の比較')
plt.legend()
plt.grid(True)
plt.show()

# 2018 Q2(2)  2024.10.6

import numpy as np
import matplotlib.pyplot as plt

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

# 分散 V[X_i] のグラフ描画
plt.figure(figsize=(10, 6))
for N in N_values:
    M_values = np.arange(2, N + 1, M_steps)
    
    # 理論値の計算
    V_Xi_theory = M_values * (N - M_values) / N**2
    
    # シミュレーションによる分散の計算
    V_Xi_simulation = []
    
    for M in M_values:
        Xi_values = []
        for _ in range(n_trials):
            balls = [1] * M + [0] * (N - M)
            drawn_ball = np.random.choice(balls, size=1, replace=False)
            Xi_values.append(drawn_ball[0])
        
        V_Xi_simulation.append(np.var(Xi_values))
    
    # グラフに追加描画
    plt.plot(M_values, V_Xi_theory, label=f'理論値 (N={N})', linestyle='--')
    plt.plot(M_values, V_Xi_simulation, label=f'シミュレーション (N={N})', marker='o', linestyle='None')

plt.xlabel('赤球の数 (M)')
plt.ylabel('分散 V[X_i]')
plt.title('異なる N における 分散 V[X_i] の理論値とシミュレーション結果の比較')
plt.legend()
plt.grid(True)
plt.show()

# 2018 Q2(2)  2024.10.6

import numpy as np
import matplotlib.pyplot as plt

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

# 共分散 Cov[X_i, X_j] のグラフ描画
plt.figure(figsize=(10, 6))
for N in N_values:
    M_values = np.arange(2, N + 1, M_steps)
    
    # 理論値の計算
    Cov_Xi_Xj_theory = -(M_values * (N - M_values)) / (N**2 * (N - 1))
    
    # シミュレーションによる共分散の計算
    Cov_Xi_Xj_simulation = []
    
    for M in M_values:
        Xi_values = []
        Xj_values = []
        for _ in range(n_trials):
            balls = [1] * M + [0] * (N - M)
            drawn_balls = np.random.choice(balls, size=2, replace=False)
            Xi_values.append(drawn_balls[0])
            Xj_values.append(drawn_balls[1])
        
        Cov_Xi_Xj_simulation.append(np.cov(Xi_values, Xj_values)[0, 1])
    
    # グラフに追加描画
    plt.plot(M_values, Cov_Xi_Xj_theory, label=f'理論値 (N={N})', linestyle='--')
    plt.plot(M_values, Cov_Xi_Xj_simulation, label=f'シミュレーション (N={N})', marker='o', linestyle='None')

plt.xlabel('赤球の数 (M)')
plt.ylabel('共分散 Cov[X_i, X_j]')
plt.title('異なる N における 共分散 Cov[X_i, X_j] の理論値とシミュレーション結果の比較')
plt.legend()
plt.grid(True)
plt.show()

# 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()

# 2018 Q2(1)-1  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(1, N + 1, M_steps)  # 赤球数をステップで変化
    P_Xi_1_theory = M_values / N  # 理論値
    
    # シミュレーションによる確率計算
    P_Xi_1_simulation = []
    
    for M in M_values:
        count_Xi_1 = 0
        for _ in range(n_trials):
            # 箱の中の球を準備（赤球 M 個、青球 N-M 個）
            balls = [1] * M + [0] * (N - M)
            
            # 1回の無作為抽出で赤球が引かれるかどうかを確認
            drawn_ball = np.random.choice(balls)
            if drawn_ball == 1:
                count_Xi_1 += 1
        
        # シミュレーション結果の確率
        P_Xi_1_simulation.append(count_Xi_1 / n_trials)
    
    # グラフに追加描画
    plt.plot(M_values, P_Xi_1_theory, label=f'理論値 (N={N})', linestyle='--')
    plt.plot(M_values, P_Xi_1_simulation, label=f'シミュレーション (N={N})', marker='o', linestyle='None')

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

# 2021 Q2(2)  2024.8.22

import numpy as np
from scipy.special import comb

# 尤度関数 L(N_A) を定義
def L(NA):
    return comb(NA, 4) * comb(100 - NA, 11) / comb(100, 15)

# L(N_A + 1) / L(N_A) を計算し、条件を満たす最小の N_A を探索
# max(0, N_A - 85) ≤ 4 ≤ min(15, N_A) の式に基づき、取り出した豆Aが4粒の場合の N_A の範囲を設定 (4 ≤ N_A ≤ 89)
NA_values = np.arange(4, 90)
ratios = [(L(NA + 1) / L(NA)) for NA in NA_values]
NA_optimal = NA_values[np.where(np.array(ratios) < 1)[0][0]]

print(f"L(N_A + 1) / L(N_A) < 1 となる最小の N_A は {NA_optimal} です。")

L(N_A + 1) / L(N_A) < 1 となる最小の N_A は 26 です。

# 2021 Q2(2)  2024.8.22

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import comb

# 尤度関数 L(N_A) を定義
def L(NA):
    return comb(NA, 4) * comb(100 - NA, 11) / comb(100, 15)

# N_A の範囲を設定
NA_values = np.arange(4, 90)
ratios = [(L(NA + 1) / L(NA)) for NA in NA_values]

# グラフの描画
plt.figure(figsize=(10, 6))
plt.plot(NA_values, ratios, 'bo-', label=r'$\frac{L(N_A + 1)}{L(N_A)}$', markersize=8)
plt.axhline(y=1, color='red', linestyle='--', label=r'$1$')
plt.xlabel(r'$N_A$')
plt.ylabel(r'$\frac{L(N_A + 1)}{L(N_A)}$')
plt.title(r'$N_A$ の尤度比 $\frac{L(N_A + 1)}{L(N_A)}$')
plt.legend()
plt.grid(True)
plt.show()

# 2021 Q2(2)  2024.8.22

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import comb

# 尤度関数 L(N_A) を定義
def L(NA):
    return comb(NA, 4) * comb(100 - NA, 11) / comb(100, 15)

# N_A の範囲を設定
NA_values = np.arange(4, 90)
ratios = [(L(NA + 1) / L(NA)) for NA in NA_values]

# 26付近のズームしたグラフを描画
plt.figure(figsize=(10, 6))
plt.plot(NA_values, ratios, 'bo-', label=r'$\frac{L(N_A + 1)}{L(N_A)}$', markersize=8)
plt.axhline(y=1, color='red', linestyle='--', label=r'$1$')
plt.xlabel(r'$N_A$')
plt.ylabel(r'$\frac{L(N_A + 1)}{L(N_A)}$')
plt.title(r'$N_A$ の尤度比 $\frac{L(N_A + 1)}{L(N_A)}$ (25 < $N_A$ < 30)')
plt.xlim(25, 30)  # N_A = 26 付近をズーム
plt.ylim(0.95, 1.05)  # y軸もズームして見やすく
plt.legend()
plt.grid(True)
plt.show()