「自然科学の統計学」第5章演習問題6のブラッドリー・テリーのモデルへの$\chi^2$適合度検定をPythonで解いてみた

import numpy as np 
import pandas as pd

df = pd.read_csv('Bradley_Terrys_Model.csv', index_col = 0)
display(df)

# チームまたは選手数を取得
N = len(df)

# 全チームの強さの合計値を決定
K = 100

# 初期thetaを決定
thetas = [ K/N for _ in range(N)] 

print('初期theta:', thetas)

def calc_bradly(thetas, dataframe):
    # 元の強さold_thetasと新しい強さthetasの差が0.000001未満に収束するまで続ける
    while True:
        old_thetas = thetas.copy()
        for i in range(len(thetas)):
            #  式①を計算するための変数sum_
            sum_ = 0
            for j in range(len(thetas)):
                # 同じチーム同士の場合はpass
                if i == j:
                    pass
                else:
                    # (対戦試合数)/(2チームの強さの和)
                    sum_ += (df.iloc[i, j] + df.iloc[j, i]) / (thetas[i] + thetas[j])
                    print(f"({i},{j})  : {sum_}")
            # チームiの総勝利数を計算
            y_sum = df.iloc[i, :].sum()
            # 強さthetaを更新
            thetas[i] = y_sum / sum_
        # 元の強さold_thetasと新しい強さthetasの差を計算
        diff = np.abs(np.array(old_thetas) - np.array(thetas)).sum()
        print('diff = ', diff)
        if abs(diff) <= 0.0000001:
            # 収束した時のthetas
            print("収束")
            return thetas
        else:
            # 合計が100になるように標準化
            thetas = thetas / np.sum(thetas)
            thetas *= 100
            print('thetas: ', thetas)
    
thetas = calc_bradly(thetas, df)
print(thetas)

# 推定結果を入れるためのDataFrameの作成
tmp = pd.DataFrame(index = df.columns, columns =df.columns)

for i in range(N):
    for j in range(N):
        if i == j:
            pass
        else:
            # 2つのチームの合計試合数
            sum_ = df.iloc[i, j] + df.iloc[j, i]
            # ブラッドリー・テリーのモデル
            p = thetas[i] / (thetas[i] + thetas[j])
            # 勝ち数を推定
            tmp.iloc[i, j] = p * sum_
display(tmp)
display(df)

import scipy.stats as sps
import math

chi = 0
for i in range(N):
    for j in range(N):
        if i == j:
            pass
        else:
            real = df.iloc[i, j]
            inf = tmp.iloc[i, j]
            chi += (real - inf)**2 / inf
print("カイ2乗値: ", chi)
v = math.comb(N -1, 2)
if sps.chi2.ppf(q = 0.99, df = v) <= chi:
    print('1%で統計的に有意であり、ブラッドリー・テリーのモデルは適合しない')
elif sps.chi2.ppf(q = 0.95, df = v) <= chi:
    print('5%で統計的に有意であり、ブラッドリー・テリーのモデルは適合しない')
else:
    print('5%で統計的に有意でなく、ブラッドリー・テリーのモデルは適合する')