An implementation of the Bradley-Terry model for paired comparisons.

bt.py

import numpy as np

def p(i, w, p_prime):
    def summand(w, i, j, p_prime):
        return (w[i][j] + w[j][i])/(p_prime[i] + p_prime[j])
    def denom(i, w, p_prime):
        return sum(summand(w, i, j, p_prime) for j in range(len(w.T)) if j != i) 
    return np.nansum(w[i])/denom(i, w, p_prime)

def bt(w, n):
    p_prime = np.ones(len(w))
    for iteration in range(n):
        ps = [p(i, w, p_prime) for i in range(len(w))]
        p_prime = [p_i / sum(ps) for p_i in ps]
    return p_prime

def P(ps, i, j):
    return ps[i] / (ps[i] + ps[j])

# Tests
import pytest
from itertools import combinations

@pytest.fixture
def data(): 
    return np.array([[np.nan, 2, 0, 1], 
                     [3, np.nan, 5, 0], 
                     [0, 3, np.nan, 1],
                     [4, 0, 3, np.nan]])

def test_p(data):
    p_prime = [1,1,1,1]
    assert p(0, data, p_prime) == 0.6
    assert p(1, data, p_prime) == pytest.approx(1.231, abs=1e-3)

def test_bt(data):
    assert bt(data, 20) == pytest.approx([0.138, 0.226, 0.143, 0.491], abs=1e-3)

def test_P(data):
    ps = bt(data, 20)
    assert P(ps, 0, 1) == pytest.approx(.380, abs=1e-3)
    for i,j in combinations(range(4), 2):
        assert P(ps, i, j) + P(ps, j, i) == 1