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

# 2021 Q2(1)  2024.8.21

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

# パラメータの設定
M = 100        # 全体の豆の数
N_A = 30       # 豆Aの数
n = 15         # 抽出する豆の数
x_values = np.arange(0, n+1)  # 可能な豆Aの数

# 理論的な超幾何分布のPMFを計算
rv = hypergeom(M, N_A, n)
pmf_theoretical = rv.pmf(x_values)

# 数値シミュレーションの設定
n_simulations = 10000  # シミュレーション回数
simulated_counts = []

# シミュレーションの実行
for _ in range(n_simulations):
    # 袋の中の豆を表すリスト（1が豆A、0が豆B）
    bag = np.array([1]*N_A + [0]*(M - N_A))
    # 無作為に15個抽出
    sample = np.random.choice(bag, size=n, replace=False)
    # 抽出した中の豆Aの数をカウント
    count_A = np.sum(sample)
    simulated_counts.append(count_A)

# シミュレーションから得られたPMFを計算
pmf_simulated, bins = np.histogram(simulated_counts, bins=np.arange(-0.5, n+1.5, 1), density=True)

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

# 理論的なPMFの描画
plt.plot(x_values, pmf_theoretical, 'bo-', label='理論的PMF', markersize=8)

# シミュレーション結果をヒストグラムとして描画
plt.hist(simulated_counts, bins=np.arange(-0.5, n+1.5, 1), density=True, alpha=0.5, color='red', label='シミュレーション結果')

# グラフの設定
plt.xlabel('豆Aの数')
plt.ylabel('確率')
plt.title('超幾何分布のPMF: 理論とシミュレーションの比較')
plt.legend()
plt.grid(True)
plt.show()