# 2014 Q5(2)  2025.1.16

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

# シミュレーションの設定
n = 100  # 人数
p_true = 0.6  # 真の成功確率
categories = np.arange(5)  # カテゴリ (0～4)
n_simulations = 1000  # シミュレーション回数

# 二項分布モデルに基づく q_j の計算 (忠実な設定)
qj_true = [np.math.comb(4, j) * (p_true ** j) * ((1 - p_true) ** (4 - j)) for j in categories]

# 二項分布に従わない (qj_noisy): q_j をランダムに変更してズラす
qj_noisy = np.array(qj_true) + np.random.uniform(-0.1, 0.1, size=len(qj_true))
qj_noisy = np.clip(qj_noisy, 0, None)  # 確率が0未満にならないようにする
qj_noisy = qj_noisy / np.sum(qj_noisy)  # 確率なので正規化

# カイ二乗統計量と p 値を記録するリスト
chi_squared_true = []
p_values_true = []
chi_squared_noisy = []
p_values_noisy = []

for _ in range(n_simulations):
    # 忠実なモデルの観測データ
    obs_true = np.random.multinomial(n, qj_true)
    exp_true = n * np.array(qj_true)
    chi2_stat_true = np.sum((obs_true - exp_true) ** 2 / exp_true)
    p_val_true = 1 - chi2.cdf(chi2_stat_true, df=3)  # 自由度 = カテゴリ数 - 1 - 推定パラメータ数
    chi_squared_true.append(chi2_stat_true)
    p_values_true.append(p_val_true)

    # ズレたモデルの観測データ
    obs_noisy = np.random.multinomial(n, qj_noisy)
    exp_noisy = n * np.array(qj_true)  # 忠実な期待値を使用
    chi2_stat_noisy = np.sum((obs_noisy - exp_noisy) ** 2 / exp_noisy)
    p_val_noisy = 1 - chi2.cdf(chi2_stat_noisy, df=3)
    chi_squared_noisy.append(chi2_stat_noisy)
    p_values_noisy.append(p_val_noisy)

# 平均カイ二乗統計量を計算
mean_chi2_true = np.mean(chi_squared_true)
mean_chi2_noisy = np.mean(chi_squared_noisy)

# カイ二乗分布のPDFを計算
x = np.linspace(0, 20, 500)  # x軸の範囲
pdf = chi2.pdf(x, df=3)  # 自由度3のPDF

# 結果の出力
print(f"忠実なモデルの平均カイ二乗統計量: {mean_chi2_true:.4f}")
print(f"忠実なモデルの平均p値: {np.mean(p_values_true):.4f}")
print(f"ズレたモデルの平均カイ二乗統計量: {mean_chi2_noisy:.4f}")
print(f"ズレたモデルの平均p値: {np.mean(p_values_noisy):.4f}")

# グラフの描画
plt.figure(figsize=(10, 6))
plt.plot(x, pdf, label="自由度3のカイ二乗分布", linewidth=2)
plt.axvline(mean_chi2_true, color="blue", linestyle="--", label=f"忠実モデルの平均: {mean_chi2_true:.2f}")
plt.axvline(mean_chi2_noisy, color="red", linestyle="--", label=f"ズレモデルの平均: {mean_chi2_noisy:.2f}")
plt.xlabel("$X^2$ 値", fontsize=12)
plt.ylabel("確率密度", fontsize=12)
plt.title("自由度3のカイ二乗分布と平均カイ二乗統計量", fontsize=14)
plt.legend(fontsize=12)
plt.grid(alpha=0.5)
plt.show()

忠実なモデルの平均カイ二乗統計量: 3.8112
忠実なモデルの平均p値: 0.4026
ズレたモデルの平均カイ二乗統計量: 13.4317
ズレたモデルの平均p値: 0.0236