The August 7, 2026 problems posed on the Fiddler on the Proof blog go as follows.
The Fiddler Baseball League consists of exactly two teams of equal skill: the Algebraists and the Geometers. Over the course of a season, these two teams play each other 162 times. Each team has an equal chance of winning each game, and the results of games are independent of one another.
At the end of the season, on average, how many games would you expect the team with the better record to have won? (If the teams have the same record, then you should include one of them in your calculation.)
If we let \(X\) be the number of wins of the winning team, then it is related to the binomial random variable \(Y\) with parameters \(n=162\) and \(p=1/2\) in that \(X = |Y - 81| + 81\).
So, for \(Y= 0, 1, 2, \ldots 80,81,82, \ldots 161,162\), we have \(X = 162, 161, \ldots, 82, 81, 82, \ldots, 161, 162\). Taking a sum-product of this sequence with the binomial distribution of \(Y\) should give us the answer.
X <- c(seq(162,81,-1),seq(82,162,1))
sum(X*dbinom(0:162, 162, 1/2))
## [1] 86.06988
Let’s do some simulations to see if we get a similar result.
simSeason <- function(N){
mean(abs(rbinom(N,162,1/2)-81)+81)
}
set.seed(42)
simSeason(1000000)
## [1] 86.0714
simSeason(1000000)
## [1] 86.0724
simSeason(1000000)
## [1] 86.06986
The simulations seem to confirm the answer of approximately 86.07.
After some expansion, the Fiddler Baseball League now boasts 30 teams. Over the course of a season, each team plays each other team five times. (Thus, each team plays a total of 145 games.) As before, each team has an equal chance of winning each game, and the results of games are independent of one another.
At the end of the season, on average, how many games would you expect the team with the best record to have won? (If more than one team has the same best record, then you should include one of them in your calculation.)
After messing around with this problem quite a bit, I’m pretty sure a nice closed form is too computationally difficult to achieve, but I hope I’m proved wrong.
Let’s first look at some simulations.
## This function is to simulate one season and output the max number of wins
## by a team.
season <- function(teams = 30){
S <- matrix(NA, nrow=teams, ncol = teams)
for(i in 1:teams){
for(j in 1:i){
if(j < i){
S[i,j] <- 5-S[j,i]
}
if(j==i){
S[j,j] <- 0
if(j<30){
S[j,(j+1):teams] <- rbinom(teams-j, 5, 1/2)
}
}
}
}
max(rowSums(S))
}
## Here a few season outputs.
season()
## [1] 83
season()
## [1] 83
season()
## [1] 86
Let’s get a few results from 100,000 simulations.
seasonSim <- function(N, teams=30){
mean(sapply(rep(teams,N), season))
}
seasonSim(100000, 30)
## [1] 84.97459
seasonSim(100000, 30)
## [1] 84.97824
seasonSim(100000, 30)
## [1] 84.9849
We know the win number of each of the thirty teams individually, \(W_i\), can be modeled by itself as a binomial random variable with parameters \(n=145\) and \(p=1/2\), which has mean \(\mu=72.5\) and variance \(\sigma^2 = 36.25\).
There is a small covariance between two of these variables. The 28 other games are independent, but the 5 games they play against each other are dependent. So, \(Cov(W_i,W_j) = Cov(X_i, 5-X_i)\). This is just \(0-Var(X_i)\) or -1.25.
The multivariate distribution \(W \sim \scr{N}_{30}(\mathbf{\mu}, \mathbf{\Sigma})\), where \(\Sigma\) has 36.25 down the diagonal and -1.25 everywhere else. Row sums are zero.
Using 30 independent standard normal random variables \(Z_i\) and their mean \(\overline{Z}\), note that \(\text{Var}(\overline{Z}) = 1/30\). We can approximate the wins now with \(W_i \approx 72.5 + A(Z_i - \overline{Z})\) for an appropriate choice of \(A\).
We know \(\text{Var}(W_i) = 36.25\) and \(\text{Var}(Z_i-\overline{Z}) = 1 - 1/30 = 29/30\). Thus, \(\text{Var}(W_i) = 36.25 = A^2*\frac{29}{30}\) which implies \(A = \sqrt{37.5}\). This also ties out the covariance, as \(\text{Cov}(Z_i - \overline{Z}, Z_j - \overline{Z}) = -\frac{1}{30}\), and therefore \(\text{Cov}(W_i,W_j) = 37.5 \left(-\frac{1}{30}\right) = -1.25\).
Now, looking at the max number of wins, we can estimate this by looking at \(\max_i(Z_i - \overline{Z})\) and taking an expected value. Since \(E(\overline{Z})=0\), we just want the expected value of the max over 30 standard normals, call that \(M_{30}\).
The distribution function for \(M_{30}\) is \(\Phi(z)^{30}\), and so the density of \(M_{30}\) is \(30\phi(z)\Phi(z)^{29}\). Let’s multiply by \(z\) and integrate to see what that expected max is.
M30_Ez <- function(z){
30*z*dnorm(z)*pnorm(z)^29
}
(maxZ <- integrate(M30_Ez, -Inf, Inf)$value)
## [1] 2.042761
Now, using that 2.042761, let’s find an estimated maximum.
72.5+ sqrt(37.5)*maxZ
## [1] 85.0093
The estimated maximum is a little above 85. Simulations have it a little lower than this. From the simulations and the math estimation, I would say 84.98 is as good a guess as I can make.