# 2018 Q1(2)  2024.10.2

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
sigma = 2   # 母分散の標準偏差
sigma_squared = sigma ** 2  # 母分散
n = 30      # サンプルサイズ
num_simulations = 10000  # シミュレーションの回数

# 理論上の標本分散の分散 Var(S^2)
theoretical_variance = (2 * sigma_squared ** 2) / (n - 1)

# シミュレーション結果を保存するリスト
sample_variances = []

# シミュレーションを繰り返す
for _ in range(num_simulations):
    # 正規分布 N(0, sigma^2) に従うサンプルを生成
    sample = np.random.normal(0, sigma, n)
    
    # 標本平均
    sample_mean = np.mean(sample)
    
    # 標本分散 (S^2)
    sample_variance = np.sum((sample - sample_mean) ** 2) / (n - 1)
    
    # 計算した標本分散をリストに追加
    sample_variances.append(sample_variance)

# サンプル分散の分散を計算
empirical_variance = np.var(sample_variances)

# 結果表示
print(f"理論的な標本分散の分散: {theoretical_variance:.4f}")
print(f"シミュレーションで得られた標本分散の分散: {empirical_variance:.4f}")

# ヒストグラムを描画して標本分散の分布を確認
plt.hist(sample_variances, bins=50, alpha=0.7, color='b', edgecolor='black')
plt.axvline(np.mean(sample_variances), color='r', linestyle='dashed', linewidth=2, label='平均分散')
plt.title('標本分散の分布')
plt.xlabel('標本分散 S^2')
plt.ylabel('頻度')
plt.legend()
plt.grid(True)
plt.show()

理論的な標本分散の分散: 1.1034
シミュレーションで得られた標本分散の分散: 1.1286