# 2014 Q4(2)  2025.1.10

import numpy as np

# シミュレーション設定
np.random.seed(42)  # 再現性のための乱数シード
n_simulations = 10000  # シミュレーション回数
sigma_squared = 1  # 誤差分散
true_theta = np.array([10, 20, 30, 40, 50])  # 真の値 (任意設定)

# (1) 各物体を2回測定する場合
X1 = np.array([
    [1, 0, 0, 0, 0],
    [1, 0, 0, 0, 0],
    [0, 1, 0, 0, 0],
    [0, 1, 0, 0, 0],
    [0, 0, 1, 0, 0],
    [0, 0, 1, 0, 0],
    [0, 0, 0, 1, 0],
    [0, 0, 0, 1, 0],
    [0, 0, 0, 0, 1],
    [0, 0, 0, 0, 1]
])

# (2) 2つずつ異なる組み合わせで測定する場合
X2 = np.array([
    [1, 1, 0, 0, 0],
    [1, 0, 1, 0, 0],
    [1, 0, 0, 1, 0],
    [1, 0, 0, 0, 1],
    [0, 1, 1, 0, 0],
    [0, 1, 0, 1, 0],
    [0, 1, 0, 0, 1],
    [0, 0, 1, 1, 0],
    [0, 0, 1, 0, 1],
    [0, 0, 0, 1, 1]
])

# 理論分散
theoretical_variance_1 = sigma_squared / 2
theoretical_variance_2 = 7 * sigma_squared / 24

# シミュレーション結果を格納するリスト
variances_1 = []
variances_2 = []

for _ in range(n_simulations):
    # 誤差を生成
    epsilon = np.random.normal(0, np.sqrt(sigma_squared), 10)

    # 測定値を生成
    x1 = X1 @ true_theta + epsilon  # (1) の場合
    x2 = X2 @ true_theta + epsilon  # (2) の場合

    # 最小二乗推定量を計算
    theta_hat_1 = np.linalg.inv(X1.T @ X1) @ X1.T @ x1
    theta_hat_2 = np.linalg.inv(X2.T @ X2) @ X2.T @ x2

    # 推定量の分散を計算
    variances_1.append(np.mean((theta_hat_1 - true_theta)**2))
    variances_2.append(np.mean((theta_hat_2 - true_theta)**2))

# シミュレーション結果の平均
simulated_variance_1 = np.mean(variances_1)
simulated_variance_2 = np.mean(variances_2)

# 結果の出力
print("(1) 理論分散:", theoretical_variance_1)
print("(1) シミュレーション分散:", simulated_variance_1)
print("(2) 理論分散:", theoretical_variance_2)
print("(2) シミュレーション分散:", simulated_variance_2)

(1) 理論分散: 0.5
(1) シミュレーション分散: 0.5003499949049185
(2) 理論分散: 0.2916666666666667
(2) シミュレーション分散: 0.2926572303618685