spb/llmindex Public

1# llmindex.io — Bradley-Terry fit (MM algorithm, ties as half-wins)2# Author:  Simon-Pierre Boucher3# Contact: contact@spboucher.ai4# License: Proprietary — © Simon-Pierre Boucher, all rights reserved56from __future__ import annotations78from dataclasses import dataclass910import numpy as np111213@dataclass14class BTResult:15    log_strength: np.ndarray  # (M,) log-strengths, mean 016    log_strength_se: np.ndarray  # (M,) SEs from the Fisher information17    iterations: int18    converged: bool192021def fit_bradley_terry(22    wins: np.ndarray,23    max_iterations: int = 1000,24    tolerance: float = 1e-10,25    damping: float = 0.1,26) -> BTResult:27    """Fit Bradley-Terry strengths from a wins matrix.2829    wins[i, j] = (possibly fractional) number of wins of i over j.30    Ties should be pre-encoded as 0.5 win to each side.31    `damping` adds a tiny uniform prior so isolated models stay finite.32    """33    W = np.asarray(wins, dtype=float)34    if W.ndim != 2 or W.shape[0] != W.shape[1]:35        raise ValueError("wins must be a square matrix")36    M = W.shape[0]37    # Regularization: everyone gets `damping` phantom wins vs everyone else.38    W = W + damping * (np.ones((M, M)) - np.eye(M))39    N = W + W.T  # total comparisons between each pair40    p = np.ones(M)41    converged = False42    it = 043    for it in range(1, max_iterations + 1):44        denom = (N / (p[:, None] + p[None, :] + 1e-300)).sum(1) - np.diag(45            N / (p[:, None] + p[None, :] + 1e-300)46        )47        new_p = W.sum(1) / np.maximum(denom, 1e-300)48        new_p = new_p / np.exp(np.log(new_p + 1e-300).mean())  # geometric-mean normalize49        if np.max(np.abs(new_p - p)) < tolerance:50            p = new_p51            converged = True52            break53        p = new_p54    log_strength = np.log(p + 1e-300)55    log_strength -= log_strength.mean()56    # Fisher information of log-strengths: I_ii = sum_j N_ij * p_ij * (1 - p_ij)57    # where p_ij = p_i / (p_i + p_j); SE_i = 1 / sqrt(I_ii).58    P = p[:, None] / (p[:, None] + p[None, :] + 1e-300)59    info = (N * P * (1 - P)).sum(1) - np.diag(N * P * (1 - P))60    se = 1.0 / np.sqrt(np.maximum(info, 1e-9))61    return BTResult(log_strength=log_strength, log_strength_se=se, iterations=it, converged=converged)62