# 2021 Q2(4)  2024.8.24

import numpy as np
from scipy.special import comb

# 事前分布 P(N_A)
def prior(NA):
    return NA + 1

# 尤度関数 P(X = 4 | N_A)
def likelihood(NA):
    return comb(NA, 4) * comb(100 - NA, 11) / comb(100, 15)

# 事後分布 P(N_A | X = 4) (正規化定数は省略)
def posterior(NA):
    return prior(NA) * likelihood(NA)

# N_A の範囲
NA_values = np.arange(4, 90)

# 事後分布の計算
posterior_values = [posterior(NA) for NA in NA_values]

# 最大値を取る N_A (事後モード) を探索
NA_mode = NA_values[np.argmax(posterior_values)]

NA_mode

30

# 2021 Q2(4)  2024.8.24

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

# 事前分布 P(N_A)
def prior(NA):
    return NA + 1

# 正規化定数
C = 1 / 5151

# 尤度関数 P(X = 4 | N_A)
def likelihood(NA):
    return comb(NA, 4) * comb(100 - NA, 11) / comb(100, 15)

# 事後分布 P(N_A | X = 4)
def posterior(NA):
    return prior(NA) * likelihood(NA)

# N_A の範囲を0から100まで（事前分布用）
NA_values_full = np.arange(0, 101)

# 事前分布を計算（0から100まで）
prior_values_full = [C * prior(NA) for NA in NA_values_full]

# N_A の範囲を4から89まで（事後分布用）
NA_values = np.arange(4, 90)

# 事後分布の計算
posterior_values = [posterior(NA) for NA in NA_values]

# 事後分布の正規化
posterior_sum = sum(posterior_values)
posterior_values_normalized = [value / posterior_sum for value in posterior_values]

# 事前分布が存在しない場合の事後分布（尤度のみを正規化）
likelihood_values = [likelihood(NA) for NA in NA_values]
likelihood_sum = sum(likelihood_values)
likelihood_values_normalized = [value / likelihood_sum for value in likelihood_values]

# 3つの分布を重ねて表示
plt.figure(figsize=(10, 6))

# 事前分布のプロット（0から100まで）
plt.plot(NA_values_full, prior_values_full, 'g-', label=r'事前分布 $P(N_A)$', markersize=5)

# 事後分布のプロット（4から89まで）
plt.plot(NA_values, posterior_values_normalized, 'bo-', label=r'事後分布 $P(N_A | X = 4)$', markersize=5)

# 事前分布が存在しない場合の事後分布（尤度のみ）
plt.plot(NA_values, likelihood_values_normalized, 'r-', label=r'事前分布なし $P(N_A | X = 4)$', markersize=5)

# 事後モードのプロット
NA_mode = NA_values[np.argmax(posterior_values_normalized)]
plt.axvline(NA_mode, color='blue', linestyle='--', label=f'事後モード: {NA_mode}')

# 事前分布なしの場合の最尤推定値のプロット
NA_mode_likelihood = NA_values[np.argmax(likelihood_values_normalized)]
plt.axvline(NA_mode_likelihood, color='red', linestyle=':', label=f'最尤推定値 (事前分布なし): {NA_mode_likelihood}')

# グラフの設定
plt.xlabel(r'$N_A$')
plt.ylabel(r'確率')
plt.title(r'事前分布ありとなしの比較')
plt.grid(True)
plt.legend()
plt.show()