# 2021 Q2(3)  2024.8.23

from sympy import symbols, summation, Eq, solve

# 記号の定義
n = symbols('n', integer=True)
C = symbols('C', real=True)

# 和を計算
sum_expr = summation(C * (n + 1), (n, 0, 100))

# 正規化条件を設定
normalization_condition = Eq(sum_expr, 1)

# C を解く
C_value = solve(normalization_condition, C)[0]
C_value

# 2021 Q2(3)  2024.8.23

import matplotlib.pyplot as plt
import numpy as np

# 正規化定数 C
C_value = 1 / 5151

# n の範囲
n_values = np.arange(0, 101)

# 確率質量関数 P(N_A = n)
pmf_values = C_value * (n_values + 1)

# グラフの描画
plt.figure(figsize=(10, 6))
plt.plot(n_values, pmf_values, 'bo-', label=r'$P(N_A = n) = \frac{1}{5151} \cdot (n + 1)$', markersize=5)
plt.xlabel(r'$n$')
plt.ylabel(r'$P(N_A = n)$')
plt.title(r'確率質量関数 $P(N_A = n)$ のグラフ')
plt.grid(True)
plt.legend()
plt.show()

# 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}

# 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}

# 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}