# 2015 Q4(1)  2024.12.17

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

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

# 総サンプルサイズ N の計算
N = observed_counts.sum()

# セル確率 pij の計算 (観測度数を総数で割る)
p_ij = observed_counts / N

# 観測度数とセル確率
row_labels = ["Highest", "Second", "Third", "Lowest"]
col_labels = ["Highest", "Second", "Third", "Lowest"]

# 観測度数
df_counts = pd.DataFrame(observed_counts, columns=col_labels, index=row_labels)
# セル確率
df_probs = pd.DataFrame(np.round(p_ij, 6), columns=col_labels, index=row_labels)

print("観測度数表 (x_ij):")
print(df_counts)

print("\nセル確率表 (p_ij):")
print(df_probs)

# 観測度数が得られる確率を計算
x_flat = observed_counts.flatten()
p_flat = p_ij.flatten()

# 多項分布の確率を計算
prob = multinomial.pmf(x_flat, n=N, p=p_flat)
print(f"\nこの観測度数表が得られる確率: {prob:.10e}")

観測度数表 (x_ij):
         Highest  Second  Third  Lowest
Highest     1520     266    124      66
Second       234    1512    432      78
Third        117     362   1772     205
Lowest        36      82    179     492

セル確率表 (p_ij):
          Highest    Second     Third    Lowest
Highest  0.203290  0.035576  0.016584  0.008827
Second   0.031296  0.202220  0.057777  0.010432
Third    0.015648  0.048415  0.236993  0.027417
Lowest   0.004815  0.010967  0.023940  0.065802

この観測度数表が得られる確率: 7.0375704386e-24