Estimating the parameters

The most common application of the Bradley–Terry model is to infer the values of the parameters ${\displaystyle p_{i}}$ given an observed set of outcomes ${\displaystyle i>j}$, such as wins and losses in a competition. The simplest way to estimate the parameters is by maximum likelihood estimation, i.e., by maximizing the likelihood of the observed outcomes given the model and parameter values.

Suppose we know the outcomes of a set of pairwise competitions between a certain group of individuals, and let w_ij be the number of times individual i beats individual j. Then the likelihood of this set of outcomes within the Bradley–Terry model is ${\displaystyle \prod _{ij}[P(i>j)]^{w_{ij}}}$ and the log-likelihood of the parameter vector p = [p₁, ..., p_n] is^[1]

${\displaystyle {\begin{aligned}L(\mathbf {p} )&=\ln \prod _{ij}[P(i>j)]^{w_{ij}}=\sum _{i=1}^{n}\sum _{j=1}^{n}\ln \left({\frac {p_{i}}{p_{i}+p_{j}}}\right)^{w_{ij}}\\[6pt]&=\sum _{ij}w_{ij}\ln {\frac {p_{i}}{p_{i}+p_{j}}}=\sum _{ij}[w_{ij}\ln p_{i}-w_{ij}\ln(p_{i}+p_{j})].\end{aligned}}}$

Zermelo^[5] showed that this expression has only a single maximum, which can be found by differentiating with respect to ${\displaystyle p_{i}}$ and setting the result to zero, which leads to

${\displaystyle p_{i}={\frac {\sum _{j}w_{ij}}{\sum _{j}(w_{ij}+w_{ji})/(p_{i}+p_{j})}}.}$

(2)

This equation has no known closed-form solution, but Zermelo suggested solving it by simple iteration. Starting from any convenient set of (positive) initial values for the ${\displaystyle p_{i}}$, one iteratively performs the update

${\displaystyle p_{i}'={\frac {\sum _{j}w_{ij}}{\sum _{j}(w_{ij}+w_{ji})/(p_{i}+p_{j})}}}$

(3)

for all i in turn. The resulting parameters are arbitrary up to an overall multiplicative constant, so after computing all of the new values they should be normalized by dividing by their geometric mean thus:

${\displaystyle p_{i}\leftarrow {\frac {p'_{i}}{\left(\prod _{j=1}^{n}p'_{j}\right)^{1/n}}}.}$

(4)

This estimation procedure improves the log-likelihood on every iteration, and is guaranteed to eventually reach the unique maximum.^[5][10] It is, however, slow to converge.^[1][11] More recently it has been pointed out^[12] that equation (2) can also be rearranged as

${\displaystyle p_{i}={\frac {\sum _{j}w_{ij}p_{j}/(p_{i}+p_{j})}{\sum _{j}w_{ji}/(p_{i}+p_{j})}},}$

which can be solved by iterating

${\displaystyle p_{i}'={\frac {\sum _{j}w_{ij}p_{j}/(p_{i}+p_{j})}{\sum _{j}w_{ji}/(p_{i}+p_{j})}},}$

(5)

again normalizing after every round of updates using equation (4). This iteration gives identical results to the one in (3) but converges much faster and hence is normally preferred over (3).^[12]

Worked example of solution procedure

Consider a sporting competition between four teams, who play a total of 22 games among themselves. Each team's wins are given in the rows of the table below and the opponents are given as the columns:

Results

	A	B	C	D
A	–	2	0	1
B	3	–	5	0
C	0	3	–	1
D	4	0	3	–

For example, Team A has beat Team B twice and lost to team B three times; not played team C at all; won once and lost four times against team D.

We would like to estimate the relative strengths of the teams, which we do by calculating the parameters ${\displaystyle p_{i}}$, with higher parameters indicating greater prowess. To do this, we initialize the four entries in the parameter vector p arbitrarily, for example assigning the value 1 to each team: [1, 1, 1, 1]. Then we apply equation (5) to update ${\displaystyle p_{1}}$, which gives

$${\displaystyle p_{1}={\frac {\sum _{j(\neq 1)}w_{1j}p_{j}/(p_{1}+p_{j})}{\sum _{j(\neq 1)}w_{j1}/(p_{1}+p_{j})}}={\frac {2{\frac {1}{1+1}}+0{\frac {1}{1+1}}+1{\frac {1}{1+1}}}{3{\frac {1}{1+1}}+0{\frac {1}{1+1}}+4{\frac {1}{1+1}}}}=0.429.}$$

Now, we apply (5) again to update ${\displaystyle p_{2}}$, making sure to use the new value of ${\displaystyle p_{1}}$ that we just calculated:

$${\displaystyle p_{2}={\frac {\sum _{j(\neq 2)}w_{2j}p_{j}/(p_{2}+p_{j})}{\sum _{j(\neq 2)}w_{j2}/(p_{2}+p_{j})}}={\frac {3{\frac {0.429}{1+0.429}}+5{\frac {1}{1+1}}+0{\frac {1}{1+1}}}{2{\frac {1}{1+0.429}}+3{\frac {1}{1+1}}+0{\frac {1}{1+1}}}}=1.172}$$

Similarly for ${\displaystyle p_{3}}$ and ${\displaystyle p_{4}}$ we get

$${\displaystyle p_{3}={\frac {\sum _{j(\neq 3)}w_{3j}p_{j}/(p_{3}+p_{j})}{\sum _{j(\neq 3)}w_{j3}/(p_{3}+p_{j})}}={\frac {0{\frac {0.429}{1+0.429}}+3{\frac {1.172}{1+1.172}}+1{\frac {1}{1+1}}}{0{\frac {1}{1+0.429}}+5{\frac {1}{1+1.172}}+3{\frac {1}{1+1}}}}=0.557}$$

$${\displaystyle p_{4}={\frac {\sum _{j(\neq 4)}w_{4j}p_{j}/(p_{4}+p_{j})}{\sum _{j(\neq 4)}w_{j4}/(p_{4}+p_{j})}}={\frac {4{\frac {0.429}{1+0.429}}+0{\frac {1.172}{1+1.172}}+3{\frac {0.557}{1+0.557}}}{1{\frac {1}{1+0.429}}+0{\frac {1}{1+1.172}}+1{\frac {1}{1+0.557}}}}=1.694}$$

Then we normalize all the parameters by dividing by their geometric mean ${\displaystyle (0.429\times 1.172\times 0.557\times 1.694)^{1/4}=0.830}$ to get the estimated parameters p = [0.516, 1.413, 0.672, 2.041].

To improve the estimates further, we repeat the process, using the new p values. For example,

$${\displaystyle p_{1}={\frac {2\cdot {\frac {1.413}{0.516+1.413}}+0\cdot {\frac {0.672}{0.516+0.672}}+1\cdot {\frac {2.041}{0.516+2.041}}}{3\cdot {\frac {1}{0.516+1.413}}+0\cdot {\frac {1}{0.516+0.672}}+4\cdot {\frac {1}{0.516+2.041}}}}=0.725.}$$

Repeating this process for the remaining parameters and normalizing, we get p = [0.677, 1.034, 0.624, 2.287]. Repeating a further 10 times gives rapid convergence toward a final solution of p = [0.640, 1.043, 0.660, 2.270]. This indicates that Team D is the strongest and Team B the second strongest, while Teams A and C are nearly equal in strength but below Teams B and D. In this way the Bradley–Terry model lets us infer the relationship between all four teams, even though not all teams have played each other.