# 2022 Q5(5)  2024.8.17

import numpy as np

# シミュレーションの設定
np.random.seed(42)
n = 5  # サンプル数
X = np.random.normal(10, 2, n)  # Xの観測値

# Y1が非常に大きな値を取るケースを考慮
Y = np.random.normal(15, 2, n)
Y[0] = 1000  # Y1が非常に大きな値を取る設定

# 全体平均とY2:nの平均
Y_bar = np.mean(Y)  # 全体の平均
Y_bar_2n = np.mean(Y[1:])  # Y2:n の平均
X_bar_2n = np.mean(X[1:])  # X2:n の平均

# θの推定値
theta_hat = Y_bar_2n - X_bar_2n

# SAの計算（Y1を含める場合）
SA_with_Y1 = np.sum(((X[1:] + Y[1:] - theta_hat) / 2 - X_bar_2n) ** 2) + (Y[0] - Y_bar) ** 2
SA_with_Y1 += np.sum(((X[1:] + Y[1:] + theta_hat) / 2 - Y_bar) ** 2)

# SAの計算（Y1を除く場合）
SA_without_Y1 = np.sum(((X[1:] + Y[1:] - theta_hat) / 2 - X_bar_2n) ** 2)
SA_without_Y1 += np.sum(((X[1:] + Y[1:] + theta_hat) / 2 - Y_bar_2n) ** 2)

# 結果を表示
{
    "Y1を含めた場合の要因A平方和 (SA)": SA_with_Y1,
    "Y1を除いた場合の要因A平方和 (SA)": SA_without_Y1
}

{'Y1を含めた場合の要因A平方和 (SA)': 774276.0933538243,
 'Y1を除いた場合の要因A平方和 (SA)': 1.6720954914667048}