Riddler Classic 2022-12-09

Riddler Riddler s/o Probability Simulation

2022 FIFA World Cup – Can You Win The Riddler Football Playoff?

Ryan McShane https://ryanmcshane.com
2022-12-09

The question from fivethirtyeight.com:

Speaking of “football,” the Riddler Football Playoff (RFP) consists of four teams. Each team is assigned a random real number between 0 and 1, representing the “quality” of the team. If team A has quality a and team B has quality b, then the probability that team A will defeat team B in a game is a/(a+b). And the probability that team B will defeat team A is b/(a+b). There are no ties!

In the semifinal games of the playoff, the team with the highest quality (the “1 seed”) plays the team with the lowest quality (the “4 seed”), while the other two teams play each other as well. The two teams that win their respective semifinal games then play each other in the final.

On average, what is the quality of the RFP champion?

Definitions

The strength ratings, before they are drawn and ordered, are distributed \(\text{Uniform}(0,1)\).

Let team \(A\) have rating \(a\), team \(B\) have rating \(b\), team \(C\) have rating \(c\), team \(D\) have rating \(d\), and \(a \ge b \ge c \ge d\) so that team \(A\) is ranked first and team \(D\) is ranked last. Thus, the first round will be \(A\) vs \(D\) and \(B\) vs \(C\). We assume that each game is an independent event, so that the outcome of one game does not affect any other game.

Let \(A \rightarrow D\) indicate that team \(A\) beat team \(D\).

We know that if \(X_1, X_2, X_3, X_4 \stackrel{iid}{\sim} \text{Uniform}(0,1)\) and if \(a = X_{(1)}, b = X_{(2)}, c = X_{(3)}, d = X_{(4)}\) are the order statistics, then each of \(a, b, c,\) and \(d\) have a \(\text{Beta}\) distribution. That is,

\[a \sim \text{Beta} (\alpha = 4, \beta = 1),\] \[b \sim \text{Beta} (\alpha = 3, \beta = 2),\] \[c \sim \text{Beta} (\alpha = 2, \beta = 3),\] \[d \sim \text{Beta} (\alpha = 1, \beta = 4),\]

where the general form of the Beta distribution is

\[f_Y(y) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}y^{\alpha - 1}(1-y)^{\beta - 1}\mathbb{I}_{[0<y<1]},\]

and the gamma function evaluates to \(\Gamma(n) = (n-1)!\) when \(n\) is an integer. We’ll need these distributions as priors for our later probability calculations. We also note that

\[E[Y] = \frac{\alpha}{\alpha + \beta} \;\; \text{ and } \;\; Var(Y) = \frac{\alpha\beta}{(\alpha + \beta)^2(\alpha + \beta + 1)}.\]

Simulation Approach

So, we can simulate \(M\) times using the sampling scheme described in the definitions – draw four \(\text{Uniform}(0,1)\), sort high to low, and assign to \(a, b, c,\) and \(d\).

M = 10000000
# M = 100
set.seed(538)
group_ratings = matrix(runif(4*M), nrow = M)
for(i in 1:nrow(group_ratings)) {
  group_ratings[i, ] = sort(group_ratings[i, ], decreasing = TRUE)
}
colnames(group_ratings) = letters[1:4]

A quick sanity check – we examine the first three simulations (to verify the randomly drawn values of \(a, b, c, d\) are \(\in [0,1]\) and ordered from best to worst).

\(a\) \(b\) \(c\) \(d\)
0.5946701 0.5004485 0.1986629 0.0034940
0.7773826 0.6322653 0.5646906 0.0867222
0.9377972 0.4526582 0.3499206 0.1221021

Perhaps more importantly, the expectations of these values match what they should be:

\(\hat E[a]\) \(\hat E[b]\) \(\hat E[c]\) \(\hat E[d]\)
0.8 0.6 0.4 0.2

And we can see the empirical densities match the Beta distributions (in yellow) we’ve described above.

Solution

We see that

\[P(A \rightarrow D, B \rightarrow C, A \rightarrow B \;|\; a, b, c, d) = \frac{a}{a + d}\frac{b}{b + c}\frac{a}{a + b} = \frac{a^2b}{(a+d)(b+c)(a+b)}.\]

Then, since \(A\) will always play \(D\) and \(B\) will always play \(C\), we know that every conditional probability will be of the form

\[\frac{1}{(a+d)(b+c)}\cdot\frac{x^2y}{x+y},\]

where \(x\) is the tournament winner’s rating and \(y\) is the runner-up’s strength. The conditional probability that team \(X\) will win is then

\[\frac{1}{(a+d)(b+c)}\left(\frac{x^2y}{x+y} + \frac{x^2z}{x+z}\right),\]

where \(z\) is the rating of the other runner-up’s strength. E.g.

\[P(A \text{ wins} \; | \; a, b, c, d) = \frac{1}{(a+d)(b+c)}\left(\frac{a^2b}{a+b} + \frac{a^2c}{a+c}\right).\]

## Function reflecting above derivation
prob_calculator = function(winner, a, b, c, d){
  s1 = ifelse(winner %in% c("a", "d"), "b", "a")
  s2 = ifelse(winner %in% c("a", "d"), "c", "d")
  s1 = case_when(s1 == "a" ~ a, s1 == "b" ~ b, s1 == "c" ~ c, s1 == "d" ~ d)
  s2 = case_when(s2 == "a" ~ a, s2 == "b" ~ b, s2 == "c" ~ c, s2 == "d" ~ d)
  winner = case_when(winner == "a" ~ a, winner == "b" ~ b, winner == "c" ~ c,
                     winner == "d" ~ d)
  pre = 1/((a + d)*(b + c))
  post = winner^2*(s1/(winner + s1) + s2/(winner + s2))
  return(pre*post)
}

Now we use this prob_calculator function on our simulated \(a, b, c,\) and \(d\) values from earlier.

exp_results = group_ratings %>%
  as.data.frame() %>%
  mutate(
    probA = prob_calculator(winner = "a", a = a, b = b, c = c, d = d),
    probB = prob_calculator(winner = "b", a = a, b = b, c = c, d = d), 
    probC = prob_calculator(winner = "c", a = a, b = b, c = c, d = d), 
    probD = prob_calculator(winner = "d", a = a, b = b, c = c, d = d)
    ) %>% 
  mutate(winner_rating = a*probA + b*probB + c*probC + d*probD, 
         check = probA + probB + probC + probD)

Results

For the purposes of the table, let \(P(A) := P(A \text{ wins }|\; a, b, c, d)\), and define \(\beta\) as the expected winning team’s rating when given \(a, b, c,\) and \(d\). That is,

\[ \begin{align} \beta = \;& aP(A \text{ wins }|\; a, b, c, d) + bP(B \text{ wins }|\; a, b, c, d) \;+ \\ & cP(C \text{ wins }|\; a, b, c, d) + dP(D \text{ wins }|\; a, b, c, d). \end{align} \]

Finally, \(\Sigma := \sum\limits_{i \in A, B, C, D} P(i \text{ wins }|\; a, b, c, d)\). As another check, we verify that this is \(1\) for every simulation. Now, we examine the first three simulations:

\(a\) \(b\) \(c\) \(d\) \(P(A)\) \(P(B)\) \(P(C)\) \(P(D)\) \(\beta\) \(\Sigma\)
0.5947 0.5004 0.1987 0.0035 0.5982 0.3294 0.0724 0.0001 0.5349 1
0.7774 0.6323 0.5647 0.0867 0.5079 0.2598 0.2196 0.0127 0.6842 1
0.9378 0.4527 0.3499 0.1221 0.6175 0.2136 0.1421 0.0268 0.7288 1 
## Indeed, all probability sums are 1
all.equal(exp_results$check, rep(1, M))
[1] TRUE

Now, we examine at the distribution of \(\beta\):

And finally, find the expected values of our simulations (remember that there were \(M = 10^7\) of them).

\(\hat E[a]\) \(\hat E[b]\) \(\hat E[c]\) \(\hat E[d]\) \(\hat P(A)\) \(\hat P(B)\) \(\hat P(C)\) \(\hat P(D)\) \(\hat E[\beta]\)
0.7999 0.6 0.4 0.2 0.5028 0.2867 0.1505 0.06 0.6735

Thus, \(E[\beta] = 0.6735\).

Update 12/18/2022

I got a shout-out on the Riddler blog for my second plot! Here it is.

© Copyright 2022 Ryan McShane, Ph.D.