# 2015 Q2(4)  2024.12.10

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

# 設定
alpha = 0.05  # 有意水準
target_power = 0.8  # 目標検出力
mu = 0.5  # 真の平均
z_alpha = norm.ppf(1 - alpha)  # 臨界値 z_alpha
z_power = norm.ppf(target_power)  # 検出力に対応する Z 値

# 理論的な検出力の計算
n_values = np.arange(1, 51, 1)  # 標本の大きさ n = 1, 2, ..., 50
powers_theoretical = [1 - norm.cdf(z_alpha - mu * np.sqrt(n)) for n in n_values]

# 実際の検出力を乱数シミュレーションで計算
num_simulations = 10000  # シミュレーション回数
n_sim_values = np.arange(5, 51, 5)  # n = 5, 10, ..., 50
powers_simulated = []

for n in n_sim_values:
    # 標本平均を生成
    samples = np.random.normal(mu, 1 / np.sqrt(n), num_simulations)
    # 標本平均が臨界値を超えた割合を計算
    rejection_rate = np.mean(samples > z_alpha / np.sqrt(n))
    powers_simulated.append(rejection_rate)

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

# 理論的な検出力の曲線
plt.plot(n_values, powers_theoretical, label="理論的な検出力 (Power)", color="blue")

# シミュレーションによる検出力
plt.scatter(n_sim_values, powers_simulated, color="red", label="シミュレーションによる検出力", zorder=5)

# 検出力の目標ラインと必要な標本サイズライン
plt.axhline(y=target_power, color="red", linestyle="--", label=f"目標検出力 Power = {target_power}")
plt.axvline(x=25, color="green", linestyle="--", label="必要な標本数 n = 25")

# グラフの設定
plt.xlabel("標本の大きさ $n$", fontsize=12)
plt.ylabel("検出力 (Power)", fontsize=12)
plt.title("理論値とシミュレーションによる検出力の比較", fontsize=14)
plt.legend()
plt.grid(True)
plt.show()

# 2015 Q2(3)  2024.12.9

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

# 定義
alpha = 0.05  # 有意水準
z_alpha = norm.ppf(1 - alpha)  # 臨界値 z_alpha
sample_sizes = [9, 16]  # n = 9, 16
true_means = [0.4, 0.8]  # 真の平均 μ = 0.4, 0.8

# 検出力の計算関数
def compute_power(mu, n, z_alpha):
    z_value = z_alpha - mu * np.sqrt(n)
    power = 1 - norm.cdf(z_value)
    return power

# 検出力を計算
mu_values = np.linspace(0, 1.5, 100)  # 真の平均 μ の範囲
powers_n9 = [compute_power(mu, n=9, z_alpha=z_alpha) for mu in mu_values]
powers_n16 = [compute_power(mu, n=16, z_alpha=z_alpha) for mu in mu_values]

# 指定された μ = 0.4, 0.8 の検出力を計算
specific_powers = {
    (n, mu): compute_power(mu, n, z_alpha)
    for n in sample_sizes
    for mu in true_means
}

# グラフの描画
plt.figure(figsize=(10, 6))
plt.plot(mu_values, powers_n9, label="検出力 (n=9)", linestyle="--", color="blue")
plt.plot(mu_values, powers_n16, label="検出力 (n=16)", linestyle="-", color="green")

# μ = 0.4, 0.8 の場合の検出力をプロット
for (n, mu), power in specific_powers.items():
    plt.scatter(mu, power, label=f"n={n}, μ={mu}, Power={power:.2f}", zorder=5)

# グラフの設定
plt.xlabel("真の平均 $\\mu$", fontsize=12)
plt.ylabel("検出力 (Power)", fontsize=12)
plt.title("検出力の挙動 ($n=9$ および $n=16$ の場合)", fontsize=14)
plt.axhline(y=alpha, color="red", linestyle=":", label=f"有意水準 α={alpha}")
plt.legend()
plt.grid(True)
plt.show()