Bradley-Terry Model › Code Implementation

Michael BrenndoerferDecember 23, 202514 min read
Language AI HandbookAll chapters
Part 37 of 67Alignment and RLHF
Chapter 274Bradley-Terry Model › Code Implementation
Reading section 8 of 12 · Code Implementation · 14 minSection 8 of 12 · 14 min

Book contents

Language AI Handbook67 parts · 513 chaptersView the whole book →
Active chapter 274 of 513

Section 8 of 12 · 14 min read

Bradley-Terry Model › Code ImplementationOne section · Sections in this chapter8 / 12 ↓
Sections in this chapter12
Reading tools Definitions & inline links

Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.

Article links

Make inline references clickable

We'll implement Bradley-Terry estimation from scratch using gradient descent, then show how it's applied in practice. The implementation closely mirrors the mathematics we derived above, which should make the connection between theory and code transparent.

In[9]:
Code
# Define our comparison data
# Each tuple: (winner_idx, loser_idx, count)
comparisons = [
    (0, 1, 7),  # A beats B, 7 times
    (1, 0, 3),  # B beats A, 3 times
    (1, 2, 8),  # B beats C, 8 times
    (2, 1, 2),  # C beats B, 2 times
    (0, 2, 9),  # A beats C, 9 times
    (2, 0, 1),  # C beats A, 1 time
]

n_items = 3  # A, B, C
item_names = ["A", "B", "C"]

Now let's define the negative log-likelihood function that we want to minimize. We negate the log-likelihood because standard optimizers minimize rather than maximize:

In[10]:
Code
def neg_log_likelihood(beta, comparisons):
    """
    Compute negative log-likelihood for Bradley-Terry model.
    beta[0] corresponds to item A, beta[1] to B, etc.
    Item C (index 2) is fixed at 0 for identifiability.
    """
    # Extend beta with the fixed anchor (beta_C = 0)
    beta_full = np.array([beta[0], beta[1], 0.0])

    nll = 0.0
    for winner, loser, count in comparisons:
        # Log probability that winner beats loser
        log_prob = np.log(sigmoid(beta_full[winner] - beta_full[loser]) + 1e-10)
        nll -= count * log_prob

    return nll

We also need the gradient for efficient optimization. The gradient tells the optimizer which direction to move in parameter space to improve the log-likelihood:

In[11]:
Code
def neg_log_likelihood_gradient(beta, comparisons):
    """
    Compute gradient of negative log-likelihood.
    """
    beta_full = np.array([beta[0], beta[1], 0.0])
    grad = np.zeros(2)  # Only for beta_A and beta_B (beta_C is fixed)

    for winner, loser, count in comparisons:
        # P(winner beats loser)
        prob = sigmoid(beta_full[winner] - beta_full[loser])

        # Gradient contribution
        # d/d(beta_winner) of log(sigma(x)) = 1 - sigma(x)
        # d/d(beta_loser) of log(sigma(x)) = -(1 - sigma(x))
        if winner < 2:  # Only update if not the anchor
            grad[winner] -= count * (1 - prob)
        if loser < 2:
            grad[loser] += count * (1 - prob)

    return grad

Now let's optimize to find the maximum likelihood estimates. We use L-BFGS-B, a quasi-Newton method that approximates the Hessian matrix and converges much faster than simple gradient descent for smooth convex objectives like ours:

In[12]:
Code
# Initial guess
beta_init = np.array([0.0, 0.0])

# Optimize using L-BFGS-B
result = minimize(
    neg_log_likelihood,
    beta_init,
    args=(comparisons,),
    method="L-BFGS-B",
    jac=neg_log_likelihood_gradient,
)

beta_A, beta_B = result.x
beta_C = 0.0  # Fixed anchor
Out[13]:
Console
Maximum Likelihood Estimates:
  β_A = 2.2156
  β_B = 1.3761
  β_C = 0.0000 (fixed anchor)

Optimization converged: True
Final negative log-likelihood: 14.3638

The estimated coefficients indicate that Item A is the strongest, followed by Item B, with Item C as the baseline. The optimization successfully converged to these maximum likelihood estimates. Item A's log-strength of roughly 2 units above C corresponds to a predicted 88% win rate against C, which is close to (but not exactly) the observed 90%. The Bradley-Terry model smooths the observed frequencies slightly, because it must simultaneously satisfy constraints from all three matchups rather than fitting each one independently.

Let's verify these estimates by computing the predicted win probabilities:

In[14]:
Code
# Compute predicted probabilities
beta_full = np.array([beta_A, beta_B, beta_C])

prob_A_beats_B = sigmoid(beta_A - beta_B)
prob_B_beats_C = sigmoid(beta_B - beta_C)
prob_A_beats_C = sigmoid(beta_A - beta_C)

# Calculate observed win rates dynamically from data
# comparisons is list of (winner, loser, count)
counts = {(w, l): c for w, l, c in comparisons}

obs_A_beats_B = counts[(0, 1)] / (counts[(0, 1)] + counts[(1, 0)])
obs_B_beats_C = counts[(1, 2)] / (counts[(1, 2)] + counts[(2, 1)])
obs_A_beats_C = counts[(0, 2)] / (counts[(0, 2)] + counts[(2, 0)])
Out[15]:
Console
Predicted Win Probabilities:
  P(A > B) = 0.698  (observed: 0.700)
  P(B > C) = 0.798  (observed: 0.800)
  P(A > C) = 0.902  (observed: 0.900)

The model closely matches the observed frequencies. Now let's convert these to Elo-style ratings for interpretability. The Elo scale is more familiar to practitioners and gives an intuitive sense of the magnitude of the quality differences:

In[16]:
Code
# Convert to Elo scale (where 400 points = 10:1 odds)
elo_scale = 400 / np.log(10)
elo_ratings = beta_full * elo_scale

# Normalize so average is 1500 (conventional Elo baseline)
elo_ratings = elo_ratings - np.mean(elo_ratings) + 1500
Out[17]:
Console
Elo-Style Ratings:
  A: 1677
  B: 1531
  C: 1292

The ratings align with our expectations: Item A has the highest rating. This reflects its dominance in the pairwise comparisons. A difference of approximately 200 points between A and B suggests a strong win probability for A, corresponding to roughly 76% expected wins in head-to-head matchups.

Visualizing the Preference Surface

Let's visualize how the Bradley-Terry probability varies with the strength difference, and where our estimated matchups fall on this curve:

Out[18]:
Visualization
S-shaped curve showing preference probability increasing from 0 to 1 as strength difference increases, with colored scatter points marking estimated differences for the A vs B, B vs C, and A vs C matchups.
Bradley-Terry preference probability curve as a function of strength difference $\beta_i - \beta_j$. The sigmoid function produces a logistic response where probabilities change most rapidly near zero difference and saturate asymptotically at the extremes. Colored points mark the estimated strength differences for the three matchups in our worked example.

The curve shows the characteristic sigmoid shape. When strengths are equal (βi−βj=0\beta_i - \beta_j = 0), the win probability is exactly 0.5. Small strength differences near zero produce probability changes close to linear, while large differences saturate toward 0 or 1. Notice that the A-vs-C point sits further to the right than the A-vs-B point (larger strength difference) and therefore has a higher predicted win probability.

Out[19]:
Visualization
Filled contour plot of the Bradley-Terry log-likelihood surface over a grid of beta_A and beta_B values, with a red star marking the maximum likelihood estimate and a viridis color scale showing log-likelihood magnitude.
Log-likelihood surface for the three-item worked example with $\beta_C=0$. The contours show how model fit varies with parameters $\beta_A$ and $\beta_B$, with the maximum likelihood estimate marked by the red star. The elongated shape reflects the correlation between parameters arising from the relative nature of pairwise comparisons.
Out[20]:
Visualization
Grouped bar chart comparing observed win frequencies (blue) and Bradley-Terry predicted probabilities (red) for the A vs B, B vs C, and A vs C matchups, with annotated values above each bar.
Observed win frequencies versus Bradley-Terry model predictions for the three pairwise matchups. The close alignment between empirical win rates (blue) and model probabilities (red) confirms the model's ability to fit the underlying preference structure. Small discrepancies arise because the model must be globally consistent rather than fitting each matchup independently.

Handling Real Preference Data

In practice, preference data comes as individual comparison records rather than aggregated counts. Large-scale LLM evaluation pipelines might receive thousands of individual annotation records, each showing which of two responses a specific annotator preferred to a specific prompt. Let's implement a more realistic version that handles raw comparison logs and scales to many items:

In[21]:
Code
class BradleyTerryModel:
    """
    Bradley-Terry model for pairwise comparisons.
    Estimates latent strength parameters from preference data.
    The last item (index n_items-1) is used as the fixed anchor at beta=0.
    """

    def __init__(self, n_items):
        self.n_items = n_items
        self.beta = np.zeros(n_items)  # Log-strength parameters
        self.fitted = False

    def _neg_log_likelihood(self, beta, winners, losers):
        """Compute NLL for optimization."""
        # Fix last item at 0 for identifiability
        beta_full = np.concatenate([beta, [0.0]])

        diffs = beta_full[winners] - beta_full[losers]
        # Use log1p for numerical stability: log(sigmoid(x)) = -log(1 + exp(-x))
        log_probs = -np.log1p(np.exp(-diffs))

        return -np.sum(log_probs)

    def _gradient(self, beta, winners, losers):
        """Compute gradient of NLL."""
        beta_full = np.concatenate([beta, [0.0]])
        diffs = beta_full[winners] - beta_full[losers]
        probs = sigmoid(diffs)

        grad = np.zeros(self.n_items - 1)
        for i in range(self.n_items - 1):
            # Contribution from comparisons where item i won
            win_mask = winners == i
            grad[i] -= np.sum(1 - probs[win_mask])

            # Contribution from comparisons where item i lost
            lose_mask = losers == i
            grad[i] += np.sum(1 - probs[lose_mask])

        return grad

    def fit(self, winners, losers):
        """
        Fit model to comparison data.

        Args:
            winners: Array of winner indices for each comparison
            losers: Array of loser indices for each comparison
        """
        winners = np.asarray(winners)
        losers = np.asarray(losers)

        beta_init = np.zeros(self.n_items - 1)

        result = minimize(
            self._neg_log_likelihood,
            beta_init,
            args=(winners, losers),
            method="L-BFGS-B",
            jac=lambda b, w, l: self._gradient(b, w, l),
        )

        self.beta = np.concatenate([result.x, [0.0]])
        self.fitted = True
        return self

    def predict_proba(self, item_i, item_j):
        """Predict probability that item_i beats item_j."""
        if not self.fitted:
            raise ValueError("Model must be fitted before prediction")
        return sigmoid(self.beta[item_i] - self.beta[item_j])

    def get_rankings(self):
        """Return items sorted by strength (best first)."""
        return np.argsort(-self.beta)

This implementation uses vectorized operations for efficiency. The _neg_log_likelihood method uses np.log1p(np.exp(-diffs)) rather than np.log(1 + np.exp(-diffs)), which is numerically more stable for large negative differences where np.exp(-diffs) might overflow. This kind of numerical care matters at scale, where comparisons can span items with very different quality levels.

Let's test this implementation with simulated preference data where we know the true underlying parameters:

In[22]:
Code
# Simulate data: 5 items with true strengths
np.random.seed(42)
n_items = 5
true_beta = np.array([2.0, 1.0, 0.0, -0.5, -1.5])

# Generate 500 random comparisons
n_comparisons = 500
winners = []
losers = []

for _ in range(n_comparisons):
    # Pick two random items
    i, j = np.random.choice(n_items, 2, replace=False)

    # Determine winner based on true probabilities
    prob_i_wins = sigmoid(true_beta[i] - true_beta[j])
    if np.random.random() < prob_i_wins:
        winners.append(i)
        losers.append(j)
    else:
        winners.append(j)
        losers.append(i)

winners = np.array(winners)
losers = np.array(losers)
In[23]:
Code
# Fit the model
model = BradleyTerryModel(n_items)
model.fit(winners, losers)

# Compare estimated to true parameters (both anchored at last item)
true_beta_anchored = true_beta - true_beta[-1]
Out[24]:
Console
Parameter Recovery:
Item   True β     Estimated β  Error   
------------------------------------
0      3.500      3.517        +0.017  
1      2.500      2.239        -0.261  
2      1.500      1.560        +0.060  
3      1.000      0.560        -0.440  
4      0.000      0.000        +0.000

The estimated parameters closely recover the true values. This shows that the Bradley-Terry model can reliably learn latent strengths from noisy pairwise comparisons. The errors are small relative to the true parameter gaps, meaning the model correctly recovers the ranking and approximate win probabilities even with just 500 comparisons across 5 items.

Visualizing Parameter Recovery

Out[25]:
Visualization
Scatter plot comparing true and estimated Bradley-Terry parameters for five items, with points clustered tightly along a dashed diagonal reference line showing near-perfect parameter recovery.
True versus estimated Bradley-Terry log-strength parameters for simulated data with 500 comparisons across 5 items. The tight clustering along the diagonal confirms accurate parameter recovery. Small deviations from the diagonal reflect the statistical uncertainty inherent in any finite sample.
Out[26]:
Visualization
Log-log error plot showing mean absolute error of Bradley-Terry parameter estimates declining from roughly 0.65 at 50 comparisons to below 0.1 at 2000 comparisons, with error bars showing trial-to-trial variability.
Mean absolute error of parameter estimates as a function of sample size, showing how estimation accuracy improves with more comparison data. The downward trend on the log-log scale indicates power-law convergence typical of maximum likelihood estimators, with error approximately halving each time the sample size quadruples.

The sample size plot reveals an important practical consideration for LLM evaluation. At 50 comparisons, the mean absolute error in log-strength is above 0.3, which translates to meaningful uncertainty in predicted win probabilities. At 500 comparisons, the error drops below 0.1, giving much more reliable estimates. At 2000 comparisons, the model has high accuracy. For a real evaluation platform with dozens or hundreds of models, getting 200 or more comparisons per model pair requires substantial annotation effort, which is why platforms like Chatbot Arena aggregate crowd-sourced votes from many users over extended time periods.