probability - Book Proofs

Can you hop to the lily pad?

This week’s Fiddler is about hopping back and forth.

You are a frog in a pond with an infinite number of lily pads in a line, marked “1,” “2,” “3,” etc. You are currently on pad 2, and your goal is to make it to pad 1. From any given pad, there are specific probabilities that you’ll jump to another pad: Whenever you are on pad $k$, you will hop to pad $k−1$ with probability $1/k$, and you will hop to pad $k+1$ with probability $(k−1)/k$.

What is the probability that you will ultimately make it to pad 1?

My solution:
[Show Solution]

The general finite case

We will first solve the most general form of this problem: arbitrary probabilities and an arbitrary starting point! Suppose there are $n$ pads numbered $1,2,\dots,n$. Let $a_k$ be the probability that we will make it to pad 1 starting from pad $k$. Further assume that when we’re at pad $k$, we transition to $k-1$ with probability $p_k$ and to $k+1$ with probability $q_k$. We have $p_k+q_k=1$ so these are the only two possibilities. This leads to the recurrence relation:
\begin{align*}
a_1 &= 1 \\
a_k &= p_k a_{k-1} + q_k a_{k+1}\qquad\text{for }k=2,3,\dots,n-1\\
a_n &= 0
\end{align*} Start by rewriting the recurrence relation using the fact that $p_k+q_k=1$:
\begin{align*}
a_1 &= 1 \\
p_k (a_{k}-a_{k-1}) &= q_k(a_{k+1}-a_k) \qquad \text{for }k=2,3,\dots,n-1\\
a_n &= 0
\end{align*}Iterating the recurrence starting from $k-1$, we obtain
\begin{align*}
(a_k-a_{k-1}) &= \frac{p_{k-1}}{q_{k-1}}(a_{k-1}-a_{k-2}) \\
&= \frac{p_{k-1}}{q_{k-1}} \cdot \frac{p_{k-2}}{q_{k-2}} (a_{k-2}-a_{k-3}) \\
&\;\;\vdots \\
&= (a_2-a_1)\prod_{i=2}^{k-1} \frac{p_{i}}{q_{i}}
\end{align*}Now sum from $k=2$ to $k=n$ and the left-hand side will telescope.
\[
a_n-a_1 = \sum_{k=2}^n (a_k-a_{k-1}) = (a_2-a_1)\sum_{k=2}^n\prod_{i=2}^{k-1} \frac{p_{i}}{q_{i}}
\]Letting $a_1=1$ and $a_n=0$, we can solve for $a_2$ and obtain:
\[
a_2 = 1-\frac{1}{\sum_{k=2}^n \prod_{i=2}^{k-1} \frac{p_i}{q_i}}
\]Now instead of summing up to $k=n$ in the telescoping sum, sum to $k=m$ to find a general term, and obtain:
\[
a_m-a_1 = (a_2-a_1)\sum_{k=2}^m\prod_{i=2}^{k-1} \frac{p_{i}}{q_{i}}
\]Substituting the values we found for $a_1$ and $a_2$, we obtain the following general formula for any $a_m$:

$\displaystyle
a_m = 1-\frac{\sum_{k=2}^m \prod_{i=2}^{k-1} \frac{p_i}{q_i}}{\sum_{k=2}^n \prod_{i=2}^{k-1} \frac{p_i}{q_i}},\qquad \text{for }m=1,2,\dots,n
$

The formula even works at the boundaries. When $m=1$, the numerator is zero (empty sum), and we recover $a_1=1$. When $m=n$, the numerator matches the denominator and we recover $a_n=0$.

This problem can also be viewed as that of finding the stationary distribution of a Markov chain. However, this approach leads to a formula that is more difficult to simplify involving the inverse of an $n\times n$ matrix, which is why I opted for the approach above.

The general infinite case

We now consider the case where there are infinitely many lily pads. The recurrence relation now looks like:
\begin{align*}
a_1 &= 1 \\
a_k &= p_k a_{k-1} + q_k a_{k+1}\qquad\text{for }k=2,3,\dots
\end{align*}Such recurrence relations are also called one-dimensional random walks.

This case is a bit trickier than the finite case. We are dealing with a second-order difference equation, so two boundary conditions are required to determine a unique solution. This was no problem in the finite case we previously solved, because we had the boundary conditions $a_1=1$ and $a_n=0$. But in this infinite case, we only have the condition $a_1=1$. This means there will be infinitely many possible solutions, each with a different value at infinity; the value of $a_\infty = \lim_{k\to\infty}a_k$. In order to have a unique solution, we must use the correct “boundary condition at infinity”.

To illustrate this fact, observe that we can satisfy the recurrence relation by setting $a_k=1$ for all $k$. This scenario has the boundary condition $a_\infty=1$ and is known as gambler’s ruin (the pads are viewed as amounts of money, and each turn we bet and either gain more money or lose money). It’s called “gambler’s ruin” because we always eventually go broke (return to pad 1). One example of gambler’s ruin is the case $p_k=q_k=\frac{1}{2}$; i.e., we flip a fair coin at every pad to see if we move forward or backward.

In cases where $q_k \gt p_k$ (we are likelier to move to a larger pad than to a smaller one), there is a non-zero probability that we will never return to pad 1. This cases is characterized by the property that the farther we get from pad 1, the less likely we are to return. In other words, our boundary condition at infinity should be $a_\infty = 0$.

To solve this case, we can use the exact same approach as in the finite case, except instead of summing up to $k=n$ and using $a_n=0$ to find $a_2$, we sum up to $k=\infty$ and use $a_\infty=0$ to find $a_2$. This results in the similar-looking formula:

$\displaystyle
a_m = 1-\frac{\sum_{k=2}^m \prod_{i=2}^{k-1} \frac{p_i}{q_i}}{\sum_{k=2}^{\infty} \prod_{i=2}^{k-1} \frac{p_i}{q_i}},\qquad \text{for }m=1,2,\dots
$

Note: To obtain the set of all possible solutions to the recurrence relation, we need only take affine combinations of the two fundamental solutions (the all-ones gambler’s ruin solution and the solution above). In other words, the general solution is given by:
\[
\begin{bmatrix}
a_1^\text{gen} \\ a_2^\text{gen} \\ \vdots
\end{bmatrix}
=
\alpha
\begin{bmatrix}
1 \\ 1 \\ \vdots
\end{bmatrix}
+ (1-\alpha)
\begin{bmatrix}
a_1 \\ a_2 \\ \vdots
\end{bmatrix}
\]where $a_1,a_2,\dots$ is the solution we derived above. This general solution satisfies $\lim_{k\to\infty}a_k = \alpha$, so based on how we choose $\alpha\in\mathbb{R}$, we can achieve any desired boundary condition at infinity.

Solving our special case

In our special case of interest, we have $p_i=\frac{1}{i}$ and $q_i=\frac{i-1}{i}$. Substituting into our formula, we obtain:
\begin{align*}
a_m &= 1-\frac{\sum_{k=2}^m \prod_{i=2}^{k-1} \frac{1}{i-1}}{\sum_{k=2}^n \prod_{i=2}^{k-1} \frac{1}{i-1}}
= 1-\frac{\sum_{k=2}^m \frac{1}{(k-2)!}}{\sum_{k=2}^n \frac{1}{(k-2)!}}
= 1-\frac{\sum_{k=0}^{m-2} \frac{1}{k!}}{\sum_{k=0}^{n-2} \frac{1}{k!}}
\end{align*}We can also write this as a single fraction:
\[
a_m = \frac{\frac{1}{(m-1)!} + \cdots + \frac{1}{(n-2)!}}{\frac{1}{0!}+\frac{1}{1!}+\frac{1}{2!}+\cdots+\frac{1}{(n-2)!}}\qquad\text{for }m=1,2,\dots,n
\]

Example 1: If we have $n=4$ pads and we start on pad $m=2$, then the probability we eventually end up on pad 1 is given by:
\[
a_2 = \frac{\frac{1}{1!}+\frac{1}{2!}}{\frac{1}{0!}+\frac{1}{1!}+\frac{1}{2!}}
=\frac{3}{5} = 60\%.
\]

Example 2: If we have $n=\infty$, we recognize the sum in the denominator as the infinite series representation for Euler’s constant $e$, and the sum in the numerator is one less. Therefore,
\[
a_2 = \frac{\frac{1}{1!}+\frac{1}{2!}+\cdots}{\frac{1}{0!}+\frac{1}{1!}+\frac{1}{2!}+\cdots}
=\frac{e-1}{e} \approx 63.212\%.
\]

If instead of starting on pad 2 we start on a larger-numbered pad, the numerator will continue to lose terms in its sum and gradually degrade to zero the farther we start.

Halloween Puzzle

This week’s Fiddler is about rounding!

You are presented with a bag of treats, which contains $n \geq 3$ peanut butter cups and some unknown quantity of candy corn kernels (with any amount being equally likely). You reach into the bag $k$ times, with $3 \leq k \leq n$, and pull out a candy at random. Each time, it’s a peanut butter cup! How many candy kernels do you expect to be in the bag?

My solution:
[Show Solution]

Define the following random variables:

$m$: number of candy corn kernels in the bag
$K$: whether we draw $k$ peanut butter cups in a row (true/false)

We are asked to calculate the expected number of corn kernels in the bag given that we draw $k$ peanut butter cups in a row. In other words, we want to find $\mathbb{E}(m\mid K)$.

An easier quantity to calculate is the probability that we draw $k$ peanut butter cups in a row given that there are $m$ candy corn kernels in the bag. This is given by:
\begin{align*}
\mathbb{P}(K\mid m)
&= \frac{n}{n+m}\cdot \frac{n-1}{n+m-1}\cdots\frac{n-k+1}{n+m-k+1} \\
&= \frac{n!(n+m-k)!}{(n-k)!(n+m)!} \\
&= \binom{n}{k}\binom{n+m}{k}^{-1}
\end{align*}Using Bayes’ rule, we can express the desired expectation in terms of the conditional probability above:
\begin{align*}
\mathbb{E}(m\mid K)
&= \sum_{m=0}^\infty m\cdot \mathbb{P}(m\mid K)
= \sum_{m=0}^\infty m\cdot \frac{\mathbb{P}(K \mid m)\mathbb{P}(m)}{\mathbb{P}(K)}
\end{align*} Now substitute: $\mathbb{P}(K) = \sum_{m=0}^\infty \mathbb{P}(K\mid m)\mathbb{P}(m)$ and obtain:
\begin{align*}
\mathbb{E}(m\mid K)
&= \frac{\sum_{m=0}^\infty m\cdot \mathbb{P}(K \mid m)\mathbb{P}(m)}{\sum_{m=0}^\infty \mathbb{P}(K \mid m)\mathbb{P}(m)} \\
&= \frac{\sum_{m=0}^\infty m\cdot \mathbb{P}(K \mid m)}{\sum_{m=0}^\infty \mathbb{P}(K \mid m)}
\end{align*}where in the last step, we used the fact that $\mathbb{P}(m)$ is the same for all $m$ so we canceled it from the numerator and denominator. Substituting our expression for $\mathbb{P}(K\mid m)$, we obtain:
\begin{align*}
\mathbb{E}(m\mid K)
&= \frac{\sum_{m=0}^\infty m\binom{n}{k}\binom{n+m}{k}^{-1}}{\sum_{m=0}^\infty \binom{n}{k}\binom{n+m}{k}^{-1}}
= \frac{\sum_{m=0}^\infty m\binom{n+m}{k}^{-1}}{\sum_{m=0}^\infty \binom{n+m}{k}^{-1}}
%= \frac{n-k+1}{k-2}
\end{align*}
To simplify this, define the following quantity:
\[
f(k,a,b) := \sum_{i=a}^b \binom{i}{k}^{-1}
\]and observe that we have:
\[
\mathbb{E}(m\mid K) = \frac{\sum_{j=1}^\infty f(k,n+j,\infty)}{f(k,n,\infty)}
\]It turns out we can evaluate $f$, and that it has the nice closed-form expression:
\[
f(k,a,b) = \frac{k}{k-1}\left( \binom{a-1}{k-1}^{-1}-\binom{b}{k-1}^{-1}\right)
\]This formula can be proved by induction and is related to something called the German tank problem.

With this formula in hand, we can take the limit $b\to\infty$ and see that:
\[
f(k,a,\infty) = \frac{k}{k-1}\binom{a-1}{k-1}^{-1}
\]Substituting into our expression for the desired expectation, we obtain:
\begin{align*}
\mathbb{E}(m\mid K)
&= \frac{\sum_{j=1}^\infty f(k,n+i,\infty)}{f(k,n,\infty)} \\
&= \frac{\sum_{j=1}^\infty \frac{k}{k-1}\binom{n+j-1}{k-1}^{-1}}{\frac{k}{k-1}\binom{n-1}{k-1}^{-1}} \\
&= \binom{n-1}{k-1}\sum_{j=1}^\infty \binom{n+j-1}{k-1}^{-1} \\
&= \binom{n-1}{k-1} f(k-1,n,\infty) \\
&= \binom{n-1}{k-1} \frac{k-1}{k-2} \binom{n-1}{k-2}^{-1} \\
&= \frac{(n-1)!}{(n-k)!(k-1)!} \frac{k-1}{k-2} \frac{(n-k+1)!(k-2)!}{(n-1)!} \\
&= \frac{n-k+1}{k-2}
\end{align*}

Therefore, our final answer is:

$\displaystyle
\mathbb{E}(m\mid K) = \frac{n-k+1}{k-2}
$

I undoubtedly found the most complicated possible solution for this problem… So if somebody can show me a more elegant solution (perhaps a counting argument?) then I would be much obliged!

Round, round, get a round

This week’s Fiddler is about rounding!

Let $\text{round}(x)$ be the value of $x$ rounded to the nearest integer. Suppose $x_1,\dots,x_n$ are independent uniformly distributed random variables in $[0,1]$. Find the probability that
\[
\text{round}(x_1+\cdots+x_n) = \text{round}(x_1)+\cdots+\text{round}(x_n)
\]

My solution:
[Show Solution]

Let’s call the probability we seek $p(n)$. The values of the $x_i$ determine what sums are even possible, so let’s consider some cases based on the possible values of $\text{round}(x_1)+\cdots+\text{round}(x_n)$. Suppose that $k$ of the variables are in in the interval $[\tfrac{1}{2},1]$ and the remaining $n-k$ variables are in $[0,\tfrac{1}{2}]$. The probability of this occurring is $\tfrac{1}{2^n}\binom{n}{k}$. In this case,
\[
\sum_{i=1}^n \text{round}(x_i) = k
\]Note: We can ignore the issue of whether $\tfrac{1}{2}$ rounds to $0$ or $1$ since the probability of $\tfrac{1}{2}$ occurring is zero. Define a set of rescaled variables $y_i$ in $[0,1]$ as follows:
\[
y_i = \begin{cases}
2x_i-1 & \text{if }i=1,\dots,k\\
2x_i & \text{if }i=k+1,\dots,n
\end{cases}
\]Now let’s calculate the probability that the sum of the $x_i$ also rounds to $k$ and express it in terms of the $y_i$:
\begin{align*}
\mathrm{Prob}\left( \text{round}\biggl(\sum_{i=1}^n x_i\biggr) = k\right)
&= \mathrm{Prob}\left( k-\tfrac{1}{2} \leq \sum_{i=1}^n x_i \leq k+\tfrac{1}{2} \right) \\
&= \mathrm{Prob}\left( k-1 \leq \sum_{i=1}^n y_i \leq k+1 \right)
\end{align*}The sum of the $y_i$, which is a sum of $n$ independent variables in $[0,1]$ follows a so-called Irwin-Hall distribution, whose CDF is given by:
\[
F(x) = \text{Prob}\bigl( y_1+\cdots+y_n \leq x \bigr) = \frac{1}{n!}\sum_{k=0}^{\lfloor x \rfloor} (-1)^k \binom{n}{k}(x-k)^n
\]Therefore, the probability we seek is $F(k+1)-F(k-1)$, or:
\[
\frac{1}{n!}\sum_{m=0}^{k+1} (-1)^m \binom{n}{m}(k+1-m)^n
-\frac{1}{n!}\sum_{m=0}^{k-1} (-1)^m \binom{n}{m}(k-1-m)^n
\]Rearranging terms a bit, we obtain:
\[
\frac{(-1)^{n+k}}{n!} \binom{n}{k}+
\frac{1}{n!}\sum_{m=0}^{k} (-1)^m \binom{n}{m}\bigl( (k+1-m)^n-(k-1-m)^n\bigr)
\]By writing out the terms carefully, we can group like $n^\text{th}$ powers and we obtain a simpler equivalent expression:
\[
\frac{1}{n!}\sum_{m=0}^k(-1)^m\left[ \binom{n+1}{m}-\binom{n+1}{m-1}\right](k+1-m)^n
\]Remember this was the probability that the two roundings are the same when $k$ of the $x_i$’s round to $1$. Multiplying each such probability by its prior probability and summing, we get the total probability that the two roundings are the same:
\[
p(n) = \frac{1}{2^n\,n!}\sum_{k=0}^n \binom{n}{k}\sum_{m=0}^k(-1)^m\left[ \binom{n+1}{m}-\binom{n+1}{m-1}\right](k+1-m)^n
\]This can be further simplified by grouping $n^\text{th}$ powers once again and observing that all odd terms cancel out. After quite a bit of algebra, we obtain the simplest expression possible, which is:

$\displaystyle
p(n) = \frac{1}{2^n\, n!}\sum _{k=0}^{\left\lfloor \frac{n}{2}\right\rfloor } (-1)^k \binom{n+1}{k} (n+1-2k)^n
$

We can evaluate this for different values of $n$ and we obtain:
\[
\begin{array}{c|cccccccccc}
n & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10\\ \hline
p(n) & 1 & \frac{3}{4} & \frac{2}{3} & \frac{115}{192} & \frac{11}{20} & \frac{5887}{11520} & \frac{151}{315} & \frac{259723}{573440} & \frac{15619}{36288} & \frac{381773117}{928972800}
\end{array}
\]

The sequence $2^n \, n!\, p(n)$ is actually a well-known sequence. It is A261398 in OEIS (the Online Encyclopedia of Integer Sequences), and according to OEIS it has no known closed-form formula and the one above is the simplest one available. Good enough for me!

Visualization in 2D and beyond

We can visualize the probabilities for the case $n=2$ by coloring in the set of points $(x_1,x_2)$ for which the sum of the rounded values equals the rounded value of the sum. This yields:

We can check the shaded region is $\frac{3}{4}$ of the total area, as anticipated. We can do the same for $n=3$. This time, we plot the set of points $(x_1,x_2,x_3)$ and the probability is the shaded volume:

Again, we can check the volume is $\frac{2}{3}$. Unfortunately, we can’t visualize the case $n=4$, but we can try… This should be a 4D hypercube. In the same way that each slice of a 3D cube is a square, each slice of a 4D hypercube is a 3D cube. Here is an animation showing what happens as we vary $x_4 \in [0,1]$ and we plot the cube resulting from the valid choices of $(x_1,x_2,x_3)$ as a function of $x_4$. The abrupt change occurs when $x_4=\tfrac{1}{2}$ as this causes the sum of rounded values to jump by $1$.

So now, if you can imagine, the probability we’re looking for is the “hyper-volume” of this 4D shape, which is given by integrating the volume above as we vary $x_4$, and as it turns out,
\[
p(4) = \int_{0}^1 \text{Volume}(x_4)\,\mathrm{d}x_4 = \frac{115}{192}.
\]

Approximation

Recall that the solution we found above is given by
\[
p(n) = \sum_{k=0}^n \frac{1}{2^n}\binom{n}{k}\bigl( F(k+1)-F(k-1) \bigr)
\]where $F(x)$ is the CDF of the Irwin-Hall distribution. The formula we found for this is quite messy, so let’s see if we can approximate it.

First, the Irwin-Hall distribution can be well-approximated by a normal distribution. This makes sense, since summing a large number of identically distributed random variables tends to a normal distribution by the central limit theorem. In this case, the limiting distribution has mean $\mu=\frac{n}{2}$ and variance $\sigma^2=\frac{n}{12}$. Therefore, we can approximate:
\begin{align*}
F(k+1)-F(k-1) &\approx 2F'(k) \\
&\approx \frac{2}{\sqrt{2\pi \sigma^2}} \exp\biggl( -\frac{(x-\mu)^2}{2\sigma^2} \biggr) \\
&= \frac{1}{\sqrt{\pi n/24}}\exp\biggl( -\frac{(k-\tfrac{n}{2})^2}{n/6} \biggr)
\end{align*}Next, the weighted binomial sum is actually a binomial distribution, which we can also approximate using a normal distribution. In general, we have $B(n,p) \approx \mathcal{N}(\mu,\sigma^2)$, with $\mu=np$ and $\sigma^2=np(1-p)$. Therefore, if $g$ is any function, we can write the binomial sum as an expectation:
\begin{align*}
\sum_{k=0}^n \frac{1}{2^n}\binom{n}{k} g(k)
&= \mathbb{E}_{k \sim B(n,\tfrac{1}{2})}\bigl( g(k) \bigr) \\
&\approx \mathbb{E}_{k \sim \mathcal{N}(\tfrac{n}{2},\tfrac{n}{4})}\bigl( g(k) \bigr) \\
&=\frac{1}{\sqrt{\pi n/2}} \int_{-\infty}^\infty \exp\biggl(-\frac{(k-\tfrac{n}{2})^2}{n/2} \biggr)g(k)\,\mathrm{d}k
\end{align*}Setting $g(k)=F(k+1)-F(k-1)$ and substituting our previous approximation, we find:
\begin{align*}
p(n) &\approx \frac{1}{\sqrt{\pi n/2}}\frac{1}{\sqrt{\pi n/24}}\int_{-\infty}^\infty \exp\biggl( -\frac{(k-\tfrac{n}{2})^2}{n/6} \biggr)\exp\biggl(-\frac{(x-\tfrac{n}{2})^2}{n/2} \biggr)\,\mathrm{d}k \\
&= \frac{\sqrt{48}}{\pi n} \int_{-\infty}^\infty \exp\biggl(-\frac{(k-\tfrac{n}{2})^2}{n/8} \biggr) \,\mathrm{d}k \\
&= \frac{\sqrt{48}}{\pi n} \sqrt{\pi n/8} \\
\end{align*}Therefore, we have:

$\displaystyle
p(n) \approx \sqrt{\frac{6}{\pi n}}
$

To check our approximation, we can plot it together with the true values:

Not bad! So in conclusion, we have an exact formula and an asymptotically exact approximation for the probability that
\[
\text{round}(x_1+\cdots+x_n) = \text{round}(x_1)+\cdots+\text{round}(x_n)
\]They are given by:

$\displaystyle
p(n) = \frac{1}{2^n\, n!}\sum _{k=0}^{\left\lfloor \frac{n}{2}\right\rfloor } (-1)^k \binom{n+1}{k} (n+1-2k)^n
\approx \sqrt{\frac{6}{\pi n}}
$

Showcase Showdown

This week’s Fiddler is based on “Showcase Showdown” on the game show “The Price is Right”.

Suppose we have some number of players. Player A is the first to spin a giant wheel, which spits out a real number chosen randomly and uniformly between 0 and 1. All spins are independent of each other. After spinning, A can either stick with the number they just got or spin the wheel one more time. If they spin again, their assigned number is the sum of the two spins, as long as that sum is less than or equal to 1. If the sum exceeds 1, A is immediately declared a loser.

After A is done spinning (whether once or twice), B steps up to the wheel. Like A, they can choose to spin once or twice. If they spin twice and the sum exceeds 1, they are similarly declared the loser. This continues until all players are done. Whoever has the greater value (that does not exceed 1) is declared the winner.

Assuming all players play the game optimally, what are player A’s chances of winning?

My solution:
[Show Solution]

We’ll assume there are $n$ players in total. Define the following:
\[
W_k(y) = \left\{ \begin{aligned}
&\text{Probability that the $k^\text{th}$}\\
&\text{player from the end will win}\\
&\text{assuming the highest number}\\
&\text{obtained so far is $y$.}
\end{aligned}
\right\}
\]In other words, $W_1$ is the probability the last player wins, $W_2$ is the probability the second-last player wins, etc. This is the quantity we need to solve for. It will also be helpful to define the probability that the last $k$ players all lose:
\[
L_k(y) = \left\{ \begin{aligned}
&\text{Probability that the last $k$}\\
&\text{players will be declared}\\
&\text{losers assuming the highest}\\
&\text{number obtained so far is $y$.}
\end{aligned}
\right\}
\]Let’s work our way backwards, starting from the last player. Suppose the number to beat is $y$. Since there is only one player left, $W_1(y) + L_1(y) = 1$. Let’s say they obtain $x$ on their first spin. Clearly, if $x > y$, they will stop (and win the game). Otherwise, they will try their chances with a second spin. Say they obtain $z$. They will win if they obtain a sum $x+z$ that is larger than $y$ but not larger than $1$. Therefore,
\begin{align*}
W_1(y) &= \int_y^1\,\mathrm{d}x + \int_0^y \int_{y-x}^{1-x}\, \mathrm{d}z\, \mathrm{d}x = 1-y^2
\end{align*}Consequently, $L_1(y)=y^2$. This makes sense. If $y$ (the number to beat) approaches $1$, then the probability of winning approaches zero because it becomes increasingly difficult to beat it.

Now consider the case with $k$ players remaining. Again, the number to beat is $y$. We will assume a threshold strategy: the player will spin a second time if their first spin is less than $a$, where $a$ is the threshold value. Clearly, if their first spin is less than $y$, they must spin again, otherwise they are guaranteed to lose. So $y < a < 1$. There are several cases to consider. Suppose the first spin is $x$ and the second spin, if needed, is $z$.

If $y < a < x$, then $W_k(y) = L_{k-1}(x)$.
If $y < x < a$, then if $x+z > 1$, we lose. Else, $W_k(y) = L_{k-1}(x+z)$.
If $x < y < a$, then if $x+z > 1$ or $x+z < y$, we lose. Otherwise, $W_k(y) = L_{k-1}(x+z)$.

Expressing as integrals, the probability of winning with threshold $a$ is:
\begin{align*}
W_k(y,a) &=
\int_a^1 L_{k-1}(x)\,\mathrm{d}x
+ \int_y^a \int_{0}^{1-x} L_{k-1}(x+z)\,\mathrm{d}z\, \mathrm{d}x
+ \int_0^y \int_{y-x}^{1-x} L_{k-1}(x+z)\,\mathrm{d}z\, \mathrm{d}x\\
&= \int_a^1 L_{k-1}(x)\,\mathrm{d}x
+ \int_y^a \int_{x}^{1} L_{k-1}(v)\,\mathrm{d}v\, \mathrm{d}x
+ \int_0^y \int_{y}^{1} L_{k-1}(v)\,\mathrm{d}v\, \mathrm{d}x\\
&= \int_a^1 L_{k-1}(x)\,\mathrm{d}x
+ \int_y^a \int_{x}^{1} L_{k-1}(v)\,\mathrm{d}v\, \mathrm{d}x
+ y \int_{y}^{1} L_{k-1}(v)\,\mathrm{d}v
\end{align*}This player will pick $a \in [y,1]$ to maximize this probability of winning. In other words:
\[
W_k(y) = \underset{a \in [y,1]}{\text{maximize}}\; W_k(y,a)
\]The maximum will occur when the derivative with respect to $a$ is zero, or at $a=y$. Therefore, we have:
\[
a_k^\star(y) = \begin{cases}
\text{soln. of }L_{k-1}(a) = \int_{a}^1 L_{k-1}(v)\,\mathrm{d}v & \text{if it satisfies }a > y \\
y & \text{otherwise}
\end{cases}
\]Next, we can recursively compute $L_k$. In order for all players to lose, the current player must lose, and so must all others. This happens when the first spin is below the threshold, the second spin doesn’t beat $y$, and all subsequent players lose. In other words,
\begin{align*}
L_k(y) &= \int_y^{a_k^\star(y)}
\int_{1-x}^1 L_{k-1}(y)\,\mathrm{d}z\,\mathrm{d}x
+
\int_0^{y}\left(
\int_0^{y-x} L_{k-1}(y)\,\mathrm{d}z +
\int_{1-x}^1 L_{k-1}(y)\,\mathrm{d}z\right)\mathrm{d}x \\
&= \left( \int_y^{a_k^\star(y)}
\int_{1-x}^1 \mathrm{d}z\,\mathrm{d}x
+
\int_0^{y}\left(
\int_0^{y-x} \mathrm{d}z +
\int_{1-x}^1 \mathrm{d}z\right)\mathrm{d}x\right)L_{k-1}(y) \\
&= \left( \int_y^{a_k^\star(y)}
x\,\mathrm{d}x
+
\int_0^{y} y\, \mathrm{d}x\right)L_{k-1}(y) \\
&= \tfrac{1}{2}\left(a_k^\star(y)^2 + y^2\right) L_{k-1}(y)
\end{align*}The reason this can simplify so much is that the integrands are evaluated at $y$ because when the present player loses, the “best spin” does not change for subsequent players!

Two players

Applying the recursive formulas above, we find:
\begin{align*}
W_1(y) &= 1-y^2 \\
L_1(y) &= y^2 \\
W_2(y,a) &= \tfrac{1}{12}\left(4+4 a-4 a^3-a^4-3 y^4\right)
\end{align*}Things are starting to get messy… It turns out the optimal $a$ is given by:
\[
a_2^\star(y) = \max\left\{y, 2 \cos \left(\tfrac{2 \pi }{9}\right)-1\right\} \approx \max\{y,0.532089\}
\]Finding $W_2(y)$ is even messier, but thankfully, in this case, we can just let $y=0$ because the second last player is also the first player! This gives us a probability of winning of approximately:
\[
W_2(0) = \tfrac{1}{12} \left(9+4 \sin (\tfrac{\pi }{18})+2 \cos (\tfrac{\pi }{9})-8 \cos (\tfrac{2 \pi }{9})\right) \approx 0.453802
\]So the probability that the first (of two) players wins is about 45.38%. Here is a plot of $W_2(0,a)$ that shows how the probability of the first player winning varies as a function of their threshold. The function is maximized at $a^\star_2(0)$

Three players

Continuing our computations, we have:
\begin{align*}
a_2^\star(y) &\approx \max\{y,0.532089\} \\
L_2(y) &= \tfrac{1}{2}\left( y^2 + a_2^\star(y)^2 \right)y^2 \\
%W_3(y,a) &= -\tfrac{1}{24}\left(a^4+4 a^3-4 a+3 y^4-4\right) a_2^\star(y)^2\\
%&\hspace{3cm}-\tfrac{1}{60}\left(a^6+6 a^5-6 a+5 y^6-6\right)
W_3(y,a) &= \text{(messy)}
\end{align*} At this point, I abandoned any hope of analytic solutions… Optimizing numerically, I obtained:
\begin{align*}
a_3^\star(y) &\approx \max\{y,0.648655\} \\
W_3(0) &\approx 0.305227
\end{align*}So the probability that the first (of three) players wins is about 30.52%.

Similar to the previous section, here is a plot of $W_3(0,a)$ that shows how the probability of the first player winning varies as a function of their threshold. It’s a lot harder to win when there are three players!

More players

In principle, we can keep turning the crank and generating solutions for more players. The nice thing about this problem is that the optimal strategy for a given player does not depend on the total number of players; it only depends on how many players there are after the current player.

Here is a table of the thresholds I was able to calculate:

$k$	$a_k^\star(y)$ (threshold)	$W_k(0)$ prob of win
$1$	$\max\{y,0\}$	$1.000000$
$2$	$\max\{y,0.532089\}$	$0.453802$
$3$	$\max\{y,0.648655\}$	$0.305227$
$4$	$\max\{y,0.711449\}$	$0.231060$
$5$	$\max\{y,0.752249\}$	$0.186229$

Here is how you read the table: Suppose you are the $k^\text{th}$ player from the end, say $k=3$. When your turn comes, let $y$ be the largest value obtained by any player before you. Pick whichever is largest, $y$ or $0.648655$, this becomes your threshold. Take your first spin. If you obtain something larger than the threshold, you’re currently winning, so pass your turn. Otherwise, spin again and hope for the best. The last column shows the probability of this player winning assuming they go first.

As $k$ gets larger, so does the threshold; the more players will spin after you, the more aggressive you have to be with your thresholding strategy. Also, notice that the probability of winning when you go first is a little less than $1/k$; the first player is always at a disadvantage, so we can expect their probability of winning to be below average.

Dungeon Master’s Dice

This week’s Fiddler is a probability question about a dice-rolling game.

Two people are sitting at a table together, each with their own bag of six “DnD dice”: a d4, a d6, a d8, a d10, a d12, and a d20. Here, “dX” refers to a die with X faces, numbered from 1 to X, each with an equally likely probability of being rolled. Both people randomly pick one die from their respective bags and then roll them at the same time. For example, suppose the two dice selected are a d4 and a d12. The players roll them, and let’s further suppose that both rolls come up as 3. What luck! What’s the probability of something like this happening? That is, what is the probability that both players roll the same number, whether or not they happened to pick the same kind of die?

Extra Credit
Instead of two people sitting at the table, now suppose there are three. Again, all three randomly pick one die from their respective bags and roll them at the same time. For example, suppose the three dice selected are a d4, a d20, and a d12. The players roll them, and let’s further suppose that the d4 comes out as 4, the d20 comes out as 13, and the d12 comes out as 4. In this case, there are two distinct numbers (4 and 13) among the three rolls. On average, how many distinct numbers would you expect to see among the three rolls?

My solution:
[Show Solution]

Each person performs two actions: selecting a die at random, and rolling that die. In the language of probability, the sample space (the set of all things that can occur) is $(n,d)$, where $d$ is the selected die, and $n$ is the number that is rolled. The probability of picking the die $d$ and rolling a $n$ is given by:
\[
p_{n,d} = \frac{1}{6} \times \begin{cases}
\frac{1}{d} & \text{if }1 \leq n \leq d \\
0 & \text{otherwise}
\end{cases}
\]The $\frac{1}{6}$ comes from the fact that there are six different dice, each equally likely to be picked. Arranging these probabilities in a table, we obtain:
\[
p_{n,d} = \frac{1}{6}\begin{bmatrix}
\frac{1}{4} & \frac{1}{6} & \frac{1}{8} & \frac{1}{10} & \frac{1}{12} & \frac{1}{20} \\
\frac{1}{4} & \frac{1}{6} & \frac{1}{8} & \frac{1}{10} & \frac{1}{12} & \frac{1}{20} \\
\frac{1}{4} & \frac{1}{6} & \frac{1}{8} & \frac{1}{10} & \frac{1}{12} & \frac{1}{20} \\
\frac{1}{4} & \frac{1}{6} & \frac{1}{8} & \frac{1}{10} & \frac{1}{12} & \frac{1}{20} \\
0 & \frac{1}{6} & \frac{1}{8} & \frac{1}{10} & \frac{1}{12} & \frac{1}{20} \\
0 & \frac{1}{6} & \frac{1}{8} & \frac{1}{10} & \frac{1}{12} & \frac{1}{20} \\
0 & 0 & \frac{1}{8} & \frac{1}{10} & \frac{1}{12} & \frac{1}{20} \\
0 & 0 & \frac{1}{8} & \frac{1}{10} & \frac{1}{12} & \frac{1}{20} \\
0 & 0 & 0 & \frac{1}{10} & \frac{1}{12} & \frac{1}{20} \\
0 & 0 & 0 & \frac{1}{10} & \frac{1}{12} & \frac{1}{20} \\
0 & 0 & 0 & 0 & \frac{1}{12} & \frac{1}{20} \\
0 & 0 & 0 & 0 & \frac{1}{12} & \frac{1}{20} \\
0 & 0 & 0 & 0 & 0 & \frac{1}{20} \\
0 & 0 & 0 & 0 & 0 & \frac{1}{20} \\
0 & 0 & 0 & 0 & 0 & \frac{1}{20} \\
0 & 0 & 0 & 0 & 0 & \frac{1}{20} \\
0 & 0 & 0 & 0 & 0 & \frac{1}{20} \\
0 & 0 & 0 & 0 & 0 & \frac{1}{20} \\
0 & 0 & 0 & 0 & 0 & \frac{1}{20} \\
0 & 0 & 0 & 0 & 0 & \frac{1}{20}
\end{bmatrix}
\quad \implies \quad p_n = \begin{bmatrix}
\frac{31}{240} \\ \frac{31}{240} \\ \frac{31}{240} \\ \frac{31}{240} \\
\frac{7}{80} \\ \frac{7}{80} \\ \frac{43}{720} \\ \frac{43}{720} \\
\frac{7}{180} \\ \frac{7}{180} \\
\frac{1}{45} \\ \frac{1}{45} \\
\frac{1}{120} \\ \frac{1}{120} \\ \frac{1}{120} \\ \frac{1}{120} \\ \frac{1}{120} \\ \frac{1}{120} \\ \frac{1}{120} \\ \frac{1}{120}
\end{bmatrix}
\]In the matrix on the left, each column corresponds to a different die and each row corresponds to a different number. All the numbers in this table sum to 1 (this is a joint distribution). The column vector on the right is what you get when you sum across rows. These numbers again sum to 1 and form what is called a marginal distribution. This is the probability distribution for the dice rolls regardless of which die was picked. Picking one of the dice at random, then rolling it, is the same as rolling a single weighted 20-sided die whose weights are given by the marginal distribution above. We can plot this distribution as a graph, and we get the following:

From now on, I will write $p_n$ to denote the (marginal) probability of rolling the number $n$.

The first problem

The first problem asks for the probability that both players roll the same number (regardless of which die they chose). The probability that both players roll $n$ is $p_n^2$. Therefore, the probability that they both roll the same number is
\[
\sum_{n=1}^{20} p_n^2 = \frac{3}{32} = 0.09375
\]Therefore, the probability that both players roll the same number is approximately 9%.

The extra credit

We now have 3 players, and we are asked to find the expected number of distinct rolls. We can calculate the probability of each in turn.

The probability of obtaining one distinct number (all three players roll the same number) is similar to the case with two players rolling the same number, except now it’s $p_n^3$. Therefore,
\[
q_1 = \sum_{n=1}^{20}p_n^3 = \frac{32753}{3110400} \approx 0.01053
\]
The probability of obtaining two distinct numbers is a bit more complicated. We have to count how many ways we can pick two players that will have the same number $n$, which is $\binom{3}{2}p_n^2$, and then the third player can have any of the remaining numbers, so we multiply by $(1-p_n)$. Therefore,
\[
q_2 = \sum_{n=1}^{20}\binom{3}{2}p_n^2(1-p_n) = \frac{258847}{1036800} \approx 0.24966
\]
Finally, the probability of obtaining three distinct numbers is the leftover probability after excluding the first two cases above:
\[
q_3 = 1-q_1-q_2 = \frac{1150553}{1555200} \approx 0.73981
\]

We can calculate the expected number of distinct rolls by multiplying each probability above by the number of distinct rolls and summing them up:
\begin{align}
\mathbb{E}(\text{distinct numbers})
&= q_1 + 2q_2 + 3q_3
= \frac{8489153}{3110400}
\approx 2.72928
\end{align}

Therefore, we would expect to see, on average, 2.729 distinct numbers among the three dice rolls.

Don’t flip out

This week’s Fiddler is a probability question about a coin-flipping game.

Kyle and Julien are playing a game in which they each toss their own fair coins. On each turn of the game, both players flip their own coin once. If, at any point, Kyle’s most recent three flips are Tails, Tails, and Heads (i.e., TTH), then he wins. If, at any point, Julien’s most recent three flips are Tails, Tails, and Tails (i.e, TTT), then he wins.

However, both players can’t win at the same time. If Kyle gets TTH at the same time Julien gets TTT, then no one wins, and they continue flipping. They don’t start over completely or erase their history, mind you—they merely continue flipping, so that one of them could conceivably win in the next flip or two.

What is the probability that Kyle wins this game?

Extra Credit
Kyle and Julien write down all eight possible sequences for three coin flips (HHH, HHT, HTH, THH, HTT, THT, TTH, and TTT) on eight different slips of paper. They place these slips into a hat and shake it.

They will each randomly draw slips of paper out of the hat, at which point they will play the same game as previously described, but looking for the sequence specified on the slip of paper they each selected. Kyle draws first and looks at his slip of paper. After doing some calculations, he says: “Well, at this point, it’s about as fair a match as it could possible be.”

Which slip or slips of paper might Kyle have drawn? And what are his chances of winning at this point (i.e., before Julien selects his own slip of paper)?

My solution:
[Show Solution]

The state of each player’s game is determined by their most recent 3 coin flips. So there are 8 possible states:
\[
\textsf{HHH}, \textsf{HHT}, \textsf{HTH}, \textsf{HTT}, \textsf{THH}, \textsf{THT}, \textsf{TTH}, \textsf{TTT}
\]Each turn, each player transitions from one of these states to another, depending on whether they flip Heads or Tails. Therefore, it is a Markov Chain. Here is what it looks like represented graphically:

Each edge has a weight of $\frac{1}{2}$, since we are flipping a fair coin. Mathematically, we can write $P_{ij}$ as the probability of transitioning from state $i$ to state $j$, where the states $1,2,\dots,8$ are ordered in the same order as the list above. (${\small\textsf{HHH}}:1$, ${\small\textsf{HHT}}:2$, etc.) This leads to the following transition matrix:
\[
P = \frac{1}{2}\begin{bmatrix}
1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 \\
1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 1 & 1
\end{bmatrix}
\]Since both Kyle and Julien are flipping coins independently, what we really have is a Markov chain whose state is the Cartesian product of the states of each player: $(x_\text{Kyle},x_\text{Julien})$. This larger Markov chain has $8\times 8 = 64$ states and its transition matrix is the Kronecker product of two smaller transition matrices: $\hat P = P\otimes P$. The matrix $\hat P$ has the following sparsity pattern

All nonzero entries of $\hat P$ are $\frac{1}{4}$. This is a valid transition matrix because each row sums to 1 since there are exactly four nonzero entries per row. If you’re curious, here is what the graph looks like (yikes!)

In the augmented Markov chain, Kyle wins if we end up in state $(7,j)$ where $j\neq 8$. This corresponds to states numbered:
\[
\mathcal{X}_\text{K} = \{49, 50, 51, 52, 53, 54, 55\}
\]Likewise, Julien wins if we end up in state $(i,8)$ where $i\neq 7$, which corresponds to states:
\[
\mathcal{X}_\text{J} = \{8, 16, 24, 32, 40, 48, 64\}
\]As an additional notation, let’s define the set of all transient nodes, which are the nodes that are not in $\mathcal{X}_\text{K}$ or in $\mathcal{X}_\text{J}$. These are the 50 nodes:
\begin{multline}
\mathcal{X}_\text{T} =
\{1,2,3,4,5,6,7,9,10,11,12,13,14,15,17,18,19,20, \\
21,22,23,25,26,27,28,29,30,31,33,34,35,36,37,38, \\
39,41,42,43,44,45,46,47,56,57,58,59,60,61,62,63\}
\end{multline}Note that the case where both players win at the same time corresponds to $(7,8)$ (state number 56), which is treated as just another transient node, since nothing special happens in this case.

Now that we have the appropriate mathematical framework set up, we can tackle the original problem. We would like to find the probability that Kyle wins the game. Translated to Markov-speak, we want the probability that we first land in $\mathcal{X}_\text{K}$ before first landing in $\mathcal{X}_\text{J}$. Let $v_i$ be the probability of Kyle winning if we start in node $i$. We can make the following observations:

$v_i=1$ whenever $i \in \mathcal{X}_\text{K}$. If we are in Kyle’s winning state, then Kyle must necessarily win before Julien does.
$v_i=0$ whenever $i \in \mathcal{X}_\text{J}$. If we are in Julien’s winning state, then Kyle can’t possibly win before Julien does.
If $i \in \mathcal{X}_\text{T}$, the probability of Kyle winning depends recursively on which states we transition to next:
\[
v_i = \sum_{j = 1}^{64} \hat P_{ij} v_j \qquad\text{if }i \in \mathcal{X}_\text{T}
\]

The above gives us $64$ equations to solve in the $64$ unknowns $v_1,\dots,v_{64}$. Compactly, we can write the solution as:
\begin{align}
v_\text{K} &= 1 \\
v_\text{J} &= 0 \\
v_\text{T} &= \left(I-\hat P_{\text{T},\text{T}}\right)^{-1}
\hat P_{\text{T},\text{K}} \mathbf{1}
\end{align}Since each of the 64 starting states is equally likely, the probability that Kyle wins the game is simply the mean value of the vector $v$. Here is a short Matlab script that computes the solution:

P = 1/2*kron([1;1],kron(eye(4),[1 1]));
s = ["HHH","HHT","HTH","HTT","THH","THT","TTH","TTT"];
Kwin = "TTH";
Jwin = "TTT";
XK = find(kron(s == Kwin, s ~= Jwin));
XJ = find(kron(s ~= Kwin, s == Jwin));
XT = setdiff(1:64,union(XK,XJ));
Phat = kron(P,P);
v = zeros(64,1);
v(XK) = 1;
v(XJ) = 0;
v(XT) = (eye(length(XT)) - Phat(XT,XT)) \ Phat(XT,XK) * ones(length(XK),1);
mean(v)

The result of this computation is $\frac{45667}{71106} \approx 64.22\%$. So Kyle has a 64% chance of winning!

Extra credit

Here, we consider the same game as before, except we vary the winning state for Kyle and Julien. To do this, we loop through all possible combinations of winning states and run the same script as above. This yields the following matrix of probabilities. Rows and columns correspond to Kyle and Julien’s choice of winning sequence, respectively. Entries in the matrix are Kyle’s probability of winning.
\[
\begin{bmatrix}
\frac{1}{2} & \frac{25439}{71106} & \frac{1177}{2850} & \frac{25439}{71106} & \frac{25439}{71106} & \frac{1177}{2850} & \frac{25439}{71106} & \frac{1}{2} \\
\frac{45667}{71106} & \frac{1}{2} & \frac{28675}{52162} & \frac{1}{2} & \frac{1}{2} & \frac{28675}{52162} & \frac{1}{2} & \frac{45667}{71106} \\
\frac{1673}{2850} & \frac{23487}{52162} & \frac{1}{2} & \frac{23487}{52162} & \frac{23487}{52162} & \frac{1}{2} & \frac{23487}{52162} & \frac{1673}{2850} \\
\frac{45667}{71106} & \frac{1}{2} & \frac{28675}{52162} & \frac{1}{2} & \frac{1}{2} & \frac{28675}{52162} & \frac{1}{2} & \frac{45667}{71106} \\
\frac{45667}{71106} & \frac{1}{2} & \frac{28675}{52162} & \frac{1}{2} & \frac{1}{2} & \frac{28675}{52162} & \frac{1}{2} & \frac{45667}{71106} \\
\frac{1673}{2850} & \frac{23487}{52162} & \frac{1}{2} & \frac{23487}{52162} & \frac{23487}{52162} & \frac{1}{2} & \frac{23487}{52162} & \frac{1673}{2850} \\
\frac{45667}{71106} & \frac{1}{2} & \frac{28675}{52162} & \frac{1}{2} & \frac{1}{2} & \frac{28675}{52162} & \frac{1}{2} & \frac{45667}{71106} \\
\frac{1}{2} & \frac{25439}{71106} & \frac{1177}{2850} & \frac{25439}{71106} & \frac{25439}{71106} & \frac{1177}{2850} & \frac{25439}{71106} & \frac{1}{2}
\end{bmatrix}
\]Or, in numerical form,
\[
\small
\begin{bmatrix}
0.5 & 0.357762 & 0.412982 & 0.357762 & 0.357762 & 0.412982 & 0.357762 & 0.5 \\
0.642238 & 0.5 & 0.54973 & 0.5 & 0.5 & 0.54973 & 0.5 & 0.642238 \\
0.587018 & 0.45027 & 0.5 & 0.45027 & 0.45027 & 0.5 & 0.45027 & 0.587018 \\
0.642238 & 0.5 & 0.54973 & 0.5 & 0.5 & 0.54973 & 0.5 & 0.642238 \\
0.642238 & 0.5 & 0.54973 & 0.5 & 0.5 & 0.54973 & 0.5 & 0.642238 \\
0.587018 & 0.45027 & 0.5 & 0.45027 & 0.45027 & 0.5 & 0.45027 & 0.587018 \\
0.642238 & 0.5 & 0.54973 & 0.5 & 0.5 & 0.54973 & 0.5 & 0.642238 \\
0.5 & 0.357762 & 0.412982 & 0.357762 & 0.357762 & 0.412982 & 0.357762 & 0.5 \\
\end{bmatrix}
\]Note that the $(7,8)$ cell corresponds to the case we solved in the first part. Many of the cells contain $\frac{1}{2}$; these are cases where both Kyle and Julien have equally easy (or difficult) winning states that they are trying to reach. If we don’t know ahead of time which winning sequence Julien will pick, then we should average the rows of this matrix to see Kyle’s average chance of winning the game for each possible choice of winning sequence. This yields:
\[
\begin{array}{ccc}
\text{Kyle’s} & \text{Prob that} & \text{numerical} \\[-1mm]
\text{sequence} & \text{Kyle wins} & \text{approximation} \\ \hline
\textsf{HHH} & \frac{6875419}{16887675} & 0.407126\\
\textsf{HHT} & \frac{508129861}{927257793} & 0.547992\\
\textsf{HTH} & \frac{18467111}{37165425} & 0.496890\\
\textsf{HTT} & \frac{508129861}{927257793} & 0.547992\\
\textsf{THH} & \frac{508129861}{927257793} & 0.547992\\
\textsf{THT} & \frac{18467111}{37165425} & 0.496890\\
\textsf{TTH} & \frac{508129861}{927257793} & 0.547992\\
\textsf{TTT} & \frac{6875419}{16887675} & 0.407126\\
\end{array}
\]Note that several rows are repeated in the matrices above. This is because different choices of winning state can be considered “equivalent” due to the symmetry inherent in this game. Specifically, we can define the following three equivalence classes:
\begin{align}
A &: \{\textsf{HHH},\textsf{TTT}\} \\
B &: \{\textsf{HTH},\textsf{THT}\} \\
C &: \{\textsf{HHT},\textsf{HTT},\textsf{THH},\textsf{TTH}\}
\end{align}We can simplify the table above to:
\[
\begin{array}{ccc}
\text{Equivalence} & \text{Prob that} & \text{numerical} \\[-1mm]
\text{class} & \text{Kyle wins} & \text{approximation} \\ \hline
A & \frac{6875419}{16887675} & 0.407126\\
B & \frac{18467111}{37165425} & 0.496890\\
C & \frac{508129861}{927257793} & 0.547992\\
\end{array}
\]Based on this analysis, it appears that Kyle picking the $B$ equivalence class ($\small\textsf{HTH}$ or $\small\textsf{THT}$) leads to the game being close to 50% on average. This means that if Julien picks his winning sequence at random, the game will be, on average, close to perfectly fair.

Alternative interpretations

Different notions of “close to 1/2” lead to different choices for Kyle. For example, if our notion of fairness is “equal to 1/2 as often as possible”, then Kyle should pick option $C$. In this case, if Julien picks randomly, the game will be perfectly fair 50% of the time. But in option $B$, the game is only perfectly fair 25% of the time.

Alternatively, what if Julien does not pick randomly, but rather is greedy and picks his winning sequence in an effort to maximize his chance of winning? And what if Kyle knows that Julien is doing this, and tries to maximize his chances of winning as well? And what if Julien knows that Kyle knows? (and so on…) Then it becomes a game in the mathematical sense of the word, and the solution is the Nash equilibrium. To find the Nash equilibrium, rewrite the probability of Kyle winning as a reduced $3\times 3$ matrix based on equivalence class:
\[
\begin{bmatrix}
\frac{1}{2} & \frac{1177}{2850} & \frac{25439}{71106} \\
\frac{1673}{2850} & \frac{1}{2} & \frac{23487}{52162} \\
\frac{45667}{71106} & \frac{28675}{52162} & \frac{1}{2}
\end{bmatrix}
\approx
\begin{bmatrix}
0.5 & 0.412982 & 0.357762 \\
0.587018 & 0.5 & 0.45027 \\
0.642238 & 0.54973 & 0.5
\end{bmatrix}
\]
In this game, Kyle picks a row and Julien picks a column. Kyle wants to maximize the resulting number and Julien wants to minimize it. By inspection, we can see that the best choice for both Kyle and Julien is again to pick option $C$.

So perhaps $C$ is a more sensible choice for Kyle after all!

The Likeliest Monopoly Square

This week’s Fiddler is about rolling dice in the board game Monopoly.

We have a square board with 40 individual spaces around it, numbered from 0 to 39. All players begin on space 0 (akin to the “Go” square in Monopoly) and roll a pair of dice to determine how many spaces they advance each turn. However, unlike Monopoly, there is no way to otherwise advance around the board (i.e., there’s no “Chance,” “Community Chest,” going to jail, etc.). In their first pass around the board, which space from 1 to 39 are players most likely to land on at some point (i.e., not necessarily on their first or last roll, but after any number of rolls)?

Extra Credit
The square board has 10 spaces on each side. The first side has spaces 0 through 9, the second side has spaces 10 through 19, the third side has spaces 20 through 29, and the fourth side has spaces 30 through 39. Because you’re rolling two dice, it’s impossible to land on space 1 in your first pass around the board. Several other spaces on the first side of the board are similarly unlikely. Putting that first side of the board aside, which space from 10 to 39 are players least likely to land on at some point during their first pass around the board? (Another question: What if you rolled three dice at a time instead of two?)

My solution:
[Show Solution]

An efficient way to tackle such a problem is via the use of polynomials, or more generally, generating functions. Here is the basic idea. Suppose we roll a single die. Each number $1,\dots,6$ has equal probability $\frac{1}{6}$ of occurring. We will encode this using the polynomial:
\[
q(x) = \tfrac{1}{6}x + \tfrac{1}{6}x^2 + \tfrac{1}{6}x^3 + \tfrac{1}{6}x^4 + \tfrac{1}{6}x^5 + \tfrac{1}{6}x^6
\]We can interpret the coefficient of $x^k$ as the probability that we will land on the $k^\text{th}$ square. If we roll two dice, the corresponding probability distribution is obtained by taking the product of the polynomials for each die. That is, we have:
\begin{multline}
q(x)^2 =
\tfrac{1}{36}x^2+\tfrac{1}{18}x^3+\tfrac{1}{12}x^4+\tfrac{1}{9}x^5+\tfrac{5}{36}x^6+\tfrac{1}{6}x^7\\+\tfrac{5}{36} x^8+\tfrac{1}{9}x^9+\tfrac{1}{12}x^{10}+\tfrac{1}{18}x^{11}+\tfrac{1}{36}x^{12}
\end{multline}So the probability of rolling a 6 would be $\tfrac{5}{36}$. This works because of how polynomials multiply. We can roll a 6 in several ways:
\[
(1+5), (2+4), (3+3), (4+2), (5+1)
\]Likewise, when we multiply the two polynomials, the coefficient of $x^6$ can be obtained by selecting the same corresponding pairs of terms from the first and second polynomial:
\[
(x^1\cdot x^5),
(x^2\cdot x^4),
(x^3\cdot x^3),
(x^4\cdot x^2),
(x^5\cdot x^1)
\]The probability of each of these occurring is the product of the coefficients from each polynomial, so the net probability is the sum of each such product, i.e. the coefficient of the $x^6$ term in the resulting polynomial.

Since we can also obtain a particular number by rolling 4 dice, 6 dice, 8 dice, etc., what we’re after is the coefficient of $x^n$ in the infinite polynomial:
\[
p(x) = \sum_{k=1}^\infty q(x)^{2k}
= \frac{q(x)^2}{1-q(x)^2}.
\]The compact representation comes from summing the geometric series, but in this case it’s not that helpful because we still need to extract the polynomial coefficients when it is expanded as a series. We can evaluate coefficients numerically in Mathematica and produce a plot of the coefficients of the resulting polynomial. Here is the code:

q = 1/6 (x + x^2 + x^3 + x^4 + x^5 + x^6);
BarChart[
 CoefficientList[Normal[Series[q^2/(1 - q^2), {x, 0, 40}]], x], 
 ChartLabels -> Table[k, {k, 0, 40}], GridLines -> Automatic]

and the resulting plot:

We can now answer all the questions in the original problem.

The likeliest square to land on is 7 (lucky number 7!). The associated probability is $\frac{1417}{7776} \approx 0.182227$.
Excluding squares 0-9, the least likely square to land on is 13 (unlucky number 13!). The associated probability is $\frac{22621045}{181398528} \approx 0.124704$.

If we are rolling three dice instead, we simply change $q(x)^2$ to $q(x)^3$. The resulting plot is:

The reason for this curious oscillating shape is that it comes from summing the distributions for $q(x)^3$, $q(x)^6$, $q(x)^9$, etc. Each of these tends to a normal distribution with increasing variance. Here are the distributions plotted separately (3 dice at once):

Analytic solution

The solution above gives us a way to compute the solution for any $n$ by multiplying and adding polynomials together. But can we obtain an actual formula for the probability of landing on a particular square? Yes! But it will be complicated… It involves rewriting our infinite series as a sum of simpler series using partial fraction decomposition:
\begin{align}
p(x) &= \sum_{k=1}^\infty q(x)^{2k} = \frac{q(x)^2}{1-q(x)^2} \\
&= -1 + \frac{\tfrac{1}{2}}{1-q(x)} + \frac{\tfrac{1}{2}}{1+q(x)}\\
&= -1-\frac{3}{x^6+x^5+x^4+x^3+x^2+x-6}+\frac{3}{x^6+x^5+x^4+x^3+x^2+x+6}\\
&= -1-\frac{3}{(x-1)(x^5+2 x^4+3 x^3+4 x^2+5 x+6)}+\frac{3}{x^6+x^5+x^4+x^3+x^2+x+6}\\
&= -1-\frac{1}{7(x-1)} + \frac{x^4+3 x^3+6 x^2+10 x+15}{7 \left(x^5+2 x^4+3 x^3+4 x^2+5 x+6\right)}+\frac{3}{x^6+x^5+x^4+x^3+x^2+x+6}
\end{align}We continue by breaking up the polynomials even further.
\begin{align}
&\text{Let $\{\eta_k\}$ be the 5 roots of the polynomial $x^5+2x^4+3x^3+4x^2+5x+6$}\\
&\text{Let $\{\zeta_k\}$ be the 6 roots of the polynomial $x^6+x^5+x^4+x^3+x^2+x+6$}
\end{align}Applying yet another partial fraction decomposition, we obtain:
\[
p(x) = -1-\frac{1}{7(x-1)}
+\sum_{k=1}^5 \frac{a_k}{x-\eta_k}+\sum_{k=1}^6 \frac{b_k}{x-\zeta_k}
\]where the $a_k$ and $b_k$ coefficients are given by:
\begin{align}
a_k &= \frac{\eta_k^4+3\eta_k^3+6\eta_k^2+10\eta_k+15}{7\prod_{i\neq k}(\eta_k-\eta_i)}, &
b_k &= \frac{3}{\prod_{i\neq k}(\zeta_k-\zeta_i)}
\end{align}Rearranging the terms of $p(x)$ and expanding each term using the familiar geometric series $\frac{1}{1-rx}=\sum_{n=0}^\infty r^kx^k$, we obtain
\begin{align}
p(x) &= -1 + \frac{1}{7}\left( \frac{1}{1-x} \right)
-\sum_{k=1}^5 \frac{a_k}{\eta_k} \left(\frac{1}{1-\eta_k^{-1}x}\right)
-\sum_{k=1}^6 \frac{b_k}{\zeta_k} \left( \frac{1}{1-\zeta_k^{-1}x} \right) \\
&= -1+\sum_{n=0}^\infty\underbrace{\left(
\frac{1}{7}-\sum_{k=1}^5 \frac{a_k}{\eta_k^{n+1}}-\sum_{k=1}^6\frac{b_k}{\zeta_k^{n+1}} \right)}_{q_n} x^n
\end{align}And that’s it! Our formula for the probability that a player lands on the $n^\text{th}$ square during their first pass around the board (rolling 2 dice at a time) is given by the $q_n$ coefficient in the infinite sum above.

It’s also possible to derive a formula for 3 dice (or any number of dice) at once using the same approach, but the formulas get progressively more complicated as we add more dice.

What about the limit as $n\to\infty$? In other words, what is the probability of landing on any square in the limit as we get very far from the starting point? This can be calculated easily from the formula for $q_n$. As $n\to\infty$, the complicated terms drop out, and we obtain $q_n \to \frac{1}{7}\approx 0.142857$. In the case where we roll 3 dice at once, the probability turns out to be $\frac{2}{21}\approx 0.095238$.

And what about for more dice? The probability $\lim_{n\to\infty}q_n$ if $k$ dice are rolled at once is simply the residue of $p(x)$ for the pole at $x=1$. This follows because all other poles will have roots with magnitude less than 1, so they decay to zero as $n\to\infty$. To calculate the residue, we use the fact that $q(1)=1$ and apply L’Hôpital’s rule:
\begin{align}
&\lim_{x\to 1} (1-x) p(x) \\
&= \lim_{x\to 1} \frac{(1-x)q(x)^k}{1-q(x)^k} \\
&= \lim_{x\to 1} \frac{1-x}{1-\left(\frac{x+x^2+x^3+x^4+x^5+x^6}{6}\right)^k} \\
&= \lim_{x\to 1} \frac{-1}{-k \left(\frac{x+x^2+x^3+x^4+x^5+x^6}{6}\right)^{k-1}\left( \frac{1+2x+3x^2+4x^3+5x^4+6x^5}{6} \right)} \\
&= \frac{2}{7k}
\end{align}And we can see that this formula agrees with what we found for $k=1$ and $k=2$. Now I’m wondering if there isn’t a way to arrive at this formula purely from first principles!

Picking a speaker at random

This week’s Fiddler is about selecting a speaker of the house at random. How long will it take?

There are three candidates who want the job of Speaker. All 221 members of the party vote by picking randomly from among the candidates. If one candidate earns the majority of the votes, they become the next Speaker. Otherwise, the candidate with the fewest votes is eliminated and the process is repeated with one less candidate. If two or more candidates receive the same smallest number of votes, then exactly one of them is eliminated at random. What is the average number of rounds needed to select a new Speaker?

Extra Credit
What if there were 10 candidates running for Speaker?

My solution:
[Show Solution]

We will consider a more general case; suppose there are $n$ candidates and $2m+1$ voting members in the party. We write $2m+1$ since there should be an odd number of party members to break ties. Define the following probability:
\[
p(n,m) = \left\{
\begin{array}{l}
\text{Probability that one of the $n$ candidates} \\
\text{wins a majority of the $2m+1$ votes.}
\end{array}
\right\}
\]Let’s find a formula for $p(n,m)$. For a particular candidate to win $k$ votes, $k$ members must vote for the candidate and the remaining members must vote for one of the other $n-1$ candidates. This occurs with probability $\binom{2m+1}{k}\left(\frac{1}{n}\right)^k \left(1-\frac{1}{n}\right)^{2m+1-k}$. The candidate wins a majority if they garner at least $m+1$ of the $2m+1$ votes. We multiply the final result by $n$ since we are counting the probability that any of the $n$ candidates wins a majority. Therefore,
\[
p(n,m) = n\sum_{k=m+1}^{2m+1} \binom{2m+1}{k}\left(\frac{1}{n}\right)^k \left(1-\frac{1}{n}\right)^{2m+1-k}
\]Note: We can add the probabilities of the different candidates winning as we did because the events are disjoint; it is impossible for two candidates to simultaneously win a majority of the votes. If the criteria allowed for overlap, e.g., one candidate winning at least 40% of the vote, then two candidates could win simultaneously, and we would have to account for this case to avoid double-counting.

Making the change of variables $k \mapsto 2m+1-k$, we obtain

$\displaystyle
p(n,m) = n \sum_{k=0}^m \binom{2m+1}{k} \left(1-\frac{1}{n}\right)^{k}\left(\frac{1}{n}\right)^{2m+1-k}
$

We recognize this as $n$ times the CDF of a binomial distribution with $2m+1$ trials and probability $1-\frac{1}{n}$.

What we’re really after is the expected number of rounds required for a candidate to secure a majority of the vote.
\[
f(n,m) = \left\{
\begin{array}{l}
\text{Expected number of rounds until one of the } \\
\text{that one of the $n$ candidates is elected.}
\end{array}
\right\}
\]This will happen in one round with probability $p(n,m)$. It will happen in $k$ rounds if it does not happen in the first $k-1$ rounds, and happens in the $k^\text{th}$ round. Putting this together, we get:

$\displaystyle
f(n,m) = \sum_{k=1}^{n-1} k\cdot p(n-k+1,m) \prod_{i=0}^{k-1} \bigl( 1-p(n-k+1+i,m) \bigr)
$

Unfortunately, these formulas do not simplify in any meaningful way, so we must evaluate them numerically. The first question asks us to compute $f(3,110)$ and the extra credit asks for $f(10,110)$. Calculating these, we obtain:
\begin{align}
f(3,110) &\approx 1.9999995110943353682 \\
f(10,110) &\approx 8.9999995110943296464
\end{align}These are rational numbers, of course, but $f(3,110)$ and $f(10,110)$ have 105-digit and 1357-digit denominators, respectively, when expressed as irreducible fractions, so I will not write them out here! One thing for sure is that it is very unlikely that the voting will end early. It is a virtual certainty that it will go the full $n-1$ rounds until there are only two candidates remaining.

Bonus: Changing the number of voters

The reason it is so difficult to elect a Speaker early is that there are so many voting members, and therefore it is unlikely that one candidate will obtain a majority of the votes by random chance. We can see how this probability would scale if we changed the number of voting members. Here is a plot of $f(n,m)$ for different values of $n$ (number of candidates), as a function of $2m+1$ (number of voting members). In all cases, we see that the expected number of rounds rapidly approaches $n-1$.

Pinball probability

This week’s Fiddler is a challenging probability question.

You’re playing a game of pinball that includes four lanes, each of which is initially unlit. Every time you flip the pinball, it passes through exactly one of the four lanes (chosen at random) and toggles that lane’s state. So if that lane is unlit, it becomes lit after the ball passes through. But if the lane is lit, it becomes unlit after the ball passes through. On average, how many times will you have to flip the pinball until all four lanes are lit?

Extra credit: Instead of four lanes, now suppose your pinball game has $n$ lanes. And let’s say that $E(n)$ represents the average number of pinball flips it takes until all $n$ lanes are lit up. Now, each time you increase the number of lanes by one, you find that it takes you approximately twice as long to light up all the lanes. In other words, $E(n+1)$ seems to be about double $E(n)$. But upon closer examination, you find that it’s not quite double. Moreover, there’s a particular value of $n$ where the ratio $E(n+1)/E(n)$ is at a minimum. What is this value of $n$?

My solution:
[Show Solution]

Qualitative behavior and derivation of $E(n)$

Let’s fix $n$ for now, and define $T_k$ to be the number of expected additional flips required to light up all lanes given that there are currently $k$ lanes lit up. Our goal is to find $E(n) = T_0$. Two things can happen every time we flip the pinball:

We pass through a lit lane, which causes it to become unlit. This happens with probability $\tfrac{k}{n}$, leaving $k-1$ lit lanes.
We pass through an unlit lane, which causes it to light up. This happens with probability $\tfrac{n-k}{n}$, leaving $k+1$ lit lanes.

The boundary condition is that $T_n=0$, because if all lanes are lit, we have accomplished the goal and there are no more flips required. Putting everything together, we obtain the following recursion for $T_k$:
\begin{align}
T_0 &= T_1 + 1 \\
T_k &= \tfrac{k}{n}T_{k-1} + \tfrac{n-k}{n}T_{k+1} + 1 && \text{for }k=1,\dots,n-1.\\
T_n &= 0
\end{align}This is a set of $n+1$ linear equations in the $n+1$ unknowns $\{T_0,\dots,T_n\}$, so we have enough information to solve the problem and find what we’re after, which is $T_0 = E(n)$.

If you’re thinking “hey wait a minute… isn’t this just MFPT for an Ehrenfest Markov chain?”, I got you! I wanted to solve the problem without bringing in any heavy machinery first. I added a section at the end of this post that explains a bunch of connections to Markov chains.

Ok, back to our system of equations. Let’s start by making the substitution
\[
\Delta_k := T_{k-1}-T_k\qquad\text{for }k=1,\dots,n.
\]Since $T_n=0$, we can write $T_0$ in terms of the $\Delta_k$’s as:
\[
E(n) = T_0 = \Delta_1 + \Delta_2 + \cdots + \Delta_n
\]Meanwhile, we can rewrite the recursion in terms of the $\Delta_k$’s as:
\begin{align}
&& T_k &= \tfrac{k}{n}T_{k-1} + \tfrac{n-k}{n}T_{k+1} + 1\\
\iff && T_k-T_{k+1} &= \tfrac{k}{n} \left( T_{k-1}-T_{k+1} \right) + 1 \\
\iff && \Delta_{k+1} &= \tfrac{k}{n}\left( \Delta_{k}+\Delta_{k+1} \right) + 1 \\
\iff && \left(1-\tfrac{k}{n}\right)\Delta_{k+1} &= \tfrac{k}{n}\Delta_k + 1 \\
\iff && \Delta_{k+1} &= \tfrac{k}{n-k}\Delta_k + \tfrac{n}{n-k}
\end{align}So our simplified recursion looks like:
\begin{align}
\Delta_1 &= 1 \\
\Delta_{k+1} &= \tfrac{k}{n-k}\Delta_k + \tfrac{n}{n-k} && \text{for }k=1,\dots,n-1
\end{align}We can now solve for each $\Delta_k$ by forward substitution. The first few terms look like (rearranged to highlight the pattern):
\begin{align}
\Delta_1 &= \tfrac{n}{n} \\
\Delta_2 &= \tfrac{n}{(n-1)n} + \tfrac{n}{n-1} \\
\Delta_3 &= \tfrac{2\cdot 1 n}{(n-2)(n-1)n} + \tfrac{2n}{(n-2)(n-1)} + \tfrac{n}{n-2} \\
\Delta_4 &= \tfrac{3\cdot 2 \cdot 1n}{(n-3)(n-2)(n-1)n} + \tfrac{3\cdot 2 n}{(n-3)(n-2)(n-1)} + \tfrac{3n}{(n-3)(n-2)} + \tfrac{n}{n-3}
\end{align}The pattern is clear, and we have the following formula in terms of factorials, which we can rewrite in terms of binomial coefficients:
\begin{align}
T_k &= \sum_{i=0}^{k-1} \frac{n(k-1)!/i!}{(n-i)!/(n-k)!} = \sum_{i=0}^{k-1} \frac{\binom{n}{i}}{\binom{n-1}{k-1}}\\
%&= \sum_{i=0}^{k-1} \frac{n(k-1)! (n-k)!}{(n-i)!i!} \\
%&= \sum_{i=0}^{k-1} \frac{n! (k-1)! (n-k)!}{(n-1)!(n-i)!i!} \\
\end{align}Since $E(n)$ is the sum of the $\Delta_k$’s, we can now obtain a compact formula for the average number of pinball flips it takes until all $n$ lanes are lit up.

$\displaystyle
E(n) = \sum_{k=1}^n \sum_{i=0}^{k-1} \frac{\binom{n}{i}}{\binom{n-1}{k-1}}
$

Let’s try to make some sense of this complicated formula. We can manipulate it a bit to first obtain:
\[
E(n) = 2^n \sum_{k=1}^n \underbrace{\frac{1}{\binom{n-1}{k-1}}}_{\alpha_k} \,\underbrace{\sum_{i=0}^{k-1} \frac{1}{2^n}\binom{n}{i} }_{\Phi_k}
\]So we can think of $E(n)$ as a dot product between two sequences $\{\alpha_k\}$ and $\{\Phi_k\}$, scaled by $2^n$. But what do these sequences represent?

Recall that the binomial distribution with $p=1/2$ tends to a normal distribution with mean $n/2$ and variance $n/4$ as $n$ gets large, so we have:
\[
\frac{1}{2^n} \binom{n}{i} \approx \frac{1}{\sqrt{\pi n/2}}\exp\left( -\frac{(i-n/2)^2}{n/2} \right)
\]Therefore, the term $\Phi_k$ is approximately the cumulative distribution function for a normal distribution. Meanwhile, the weights $\alpha_k$ are weights that start at $1$ when $k=1$, reach a minimum when $k=\tfrac{n+1}{2}$, and then increase back to $1$ as $k$ approaches $n-1$. Here is what these sequences look like when $n=\{8,20,50\}$:

As we can see, the weights get steeper and steeper. The last two weights are
\[
\alpha_{n-2} = \binom{n-1}{n-2}^{-1} = \frac{1}{n-1}\qquad\text{and}\qquad
\alpha_{n-1} = \binom{n-1}{n-1}^{-1} = 1
\]Therefore, as $n$ gets large, we have
\[
E(n) \approx 2^n\left(1+\frac{1}{n-1}\right) = 2^n \frac{n}{n-1}
\]
This explains why we can expect that
\[
\lim_{n\to\infty} \frac{E(n+1)}{E(n)} = 2
\]
Unfortunately, I don’t see any way to make meaningful simplifications to the expression for $E(n)$, so we’ll have to proceed numerically to answer the questions asked.

The answer

The first problem asks us to calculate $E(4)$, which is easy enough:

$\displaystyle
E(4) = \frac{64}{3} \approx 21.3333
$

We can confirm the answer is correct via Monte Carlo simulation, and here is the histogram we obtain after $10^7$ trials.

As we can see, the empirical average of $21.3291$ is in good agreement with the true value of $21.3333$ we computed. The histogram only has bars at even values $n$ since it is impossible to go from $0$ to $4$ in an odd number of steps (at every step, we toggle the parity, so since we want to go from an even number to an even number, we must move an even number of times).

Extra credit

We can plot the function $\frac{E(n+1)}{E(n)}$ as a function of $n$, and we obtain:

Sure enough, the values are close to $2$ (I started the plot at $n=4$ to see more detail), and the minimum value occurs at $n=6$. This minimum value is:

$\displaystyle
\frac{E(7)}{E(6)} = \frac{151}{78} \approx 1.9359
$

The solution method above was relatively straightforward, but it turns out this problem has a Markov chain interpretation, and it’s actually a well-studied type of model (enough that it has a name!). If you’re interested, read on…

Markov chain interpretation

We can also interpret this problem as a Markov chain. Here, the state $k\in\{0,1,\dots,n\}$ represents having $k$ lanes lit up. The transition probabilities are represented below:

and the associated transition matrix is:
\[
P = \begin{bmatrix}
0 & 1 & 0 & 0 & \cdots & 0 & 0 \\
\tfrac{1}{n} & 0 & \tfrac{n-1}{n} & 0 & \cdots & 0 & 0 \\
0 & \tfrac{2}{n} & 0 & \tfrac{n-2}{n} & \ddots & \vdots & \vdots \\
0 & 0 & \tfrac{3}{n} & 0 & \ddots & 0 & 0 \\
\vdots & \vdots & \ddots & \ddots & \ddots &\tfrac{2}{n} & 0 \\
0 & 0 & \cdots & 0 & \tfrac{n-1}{n} & 0 & \tfrac{1}{n} \\
0 & 0 & \cdots & 0 & 0 & 1 & 0
\end{bmatrix}
\]This particular transition matrix has a name! It’s called the Ehrenfest Model, and it is used to model diffusion, among other things. The standard setup is that you have two urns, and a total of $n$ marbles. Each marble has a probability of $1/2$ of spontaneously jumping from one urn to the other. In our problem, all $n$ marbles start in the first urn, and we are asked to find the expected number of jumps required before all marbles first simultaneously end up in the other urn.

The expected number of steps required to first reach a given state when starting at another given state in a Markov chain is called the Mean First Passage Time (MFPT). It is also sometimes called the First Hitting Time. The way to calculate the MFPT is to solve the system of equations
\[
\hat{P}x + \mathbf{1} = x,
\]where $\mathbf{1}$ is the vector of 1’s, and $\hat P$ is the transition matrix with rows and columns corresponding to the target set removed. In our case, the target set is the last state, so we should remove the last row and column. Then, $x$ will be the vector of MFPT’s from each of the other states. In our case, we are interested in the first component of $x$. The system of equations above is precisely the one we derived earlier, so nothing new here; just another way to get the same answer!

Optimal baseball lineup

This week’s Fiddler is a problem about how to set the optimal baseball lineup.

Eight of your nine batters are “pure contact” hitters. One-third of the time, each of them gets a single, advancing any runners already on base by exactly one base. (The only way to score is with a single with a runner on 3rd). The other two-thirds of the time, they record an out, and no runners advance to the next base. Your ninth batter is the slugger. One-tenth of the time, he hits a home run. But the remaining nine-tenths of the time, he strikes out. Your goal is to score as many runs as possible, on average, in the first inning. Where in your lineup (first, second, third, etc.) should you place your home run slugger?

Extra Credit: Instead of scoring as many runs as possible in the first inning, you now want to score as many runs as possible over the course of nine innings. What’s more, instead of just having one home run slugger, you now have two sluggers in your lineup. The other seven batters remain pure contact hitters. Where in the lineup should you place your two sluggers to maximize the average number of runs scored over nine innings?

My solution:
[Show Solution]

We will attack this problem using dynamic programming. Rather than solving for the average number of runs for a particular configuration (of outs, runners on base, current batter, etc.), we will solve for all possible configurations! To this effect define the following:
\[
V(o,r,b,i) = \left\{
\begin{array}{l}
\text{Expected number of runs at end of}\\
\text{game given that we are currently in}\\
\text{the $i^\text{th}$ inning with $o$ outs, $r$ runners}\\
\text{on base, and the $b^\text{th}$ batter is up.}
\end{array}\right\}
\]Each index can only take on certain values. Define the sets:
\begin{align}
O &= \{0,1,2,3\} && \text{(number of outs)} \\
R &= \{0,1,2,3\} && \text{(runners on base)} \\
B &= \{1,2,\dots,9\} && \text{(batter currently up)} \\
I &= \{1,2,\dots,9\} && \text{(current inning)}
\end{align}Therefore there are $4\times 4\times 9\times 9 = 1296$ possible configurations we must solve for. Finally, some subset of the batters are the sluggers. Let this subset be $S \subseteq B$. The remaining batters are the pure contact hitters, which we denote $\bar S$. Therefore $S\cap \bar S = \emptyset$ and $S \cup \bar S = B$.

Finally, define $p$ to be the on-base average for pure contact hitters, and $q$ to be the slugging percentage for the sluggers. In the problem statement, we have $p=\frac{1}{3}$ and $q=\frac{1}{10}$.

We can now write equations that determine how the various configuration are related to each other:

When we have 3 outs, clear bases and advance to next inning, or end the game if we are in the 9th inning.\begin{align}
V(3,r,b,i) &= V(0,0,b,i+1) && \forall r\in R, b\in B, i\in\{1,\dots,8\}\\
V(3,r,b,9) &= 0 && \forall r\in R, b\in B
\end{align}This is a total of $4\times 9\times 9 = 324$ equations.
When a pure contact hitter is up to bat, they either advance runners (and possibly score if bases are loaded), or they strike out.
\begin{align}
V(o,r,b,i) &= p V(o,r+1,b+1,i) + (1-p) V(o+1,r,b+1,i) && \forall o\in \{0,1,2\}, r\in \{0,1,2\}, b\in \bar S, i\in I \\
V(o,3,b,i) &= p (1+V(o,3,b+1,i)) + (1-p) V(o+1,3,b+1,i) && \forall o\in \{0,1,2\}, b\in \bar S, i \in I
\end{align}This is a total of $3\times 4\times 9\times |\bar S| = 108|\bar S|$ equations.
When a slugger is up to bat, they either hit a home run (and clear the bases), or they strike out.
\begin{align}
V(o,r,b,i) &= q (r+1+V(o,0,b+1,i)) + (1-q) V(o+1,r,b+1,i) && \forall o\in\{0,1,2\}, r\in R, b\in S, i\in I
\end{align}This is a total of $3\times 4\times 9\times |S| = 108|S|$ equations.

In all the equations above, when we write “$b+1$”, we mean that if we get to the end of the order, we should cycle back to the top of the order.
This adds up to a total of $324+108|\bar S|+108|S| = 1296$ equations. We have the same number of equations as variables, which tells us that we accounted for everything!

I used Mathematica to input all the equations programmatically and solve them, and here is what I found.

First part: one inning, one slugger

When looking at a single inning with a single slugger, we pick $I = \{1\}$, and we can choose $S = \{9\}$ for example, and then ask, which batter should bat first to maximize score? this amounts to see which of
\[
V(0,0,b,1),\qquad b=1,2,\dots,9
\]is largest. Using the hitting percentages from the problem ($p=\frac{1}{3}$ and $q=\frac{1}{10}$), we obtain the following bar chart:

So it is best to put the slugger fourth in the batting order. This leads to an expected score of $\frac{3051922249868633}{12708748019768805} \approx 0.240143$ runs scored in the first inning. This is only narrowly better than putting the slugger third, which leads to an average of $0.23937$ runs scored. Overall, here is what we find:
\[
\begin{array}{r|ccccccccc}
\text{slugger position} & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 \\ \hline
\text{expected runs} & 0.184128 & 0.210732 & 0.23937 & 0.240143 & 0.202726 & 0.156235 & 0.168336 & 0.174246 & 0.176961
\end{array}
\]

By changing $p$ and $q$, we get different solutions. For example, if we make the pure contact hitters better, say $p=\frac{2}{3}$, but we keep the slugging at $q=\frac{1}{10}$, we obtain:

So in this case, the slugger is best placed last, presumably because they get out so much more frequently than the pure contact hitters.

Second part: nine innings, two sluggers

With nine innings and two sluggers, I iterated through all possible two-element subsets $S \subseteq B$, and then computed $V(0,0,1,1)$ for each case. This led to the following picture:

This plot shows the expected number of runs after 9 innings as a function of the two sluggers’ positions in the batting order. The best outcome occurs when the sluggers are 3rd and 4th in the batting order. This is the yellow bar in the (3,4) coordinate, which corresponds to an average of 1.99771 runs per game. Note that the y-axis has a relatively small range; it only ranges from 1.85 to 2.0. This goes to show that over the course of 9 innings, the batting order does not make that much difference. Here is a table showing all the results:

Code

To solve this problem analytically, I used Mathematica to specify each of the equations above, and then called the “Solve” function on that system of equations. Each of the equations can be collected into a list using the “Table” function, which makes it easy to enumerate the equations without having to write each one down manually. You can view my Mathematica code here, or in the window below.