Bradley-Terryモデル - データサイエンス関連+αのメモ

Bradley-Terryモデル（Bradley & Terry, 1952）は、複数の項目（プレイヤー・製品・案など）を 総当たりのペア比較（pairwise comparison） にかけたときの勝敗データから、各項目の相対的な「強さ」を推定する統計モデル。

スポーツの対戦成績からチームの強さを推定する、複数のデザイン案を2つずつ比較させてどちらが良いか回答してもらう、といった場面で使われる。前章で見たとおり、これは選択型コンジョイント分析（CBC）と多項ロジットモデルにおける多項ロジットモデルを選択肢数2（ペア比較）に限定した特殊ケースに一致する。

モデル¶

$J$ 個の項目それぞれに強さのパラメータ $\pi_j > 0$ （あるいは対数強さ $\beta_j = \log \pi_j$ ）を割り当てる。項目 $i$ と項目 $j$ を比較したとき、 $i$ が勝つ確率を

P(i \succ j) = \frac{\pi_i}{\pi_i + \pi_j} = \frac{\exp(\beta_i)}{\exp(\beta_i) + \exp(\beta_j)}

(1)

と定義する。この式は選択型コンジョイント分析（CBC）と多項ロジットモデルで導出した多項ロジットモデルの選択確率

P(y = j \mid \{j,k\}) = \frac{\exp(V_j)}{\exp(V_j) + \exp(V_k)}

(2)

において $V_j = \beta_j$ としたものと完全に同じ形をしている。パラメータは全体の水準（スケール）が不定なので、通常はどれか1項目の $\beta$ を0に固定するか、 $\sum_j \beta_j = 0$ という制約をおいて識別する。

推定：最尤法¶

項目 $i$ が項目 $j$ に勝った回数を $w_{ij}$ とすると、対数尤度は

\ell(\boldsymbol{\beta}) = \sum_{i < j} \Big[ w_{ij} \log P(i \succ j) + w_{ji} \log P(j \succ i) \Big]

(3)

これを最大化する $\boldsymbol{\beta}$ は解析的に閉じた解を持たないため、通常はIRLS（反復重み付き最小二乗法）やニュートン法などの数値最適化、あるいは Zermeloのアルゴリズム（MM algorithm） と呼ばれる専用の反復更新式で求める。

対数尤度は、項目 $i$ の総勝利数を $W_i = \sum_{j \neq i} w_{ij}$ 、項目 $i,j$ 間の総比較回数を $n_{ij} = w_{ij} + w_{ji}$ とおくと（ $\pi_j = \exp(\beta_j)$ を用いて）

\ell(\boldsymbol{\pi}) = \sum_i W_i \log \pi_i - \sum_{i<j} n_{ij} \log(\pi_i + \pi_j)

(4)

と書き直せる。この式を直接最大化するのは難しいが、右辺第2項の $-\log(\cdot)$ が凹関数であることを利用すると、現在の推定値のまわりで各項を分離できる代理関数を作ることができる。

具体的には、対数関数の凹性から任意の $x, x_0 > 0$ について接線不等式

-\log x \ge -\log x_0 - \frac{x - x_0}{x_0}

(5)

が成り立つ。これを $x = \pi_i + \pi_j$ 、 $x_0 = \pi_i^{(t)} + \pi_j^{(t)}$ として各ペアに適用すると、対数尤度 $\ell(\boldsymbol\pi)$ を下から抑え、かつ $\pi_i^{(t)}$ の点で接する（＝その点で値が一致する）代理関数 $Q(\boldsymbol\pi \mid \boldsymbol\pi^{(t)})$ が得られる。この $Q$ は各 $\pi_i$ について分離した和になっているのが重要な点で、 $\ell(\pi^{(t)}) = Q(\pi^{(t)}\mid\pi^{(t)})$ であることから、 $Q$ を最大化する $\pi^{(t+1)}$ を選べば

\ell(\pi^{(t+1)}) \ge Q(\pi^{(t+1)}\mid\pi^{(t)}) \ge Q(\pi^{(t)}\mid\pi^{(t)}) = \ell(\pi^{(t)})

(6)

が成り立ち、真の対数尤度を単調に改善できる。この「扱いにくい目的関数を、下から近似する扱いやすい代理関数（minorizer）に置き換えて最大化を繰り返す」という一般原理がMM (Minorize-Maximize) algorithmであり、Zermelo(1929)がチェスの強さのレーティングを推定するために考案した反復式は、Hunter(2004)によってこのMMアルゴリズムの一例であることが示された。

$Q$ は $\pi_i$ ごとに分離しているため、各 $\pi_i$ についての1変数最大化問題を解くだけでよく、微分して0とおくと次の閉じた形の更新式が得られる。

\pi_i^{(t+1)} = \frac{W_i}{\displaystyle\sum_{j \neq i} \dfrac{n_{ij}}{\pi_i^{(t)} + \pi_j^{(t)}}}

(7)

これは後述の実装のfit_bradley_terry関数における

分子のwin_counts[i] → $W_i$ （項目 $i$ の総勝利数）
分母のn_ij / (pi[i] + pi[j])の総和 → $\sum_{j \neq i} n_{ij}/(\pi_i^{(t)}+\pi_j^{(t)})$

にそのまま対応している。

実装例¶

5つの項目（施策案A〜E）を総当たりでペア比較したデータから、Bradley-Terryモデルで各項目の強さを推定する。

import numpy as np
import pandas as pd
import itertools

rng = np.random.default_rng(0)

items = ["案A", "案B", "案C", "案D", "案E"]
true_beta = {"案A": 1.5, "案B": 0.5, "案C": 0.0, "案D": -0.5, "案E": -1.5}

# 総当たりペアそれぞれについて、n_trials回の比較を行う
n_trials = 20
records = []
for i, j in itertools.combinations(items, 2):
    p_i_wins = np.exp(true_beta[i]) / (np.exp(true_beta[i]) + np.exp(true_beta[j]))
    wins_i = rng.binomial(n_trials, p_i_wins)
    records.append({"winner": i, "loser": j, "count": wins_i})
    records.append({"winner": j, "loser": i, "count": n_trials - wins_i})

df = pd.DataFrame(records)
df.head()

Loading...

def fit_bradley_terry(df, items, n_iter=200):
    # Zermelo(MM algorithm)によるBradley-Terryモデルの推定
    pi = {item: 1.0 for item in items}  # 強さ pi_j = exp(beta_j) を1で初期化

    win_counts = {item: 0 for item in items}
    for _, row in df.iterrows():
        win_counts[row["winner"]] += row["count"]

    for _ in range(n_iter):
        new_pi = {}
        for i in items:
            numerator = win_counts[i]
            denominator = 0.0
            for j in items:
                if j == i:
                    continue
                n_ij = df.query("winner == @i and loser == @j")["count"].sum() \
                     + df.query("winner == @j and loser == @i")["count"].sum()
                denominator += n_ij / (pi[i] + pi[j])
            new_pi[i] = numerator / denominator
        # 識別のため合計を1に正規化
        total = sum(new_pi.values())
        pi = {k: v / total for k, v in new_pi.items()}

    return pi

pi_hat = fit_bradley_terry(df, items)
beta_hat = {k: np.log(v) for k, v in pi_hat.items()}

# 真値と比較するため、案Cを基準(0)にそろえる
offset = beta_hat["案C"]
beta_hat_centered = {k: v - offset for k, v in beta_hat.items()}

pd.DataFrame({"true_beta": true_beta, "estimated_beta": beta_hat_centered})

Loading...

基準を「案C」（真値0）にそろえた推定値 $\hat\beta$ は、真の強さの順序（A > B > C > D > E）とおおよその大きさを正しく復元できている。

AHPとの違い¶

Bradley-Terryモデルは各ペアの 勝敗（二値） または勝率という確率的な観測データから強さを推定する統計モデルであるのに対し、AHP（階層分析法）は各ペアについて回答者自身が「どちらがどれくらい優れているか」を1〜9段階などで直接評定した値を使い、固有値法で重みを計算する。Bradley-Terryは尤度に基づく統計的推定（標準誤差や信頼区間を計算できる）である一方、AHPは決定論的な行列計算である点が大きな違いである。