Preference Evaluation › Code Implementation

import numpy as np
import pandas as pd

np.random.seed(42)

# True underlying model strengths (Bradley-Terry beta parameters)
true_strengths = {
    "Model_A": 1.2,
    "Model_B": 0.8,
    "Model_C": 0.0,
    "Model_D": -0.5,
}

models = list(true_strengths.keys())
n_models = len(models)


def simulate_match(model_i, model_j, n_games=1):
    """
    Simulate pairwise comparisons using Bradley-Terry model.
    Returns: (wins_i, wins_j, ties)
    """
    beta_i = true_strengths[model_i]
    beta_j = true_strengths[model_j]

    # Probability i beats j
    p_win = 1 / (1 + np.exp(beta_j - beta_i))
    # Small probability of tie (simplified model)
    p_tie = 0.05
    p_win = p_win * (1 - p_tie)
    p_loss = 1 - p_win - p_tie

    outcomes = np.random.choice(
        [1, 0, 0.5], size=n_games, p=[p_win, p_loss, p_tie]
    )
    wins_i = np.sum(outcomes == 1)
    wins_j = np.sum(outcomes == 0)
    ties = np.sum(outcomes == 0.5)
    return int(wins_i), int(wins_j), int(ties)


# Generate comparison matrix
comparison_matrix = pd.DataFrame(0, index=models, columns=models, dtype=float)
n_matches_per_pair = 100

for i, model_i in enumerate(models):
    for j, model_j in enumerate(models):
        if i < j:
            wins_i, wins_j, ties = simulate_match(
                model_i, model_j, n_matches_per_pair
            )
            comparison_matrix.loc[model_i, model_j] = wins_i + 0.5 * ties
            comparison_matrix.loc[model_j, model_i] = wins_j + 0.5 * ties

print(f"Generated {n_matches_per_pair} comparisons per pair")

class EloRatingSystem:
    def __init__(self, models, initial_rating=1500, k_factor=32):
        self.ratings = {model: initial_rating for model in models}
        self.k = k_factor
        self.history = {model: [initial_rating] for model in models}

    def expected_score(self, rating_a, rating_b):
        """Calculate expected score for player A against player B."""
        return 1 / (1 + 10 ** ((rating_b - rating_a) / 400))

    def update(self, model_a, model_b, score_a):
        """
        Update ratings based on outcome.
        score_a: 1 for win, 0.5 for tie, 0 for loss
        """
        expected_a = self.expected_score(
            self.ratings[model_a], self.ratings[model_b]
        )
        expected_b = 1 - expected_a

        # Update ratings
        self.ratings[model_a] += self.k * (score_a - expected_a)
        self.ratings[model_b] += self.k * ((1 - score_a) - expected_b)

        # Record history
        self.history[model_a].append(self.ratings[model_a])
        self.history[model_b].append(self.ratings[model_b])

    def fit_comparison_matrix(self, comparison_matrix, n_iterations=20):
        """
        Fit ratings given a comparison matrix using online updates.
        comparison_matrix[i,j] = expected wins of i over j
        """
        # Create list of individual matches for simulation
        matches = []
        for i, model_i in enumerate(comparison_matrix.index):
            for j, model_j in enumerate(comparison_matrix.columns):
                if i != j:
                    expected_wins = comparison_matrix.loc[model_i, model_j]
                    n_total = n_matches_per_pair
                    # Approximate: treat expected wins as win probability
                    n_wins = int(expected_wins)
                    n_losses = n_total - n_wins
                    matches.extend([(model_i, model_j, 1.0)] * n_wins)
                    matches.extend([(model_j, model_i, 1.0)] * n_losses)

        np.random.shuffle(matches)

        for model_a, model_b, score_a in matches:
            self.update(model_a, model_b, score_a)

        return self.ratings


# Initialize and fit Elo system
elo = EloRatingSystem(models, k_factor=16)
elo_ratings = elo.fit_comparison_matrix(comparison_matrix, n_iterations=10)

import numpy as np
from scipy.optimize import minimize


def bradley_terry_mle(comparison_matrix):
    """
    Fit Bradley-Terry model using maximum likelihood estimation.
    Returns beta parameters (relative strengths).
    """
    n = len(comparison_matrix)
    model_names = comparison_matrix.index.tolist()

    # Initial guess: all zeros
    beta_init = np.zeros(n)

    def negative_log_likelihood(beta):
        # Constraint: sum of beta = 0
        beta_adj = beta - np.mean(beta)
        nll = 0
        for i in range(n):
            for j in range(n):
                if i != j:
                    w_ij = comparison_matrix.iloc[i, j]
                    if w_ij > 0:
                        prob = 1 / (1 + np.exp(beta_adj[j] - beta_adj[i]))
                        # Avoid log(0)
                        prob = np.clip(prob, 1e-10, 1 - 1e-10)
                        nll -= w_ij * np.log(prob)
        return nll

    # Optimize
    result = minimize(negative_log_likelihood, beta_init, method="L-BFGS-B")
    beta_opt = result.x - np.mean(result.x)

    return dict(zip(model_names, beta_opt))


bt_strengths = bradley_terry_mle(comparison_matrix)

def beta_to_elo(beta, scale=400 / np.log(10)):
    """Convert Bradley-Terry beta to Elo scale (centered at 1500)."""
    return 1500 + beta * scale


def bootstrap_confidence_intervals(
    comparison_matrix, n_bootstrap=1000, alpha=0.05
):
    """
    Calculate confidence intervals for rating differences via bootstrapping.
    """
    model_names = comparison_matrix.index.tolist()
    n_models = len(model_names)
    bootstrap_ratings = []

    for _ in range(n_bootstrap):
        # Resample comparisons
        boot_matrix = pd.DataFrame(0.0, index=model_names, columns=model_names)

        for i in range(n_models):
            for j in range(n_models):
                if i < j:
                    total = int(
                        comparison_matrix.iloc[i, j]
                        + comparison_matrix.iloc[j, i]
                    )
                    if total > 0:
                        p_win = comparison_matrix.iloc[i, j] / total
                        # Resample wins
                        boot_wins = np.random.binomial(total, p_win)
                        boot_matrix.iloc[i, j] = boot_wins
                        boot_matrix.iloc[j, i] = total - boot_wins

        # Fit BT model
        try:
            boot_beta = bradley_terry_mle(boot_matrix)
            boot_elo = {m: beta_to_elo(b) for m, b in boot_beta.items()}
            bootstrap_ratings.append(boot_elo)
        except:
            continue

    # Calculate percentiles
    results = {}
    for model in model_names:
        ratings = [r[model] for r in bootstrap_ratings]
        results[model] = {
            "mean": np.mean(ratings),
            "lower": np.percentile(ratings, 100 * alpha / 2),
            "upper": np.percentile(ratings, 100 * (1 - alpha / 2)),
            "std": np.std(ratings),
        }

    return results


confidence_intervals = bootstrap_confidence_intervals(comparison_matrix)