# 2022 Q5(3)  2024.8.15

import numpy as np
import pandas as pd
import scipy.stats as stats

# シミュレーションのための設定
np.random.seed(42)
n = 10  # 因子Aの水準数
mu_A = [5, 10, 15]  # 因子Aの水準ごとの平均
mu_B = [7, 12]  # 因子Bの水準ごとの平均
sigma = 2  # 標準偏差
n_samples = n * 2  # 各水準ごとのサンプル数

# データ生成
data = np.zeros((n, 2))
for i in range(n):
    data[i, 0] = np.random.normal(mu_A[i % len(mu_A)], sigma)  # 因子A
    data[i, 1] = np.random.normal(mu_B[i % len(mu_B)], sigma)  # 因子B

# データフレームにまとめる
df = pd.DataFrame(data, columns=['A', 'B'])

# 全平均
m = df.mean().mean()

# 各水準での平均
m_Ai = df['A'].mean()
m_Bj = df['B'].mean()

# 偏差平方和の計算
SA = np.sum((df['A'] - m_Ai) ** 2)
SB = np.sum((df['B'] - m_Bj) ** 2)
SE = np.sum((df.values - m_Ai - m_Bj + m) ** 2)

# 自由度
phi_A = n - 1
phi_B = 1
phi_E = n - 1

# 分散分析の統計量
F_A = (SA / phi_A) / (SE / phi_E)
F_B = (SB / phi_B) / (SE / phi_E)

# p値の計算
p_A = 1 - stats.f.cdf(F_A, phi_A, phi_E)
p_B = 1 - stats.f.cdf(F_B, phi_B, phi_E)

# 結果の表示
{
    "F_A": F_A,
    "p_A": p_A,
    "F_B": F_B,
    "p_B": p_B,
    "SA": SA,
    "SB": SB,
    "SE": SE,
    "df_A": phi_A,
    "df_B": phi_B,
    "df_E": phi_E
}

{'F_A': 0.5528514017044434,
 'p_A': 0.8047382043674255,
 'F_B': 4.023821861099081,
 'p_B': 0.07582222108466219,
 'SA': 139.17265361165713,
 'SB': 112.54902298706004,
 'SE': 251.73609614190576,
 'df_A': 9,
 'df_B': 1,
 'df_E': 9}