# 2015 Q2(2)  2024.12.8

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

# シミュレーション設定
n = 10  # 標本サイズ
mu_0 = 0  # 帰無仮説の平均
sigma = 1  # 標準偏差
alpha = 0.05  # 有意水準
num_simulations = 10000  # シミュレーション回数

# z_alpha (有意水準に対応する臨界値)
z_alpha = norm.ppf(1 - alpha)
threshold = z_alpha / np.sqrt(n)  # 棄却域の閾値

# シミュレーションで標本平均を生成
samples = np.random.normal(mu_0, sigma, (num_simulations, n))
x_bar = np.mean(samples, axis=1)

# 棄却域に入るかどうかを確認
in_rejection_region = x_bar > threshold

# P値の計算
p_values = 1 - norm.cdf(x_bar * np.sqrt(n) / sigma)

# P値が有意水準以下かどうかを確認
p_value_condition = p_values <= alpha

# 棄却域とP値の条件が一致する割合を計算
agreement_rate = np.mean(in_rejection_region == p_value_condition)

# 棄却域に落ちる割合とP値がα以下になる割合の比較
rejection_rate = np.mean(in_rejection_region)
p_value_rate = np.mean(p_value_condition)

# 結果の出力
print(f"シミュレーション回数: {num_simulations}")
print(f"有意水準 α: {alpha}")
print(f"棄却域条件とP値条件の一致率: {agreement_rate:.3f}")
print(f"棄却域に落ちる割合: {rejection_rate:.3f}")
print(f"P値がα以下になる割合: {p_value_rate:.3f}")

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

# 標本平均 x̄ の分布
plt.hist(x_bar, bins=50, color="skyblue", alpha=0.7, edgecolor="black", label="標本平均 $\\bar{X}$ の分布")
plt.axvline(threshold, color="red", linestyle="--", label=f"棄却域の閾値 $z_{{\\alpha}}/\\sqrt{{n}} = {threshold:.3f}$")
plt.fill_betweenx([0, plt.gca().get_ylim()[1]], threshold, plt.gca().get_xlim()[1], color="red", alpha=0.2, label="棄却域")
plt.xlabel("標本平均 $\\bar{X}$", fontsize=12)
plt.ylabel("度数", fontsize=12)
plt.title("標本平均 $\\bar{X}$ の分布と棄却域", fontsize=14)
plt.legend()
plt.grid(True)
plt.show()

# P値の分布
plt.figure(figsize=(12, 6))
plt.hist(p_values, bins=50, color="lightgreen", alpha=0.7, edgecolor="black", label="P値の分布")
plt.axvline(alpha, color="blue", linestyle="--", label=f"有意水準 $\\alpha = {alpha}$")
plt.fill_betweenx([0, plt.gca().get_ylim()[1]], 0, alpha, color="blue", alpha=0.2, label="P値が $\\leq \\alpha$ の領域")
plt.xlabel("P値", fontsize=12)
plt.ylabel("度数", fontsize=12)
plt.title("P値の分布と有意水準", fontsize=14)
plt.legend()
plt.grid(True)
plt.show()

シミュレーション回数: 10000
有意水準 α: 0.05
棄却域条件とP値条件の一致率: 1.000
棄却域に落ちる割合: 0.052
P値がα以下になる割合: 0.052