From the FiveThirtyEight website:
“Riddler League Baseball, also known as the RLB, consists of three teams: the Mississippi Moonwalkers, the Delaware Doubloons and the Tennessee Taters.
Each time a batter for the Moonwalkers comes to the plate, they have a 40 percent chance of getting a walk and a 60 percent chance of striking out. Each batter for the Doubloons, meanwhile, hits a double 20 percent percent of the time, driving in any teammates who are on base, and strikes out the remaining 80 percent of the time. Finally, each batter for the Taters has a 10 percent chance of hitting a home run and a 90 percent chance of striking out.
During the RLB season, each team plays an equal number of games against each opponent. Games are nine innings long and can go into extra innings just like in other baseball leagues. Which of the three teams is most likely to have the best record at the end of the season?"
When I saw this problem, my first impulse was to create a simulation, but as I thought about how to do that, I realized that it’s possible to calculate an exact answer (up to a small error caused by cutting the tails of infinite distributions). This is my solution, including all the R code needed to reproduce it. (Scroll to the bottom of this page to see my answer.)
The first step is to produce for each team the probability distribution for the number of runs they will score in one inning. The key tool here is the negative binomial distribution, which, given a number of failures \(r\) and a probability \(p\) of failure, tells us the probability of having \(k\) successes before the \(r\)th failure. This is given by R’s dnbinom function. For example, the Taters’ probability of hitting 2 home runs before their third strikeout is
dnbinom(2,3,0.9)## [1] 0.04374Since every home run results in a run, it is simple to calculate the Taters’ single inning run distribution, which we store in a data frame.
taters.one.inning <- data.frame(runs = integer(21),prob = numeric(21))
taters.one.inning$runs <- 0:20
taters.one.inning$prob <- dnbinom(taters.one.inning$runs,3,0.9)
taters.one.inning[0:11,]## runs prob
## 1 0 7.2900e-01
## 2 1 2.1870e-01
## 3 2 4.3740e-02
## 4 3 7.2900e-03
## 5 4 1.0935e-03
## 6 5 1.5309e-04
## 7 6 2.0412e-05
## 8 7 2.6244e-06
## 9 8 3.2805e-07
## 10 9 4.0095e-08
## 11 10 4.8114e-09The case of the Doubloons is slightly more complicated, because the first double does not score a run. So we first calculate the distribution for the number of doubles and from that calcuate the distribution of runs scored.
library(dplyr)
doubloons.one.inning <- data.frame(runs = integer(21),prob = numeric(21))
doubloons.one.inning$doubles <- 0:20
doubloons.one.inning$prob <- dnbinom(doubloons.one.inning$doubles,3,0.8)
doubloons.one.inning$runs <- ifelse(doubloons.one.inning$doubles <= 1,0,doubloons.one.inning$doubles-1)
doubloons.one.inning <- doubloons.one.inning %>%
group_by(runs) %>%
summarise(prob = sum(prob)) %>%
data.frame()
doubloons.one.inning[0:11,]## runs prob
## 1 0 8.192000e-01
## 2 1 1.228800e-01
## 3 2 4.096000e-02
## 4 3 1.228800e-02
## 5 4 3.440640e-03
## 6 5 9.175040e-04
## 7 6 2.359296e-04
## 8 7 5.898240e-05
## 9 8 1.441792e-05
## 10 9 3.460301e-06
## 11 10 8.178893e-07The case of the Moonwalkers is similar, but they don’t score until the fourth walk of the inning.
moonwalkers.one.inning <- data.frame(runs = integer(21),prob = numeric(21))
moonwalkers.one.inning$walks <- 0:20
moonwalkers.one.inning$prob <- dnbinom(moonwalkers.one.inning$walks,3,0.6)
moonwalkers.one.inning$runs <- ifelse(moonwalkers.one.inning$walks<=3,0,moonwalkers.one.inning$walks-3)
moonwalkers.one.inning <- moonwalkers.one.inning %>%
group_by(runs) %>%
summarise(prob = sum(prob)) %>%
data.frame()
moonwalkers.one.inning[0:11,]## runs prob
## 1 0 0.8208000000
## 2 1 0.0829440000
## 3 2 0.0464486400
## 4 3 0.0247726080
## 5 4 0.0127401984
## 6 5 0.0063700992
## 7 6 0.0031142707
## 8 7 0.0014948499
## 9 8 0.0007066563
## 10 9 0.0003297730
## 11 10 0.0001522029Before we start comparing teams, we have to calculate for each team the probability distribution of the total number of runs they score in nine innings of play. We assume the innings are independent and use the following function to calculate the distribution of the sum of two independent random variables.
sum.of.dists <- function(d1,d2){
temp <- merge(d1,d2,by=NULL)
temp$runs <- temp$runs.x + temp$runs.y
temp$prob <- temp$prob.x*temp$prob.y
temp <- temp %>%
group_by(runs) %>%
summarise(prob = sum(prob)) %>%
data.frame()
return(temp)
}We apply this function to each one-inning distribution eight times.
moonwalkers.nine.innings <- moonwalkers.one.inning
doubloons.nine.innings <- doubloons.one.inning
taters.nine.innings <- taters.one.inning
for(i in 1:8){
moonwalkers.nine.innings <- sum.of.dists(moonwalkers.nine.innings,moonwalkers.one.inning)
doubloons.nine.innings <- sum.of.dists(doubloons.nine.innings,doubloons.one.inning)
taters.nine.innings <- sum.of.dists(taters.nine.innings,taters.one.inning)
}Here are the distributions together with the average runs per game for each team. (The latter is not relevant to the analysis but it will make for an interesting observation when we’re done.)
moonwalkers.nine.innings[0:11,]## runs prob
## 1 0 0.16909709
## 2 1 0.15378935
## 3 2 0.14828531
## 4 3 0.13021207
## 5 4 0.10709491
## 6 5 0.08382566
## 7 6 0.06307747
## 8 7 0.04594866
## 9 8 0.03256496
## 10 9 0.02253967
## 11 10 0.01528054sum(moonwalkers.nine.innings$runs*moonwalkers.nine.innings$prob)## [1] 3.363757doubloons.nine.innings[0:11,]## runs prob
## 1 0 0.166153499
## 2 1 0.224307224
## 3 2 0.209353409
## 4 3 0.159258129
## 5 4 0.105854318
## 6 5 0.063737832
## 7 6 0.035534532
## 8 7 0.018610429
## 9 8 0.009250349
## 10 9 0.004396885
## 11 10 0.002010206sum(doubloons.nine.innings$runs*doubloons.nine.innings$prob)## [1] 2.358taters.nine.innings[0:11,]## runs prob
## 1 0 0.058149737
## 2 1 0.157004290
## 3 2 0.219806006
## 4 3 0.212479139
## 5 4 0.159359354
## 6 5 0.098802800
## 7 6 0.052694826
## 8 7 0.024841847
## 9 8 0.010557785
## 10 9 0.004105805
## 11 10 0.001478090sum(taters.nine.innings$runs*taters.nine.innings$prob)## [1] 3We are finally ready to do some head-to-head comparisons, which will be aided by the following function:
win.probs <- function(d1,d2){
temp <- merge(d1,d2,by=NULL)
temp$prob <- temp$prob.x*temp$prob.y
d1.win.prob <- sum(temp$prob[temp$runs.x > temp$runs.y])
d2.win.prob <- sum(temp$prob[temp$runs.x < temp$runs.y])
tie.prob <- sum(temp$prob[temp$runs.x == temp$runs.y])
return(c(d1.win.prob,d2.win.prob,tie.prob))
}
temp <- win.probs(moonwalkers.nine.innings,doubloons.nine.innings)
mw.beat.dbl.9in <- temp[1]
dbl.beat.mw.9in <- temp[2]
mw.dbl.tie.9in <- temp[3]
temp <- win.probs(moonwalkers.nine.innings,taters.nine.innings)
mw.beat.ttr.9in <- temp[1]
ttr.beat.mw.9in <- temp[2]
mw.ttr.tie.9in <- temp[3]
temp <- win.probs(doubloons.nine.innings,taters.nine.innings)
dbl.beat.ttr.9in<- temp[1]
ttr.beat.dbl.9in <- temp[2]
dbl.ttr.tie.9in <- temp[3]Of course, a baseball game cannot end in a tie. When a game is tied after nine innings, the teams continue to play one inning at a time until one team scores more runs in a single inning. Thus we need to apply our win.probs function to the one-inning run distributions we calculated initially. We also use the fact that a team’s probability of winning after the game goes into extra innings is its probability of scoring more runs in one inning conditional on the inning not being tied.
temp <- win.probs(moonwalkers.one.inning,doubloons.one.inning)
mw.beat.dbl.extra.in <- temp[1]/(temp[1]+temp[2])
dbl.beat.mw.extra.in <- temp[2]/(temp[1]+temp[2])
temp <- win.probs(moonwalkers.one.inning,taters.one.inning)
mw.beat.ttr.extra.in <- temp[1]/(temp[1]+temp[2])
ttr.beat.mw.extra.in <- temp[2]/(temp[1]+temp[2])
temp <- win.probs(doubloons.one.inning,taters.one.inning)
dbl.beat.ttr.extra.in <- temp[1]/(temp[1]+temp[2])
ttr.beat.dbl.extra.in <- temp[2]/(temp[1]+temp[2])Now we can calculate the final head-to-head win probabilities and expected season win percentages.
mw.beat.dbl <- mw.beat.dbl.9in + mw.dbl.tie.9in*mw.beat.dbl.extra.in
mw.beat.dbl## [1] 0.5848533dbl.beat.mw <- dbl.beat.mw.9in + mw.dbl.tie.9in*dbl.beat.mw.extra.in
dbl.beat.mw## [1] 0.4151428mw.beat.ttr <- mw.beat.ttr.9in + mw.ttr.tie.9in*mw.beat.ttr.extra.in
mw.beat.ttr## [1] 0.4838547ttr.beat.mw <- ttr.beat.mw.9in + mw.ttr.tie.9in*ttr.beat.mw.extra.in
ttr.beat.mw## [1] 0.5161414dbl.beat.ttr <- dbl.beat.ttr.9in + dbl.ttr.tie.9in*dbl.beat.ttr.extra.in
dbl.beat.ttr## [1] 0.3727171ttr.beat.dbl <- ttr.beat.dbl.9in + dbl.ttr.tie.9in*ttr.beat.dbl.extra.in
ttr.beat.dbl## [1] 0.6272829Since the teams play an equal number of games against each opponent, a team’s season winning percentage is the average of its win probabilites against each opponent.
0.5*(mw.beat.dbl+mw.beat.ttr)## [1] 0.5343540.5*(dbl.beat.mw+dbl.beat.ttr)## [1] 0.393930.5*(ttr.beat.mw+ttr.beat.dbl)## [1] 0.5717122The Taters are most likely to have the best record for the season. They get shut out much less often (less than 6% of their games versus more than 16% for the others), and this outweighs the fact that the Moonwalkers, on average, score more runs per 9 innings.