# 2015 Q4(3)  2024.12.19

import numpy as np
from scipy.stats import chi2

# 観測度数表
observed_counts = np.array([[1520, 266, 124, 66],
                            [234, 1512, 432, 78],
                            [117, 362, 1772, 205],
                            [36, 82, 179, 492]])

I = observed_counts.shape[0]  # 行・列のサイズ

# 尤度比統計量 G^2 の計算 (対角要素を除外)
G2 = 2 * sum(
    observed_counts[i, j] * np.log(2 * observed_counts[i, j] / (observed_counts[i, j] + observed_counts[j, i]))
    for i in range(I) for j in range(I) if i != j
)

# 自由度の計算 (非対角要素の数: I(I-1)/2)
df = (I * (I - 1)) // 2

# P値の計算
p_value = 1 - chi2.cdf(G2, df)

# 結果の表示
print(f"尤度比統計量 G^2: {G2:.4f}")
print(f"自由度: {df}")
print(f"P値: {p_value:.10f}")

# 帰無仮説の判定
if p_value < 0.01:
    print("P値は0.01より小さいため、帰無仮説 H_0 は棄却されます。")
else:
    print("P値は0.01以上のため、帰無仮説 H_0 は棄却されません。")

尤度比統計量 G^2: 19.2492
自由度: 6
P値: 0.0037628518
P値は0.01より小さいため、帰無仮説 H_0 は棄却されます。