# 2014 Q3(2)  2025.1.7

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

# Parameters
n = 30  # サンプルサイズ
beta_true = 1.5  # 真のβ
sigma_X = 1.0  # Xの標準偏差
sigma_Y = 3.0  # Yの標準偏差
mu_X = 5.0  # Xの平均
alpha = 0.05  # 有意水準
num_simulations = 100  # シミュレーション回数

# 分散未知の場合の信頼区間を計算する関数
def confidence_interval_unknown_variance(X, Y, beta_0, alpha):
    X_bar = np.mean(X)
    Y_bar = np.mean(Y)
    U = Y - beta_0 * X
    S_square = np.sum((U - np.mean(U))**2) / (n - 1)  # 標本分散
    S_star = np.sqrt(S_square)
    t_alpha = t.ppf(1 - alpha / 2, df=n - 1)
    lower_bound = Y_bar / X_bar - t_alpha * S_star / (np.sqrt(n) * X_bar)
    upper_bound = Y_bar / X_bar + t_alpha * S_star / (np.sqrt(n) * X_bar)
    return lower_bound, upper_bound

# シミュレーション
intervals_unknown = []
contains_true_beta_unknown = []

for _ in range(num_simulations):
    X = np.random.normal(mu_X, sigma_X, n)
    epsilon = np.random.normal(0, np.sqrt(sigma_Y**2 + (beta_true * sigma_X)**2), n)
    Y = beta_true * X + epsilon

    # 信頼区間を計算
    ci = confidence_interval_unknown_variance(X, Y, beta_true, alpha)
    intervals_unknown.append(ci)
    contains_true_beta_unknown.append(ci[0] <= beta_true <= ci[1])

# カバー率を計算
cover_rate_unknown = sum(contains_true_beta_unknown) / num_simulations
print("カバー率:")
print(f"  分散未知の場合: {cover_rate_unknown:.4f}")    
    
# 可視化
plt.figure(figsize=(12, 8))
for i, ci in enumerate(intervals_unknown):
    color = 'green' if contains_true_beta_unknown[i] else 'red'
    plt.plot([ci[0], ci[1]], [i, i], color=color, marker='o', label="" if i > 0 else "信頼区間")
plt.axvline(x=beta_true, color='blue', linestyle='--', label=f"真のβ: {beta_true}")
plt.title("分散未知の場合の信頼区間 (各シミュレーション)", fontsize=16)
plt.xlabel("β", fontsize=12)
plt.ylabel("シミュレーション回数", fontsize=12)
plt.legend()
plt.grid()
plt.show()

カバー率:
  分散未知の場合: 0.9500