# 2022 Q5(4)  2024.8.16

import numpy as np

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

# ケース1: Y1を含めた場合の計算
mu_0_1 = Y[0]  # Y1に対応する μ_0
mu_i_with_Y1 = (X[1:] + Y[1:]) / 2  # 残りの μ_i の計算
SE_with_Y1 = np.sum((X[1:] - mu_i_with_Y1) ** 2) + np.sum((Y[1:] - mu_i_with_Y1) ** 2)

# Y1に対する特別な処理
SE_with_Y1 += (Y[0] - mu_0_1) ** 2  # Y1 の影響を追加

# ケース2: Y1を除いた場合の計算
mu_i_without_Y1 = (X[1:] + Y[1:]) / 2  # Y1を除いた場合の μ_i の計算
SE_without_Y1 = np.sum((X[1:] - mu_i_without_Y1) ** 2) + np.sum((Y[1:] - mu_i_without_Y1) ** 2)

# 要因Bに関する平方和（SB）の計算
mu_B_with_Y1 = np.mean(Y[1:])  # Y1を含めた場合のBの平均
mu_B_without_Y1 = np.mean(Y[1:])  # Y1を除いた場合のBの平均

# Y1専用の平均を計算してY1の寄与を計算
mu_0_Y1 = Y[0]  # Y1専用の平均はY1自体
SB_with_Y1 = np.sum((Y[1:] - mu_B_with_Y1) ** 2)
SB_with_Y1 += (Y[0] - mu_0_Y1) ** 2  # Y1の寄与を専用の平均で追加

SB_without_Y1 = np.sum((Y[1:] - mu_B_without_Y1) ** 2)

# 結果を表示
{
    "Y1を含めた場合の残差平方和 (SE)": SE_with_Y1,
    "Y1を除いた場合の残差平方和 (SE)": SE_without_Y1,
    "Y1を含めた場合の要因B平方和 (SB)": SB_with_Y1,
    "Y1を除いた場合の要因B平方和 (SB)": SB_without_Y1,
}

{'Y1を含めた場合の残差平方和 (SE)': 71.28917221001626,
 'Y1を除いた場合の残差平方和 (SE)': 71.28917221001626,
 'Y1を含めた場合の要因B平方和 (SB)': 8.53547835010338,
 'Y1を除いた場合の要因B平方和 (SB)': 8.53547835010338}