Markov chains – Book Proofs

Catch the grasshopper

This week’s Riddler classic is a probability problem about a grasshopper!

You are trying to catch a grasshopper on a balance beam that is 1 meter long. Every time you try to catch it, it jumps to a random point along the interval between 20 centimeters left of its current position and 20 centimeters right of its current position. If the grasshopper is within 20 centimeters of one of the edges, it will not jump off the edge. For example, if it is 10 centimeters from the left edge of the beam, then it will randomly jump to anywhere within 30 centimeters of that edge with equal probability (meaning it will be twice as likely to jump right as it is to jump left). After many, many failed attempts to catch the grasshopper, where is it most likely to be on the beam? Where is it least likely? And what is the ratio between these respective probabilities?

My solution:
[Show Solution]

Let’s suppose the beam has length $L$ and the grasshopper can jump a maximum distance $\frac{M}{2}$. Define $p(x)$ to be the probability density function describing how likely the grasshopper is to be at any given location $x\in\left[-\frac{L}{2},\frac{L}{2}\right]$.

We can arrive at $p(x)$ via a sequential process. Suppose that at time $k$, the grasshopper has distribution $p_k(x)$. When the grasshopper jumps, its probability distribution becomes $p_{k+1}(x)$. Our goal is to find the limiting distribution $p = \lim_{k\to\infty} p_k$.

At each step, the grasshopper’s density function gets convolved with an appropriate uniform distribution corresponding with where it might land. We can write this as:
\[
p_{k+1}(x) = \int_{-L/2}^{L/2} \Psi(x,y) p_k(y) \,\mathrm{d}y.
\]Here, for each $y$, the kernel $\Psi(x,y)$ is the density function (of $x$) corresponding to a uniform distribution between $\max\bigl(-\frac{L}{2},y-\frac{M}{2}\bigr)$ and $\min\bigl(\frac{L}{2},y+\frac{M}{2}\bigr)$, since this corresponds to the range of locations a grasshopper at location $y$ can jump to. Putting all this together, we get the formula:
\[
p_{k+1}(x) = \int_{\max\bigl(-\frac{L}{2},x-\frac{M}{2}\bigr)}^{\min\bigl(\frac{L}{2},x+\frac{M}{2}\bigr)} \frac{p_k(y)\,\mathrm{d}y}{\min\bigl(\frac{L}{2},y+\frac{M}{2}\bigr)-\max\bigl(-\frac{L}{2},y-\frac{M}{2}\bigr)}.
\]Yuck! To get a better handle on what this might look like, we can simulate it using the given problem parameters ($L=1$ and $M=0.4$) using units of meters. Here is what we obtain when recursively apply the above formula, assuming $p_0$ is a uniform distribution:

The steady-state distribution is the same no matter what we use as the starting distribution. Here is what it looks like when we start the grasshopper on the left edge of the beam. I had to use different axis scaling to make the transient behavior fit in the plot, but the limiting distribution in the same:

These simulations give us a big hint: The limiting distribution appears to be piecewise-linear! We already expected a symmetric distribution due to the symmetry in the problem, but now we can posit a density function of the form:
\[
p(x) = \begin{cases}
ax+b & -\frac{L}{2}\leq x \leq -\left|\frac{L-M}{2}\right| \\
c & -\left|\frac{L-M}{2}\right| < x < \left|\frac{L-M}{2}\right| \\ -ax+b & \left|\frac{L-M}{2}\right| \leq x \leq \frac{L}{2} \end{cases} \]We have three unknown parameters $(a,b,c)$, so we need three equations. First, $p$ is a probability density, so $\int_{-L/2}^{L/2}p(x)\,\mathrm{d}x=1$. Second, $p$ is continuous at the boundaries $x=\pm \left|\frac{L-M}{2}\right|$. Finally, we must satisfy $p(x) = p_{k+1}(x) = p_k(x)$ in the original recursion in our original recursion. Solving these equations, we find three cases: \[ p(x) = \begin{cases} \begin{cases} \frac{2 (L+M+2 x)}{M (4 L-M)} & -\frac{L}{2}\leq x \leq \frac{-L+M}{2} \\ \frac{4}{4 L-M} & \frac{-L+M}{2} < x < \frac{L-M}{2} \\ \frac{2 (L+M-2 x)}{M (4 L-M)} & \frac{L-M}{2} \leq x \leq \frac{L}{2} \end{cases} & 0 < M < L \\[3mm] \begin{cases} \frac{2 (L+M+2 x)}{M (4 L-M)} & -\frac{L}{2}\leq x \leq \frac{-M+L}{2} \\ \frac{4 L}{M (4 L-M)} & \frac{-M+L}{2} < x < \frac{M-L}{2} \\ \frac{2 (L+M-2 x)}{M (4 L-M)} & \frac{M-L}{2} \leq x \leq \frac{L}{2} \end{cases} & L < M < 2L \\[3mm] \quad\frac{1}{L} & 2L < M \end{cases} \] We can also calculate the ratio of largest to smallest likelihood in all three cases. The largest likelihood occurs in the center ($x=0$) while the smallest occurs at an edge ($x=\pm\frac{L}{2}$). We then obtain: \[ \frac{p_\text{max}}{p_\text{min}} = \begin{cases} 2 & 0 < M < L \\[1mm] \frac{2L}{M} & L < M < 2L \\[1mm] 1 & 2L < M \end{cases} \] So, amazingly, the ratio of largest likelihood to the smallest likelihood is 2, a constant indepedent of the length $L$ or the grasshopper range $M$, so long as $M < L$. In particular, the likelihood ratio is $2$ for the setting given in the problem statement ($L=1$ and $M=0.4$).

Tetrahedral dice game

This week’s Riddler Classic is a game of four-sided dice:

You have four fair tetrahedral dice whose four sides are numbered 1 through 4.

You play a game in which you roll them all and divide them into two groups: those whose values are unique, and those which are duplicates. For example, if you roll a 1, 2, 2 and 4, then the 1 and 4 will go into the “unique” group, while the 2s will go into the “duplicate” group.

Next, you reroll all the dice in the duplicate pool and sort all the dice again. Continuing the previous example, that would mean you reroll the 2s. If the result happens to be 1 and 3, then the “unique” group will now consist of 3 and 4, while the “duplicate” group will have two 1s.

You continue rerolling the duplicate pool and sorting all the dice until all the dice are members of the same group. If all four dice are in the “unique” group, you win. If all four are in the “duplicate” group, you lose.

What is your probability of winning the game?

My solution:
[Show Solution]

We can view this game as a Markov chain. At any point in time, we are in a particular state of the game, and when we re-roll the dice, we transition to a different state. Although there are many possible states (one for each possible way of rolling the four dice), we can use symmetry arguments to reduce the total states to just four. Here they are:

All dice are duplicates. This includes cases like (1,1,1,1) but also (2,2,3,3). If we ever arrive at this position, the game ends and we lose.
Three dice are duplicates. For example: (1,1,1,2)
Two dice are duplicates. For example: (1,1,2,3)
No dice are duplicates. For example: (1,2,3,4). If we ever arrive at this position, the game ends and we win.

We can associate each state transition with a probability. For example, suppose we roll (1,1,2,3). We are in the “two-duplicate” state. We must re-roll the two 1’s, and four things could happen:

We roll (2,2) or (3,3) (probability 1/8). We now have 3 duplicates and the game continues.
We roll (2,3) or (3,2) (probability 1/8). We now have 4 duplicates and the game ends (we lose).
We roll (1,4) or (4,1) (probability 1/8). We now have no duplicates and the game ends (we win).
We roll anything else (probability 5/8). We have 2 duplicates again and the game continues.

If we continue in this manner and find all possible transition probabilities between all four states, we obtain the following diagram:

We can think of the “three duplicates” state as our start state. Why? Because we start the game by rolling all four dice. This is mathematically equivalent to rolling three dice, since the value of the die we don’t roll might as well have been our first roll (all numbers are equivalent by symmetry). Therefore, we can represent our transitions as follows:
\[
A = \begin{bmatrix}
1 & \frac{5}{32} & \frac{1}{8} & 0 \\
0 & \frac{3}{16} & \frac{1}{8} & 0 \\
0 & \frac{9}{16} & \frac{5}{8} & 0 \\
0 & \frac{3}{32} & \frac{1}{8} & 1
\end{bmatrix},\qquad
x_0 = \begin{bmatrix}
0 \\ 1 \\ 0 \\ 0\end{bmatrix}.
\]Notice that all columns sum to 1 since the sum of probabilities from every node along outgoing edges must sum to 1. We can view our current state distribution as a column vector that sums to 1. For example, our initial state is given by $x_0$ above, since we start in the “three-duplicate” state with probability 1. Every time we transition to a new state, our new distribution can be found by multiplying our current state by $A$. In other words:
\[
x_{k+1} = A x_k,\quad\text{for }k=0,1,\dots
\]The question is: where will we end up on average? This is equivalent to asking for the limit
\[
x_\infty = \lim_{k\to\infty} A^k x_0
\]We can find this limit by performing an eigenvalue decomposition:
\begin{align}
A &= V\Lambda V^{-1} \\
&= \begin{bmatrix}
0 & 1 & \frac{23}{21} & -1 \\
0 & 0 & -\frac{8}{21} & 30 \\
0 & 0 & -\frac{12}{7} & -30 \\
1 & 0 & 1 & 1
\end{bmatrix}
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & \frac{3}{4} & 0 \\
0 & 0 & 0 & \frac{1}{16}
\end{bmatrix}
\begin{bmatrix}
0 & \frac{9}{20} & \frac{29}{60} & 1 \\
1 & \frac{11}{20} & \frac{31}{60} & 0 \\
0 & -\frac{21}{44} & -\frac{21}{44} & 0 \\
0 & \frac{3}{110} & -\frac{1}{165} & 0
\end{bmatrix}
\end{align}The limit is therefore:
\begin{align}
A^\infty x_0 &= V \Lambda^{\infty} V^{-1}x_0 \\
&= \begin{bmatrix}
0 & 1 & \frac{23}{21} & -1 \\
0 & 0 & -\frac{8}{21} & 30 \\
0 & 0 & -\frac{12}{7} & -30 \\
1 & 0 & 1 & 1
\end{bmatrix}
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix}
0 & \frac{9}{20} & \frac{29}{60} & 1 \\
1 & \frac{11}{20} & \frac{31}{60} & 0 \\
0 & -\frac{21}{44} & -\frac{21}{44} & 0 \\
0 & \frac{3}{110} & -\frac{1}{165} & 0
\end{bmatrix} \begin{bmatrix}
0 \\ 1 \\ 0 \\ 0\end{bmatrix}\\
&= \begin{bmatrix}
0 & 1 \\
0 & 0 \\
0 & 0 \\
1 & 0
\end{bmatrix}
\begin{bmatrix}
0 & \frac{9}{20} & \frac{29}{60} & 1 \\
1 & \frac{11}{20} & \frac{31}{60} & 0
\end{bmatrix}
\begin{bmatrix}
0 \\ 1 \\ 0 \\ 0\end{bmatrix} \\
&= \begin{bmatrix}
\frac{11}{20} \\ 0 \\ 0 \\ \frac{9}{20}\end{bmatrix}
\end{align}In other words, there is a 11/20, or 55% chance that we will lose and end up in the “all duplicate” case, and a 9/20, or 45% chance that we will win and end up in the “all unique” case.

The luckiest coin

This week’s Riddler Classic is about finding the “luckiest” coin!

I have in my possession 1 million fair coins. I first flip all 1 million coins simultaneously, discarding any coins that come up tails. I flip all the coins that come up heads a second time, and I again discard any of these coins that come up tails. I repeat this process, over and over again. If at any point I am left with one coin, I declare that to be the “luckiest” coin.

But getting to one coin is no sure thing. For example, I might find myself with two coins, flip both of them and have both come up tails. Then I would have zero coins, never having had exactly one coin.

What is the probability that I will at some point have exactly one “luckiest” coin?

Here is my solution:
[Show Solution]

Let’s define $p_n$ to be the probability that we will end up with exactly one “luckiest” coin, given that we start with $n$ coins. Our terminal conditions are that $p_0=0$ and $p_1=1$, since we lose if we end up with zero and we win if we end up with one. We can define the rest of the $p_n$ recursively as follows. If we flip the $n$ coins and obtain $k$ heads, then the probability we will win is $p_k$. The probability of obtaining $k$ heads is $\binom{n}{k}\frac{1}{2^n}$. Therefore,
\[
p_n = \sum_{k=0}^n \frac{1}{2^n} \binom{n}{k} p_k\qquad\text{for }n=2,3,\dots
\]Note that $p_n$ appears on both sides of the equation, since it’s possible that we flip all Heads, and we remain with $n$ coins at the next turn. We can solve for $p_n$ by rearranging the equation, and we obtain the complete recursion:

$\displaystyle
\begin{aligned}
p_0&=1,\quad p_1=1,\\
p_n &= \frac{1}{2^n-1} \sum_{k=0}^{n-1}\binom{n}{k} p_k\qquad\text{for }n=2,3,\dots
\end{aligned}
$

Let’s start with a simple simulation to see what we’re dealing with. Here is what we get for $n$ up to 50:

Seems like we’re done! The sequences appears to settle quickly. A numerical evaluation shows the limit to be approximately 0.72135, or 72.1%.

Down the rabbit hole

Appearances can be deceiving! First, we can expect the number of coins to roughly drop by a factor of 2 at each turn, so we should expect that $p_n \approx p_{2n}$ for large $n$. This suggests that we should plot $n$ on a log scale. Let’s see what happens if we make the same plot as before, except this time go up to $2^{12} = 4{,}096$, and we ignore the initial large swings in values to just focus on the tail end.

Would you look at that! The plot tends to a sinusoid. And although the amplitude decays rapidly at first, it does not decay to zero! If you have been following the Riddler for some time, this may seem familiar to you. In fact, a very similar pattern of decaying sinusoids on a log scale that don’t quite decay to zero was observed in the puzzle of The lonesome king, published about five and a half years ago!

It might be possible to actually compute all the way up to $10^6$, but it would likely take quite awhile (on my laptop, anyway). Computing up to $n=5000$ already took over one minute. And the scaling is not linear! Computing each $p_n$ involves all previous $p_k$, so the computation scales like $n^2$. Moreover, the binomial coefficients become so large that computing them exactly would take a lot of time and memory, and approximating them would further complicate the solution.

To approximate the probability when $n=10^6$, the easiest way is to leverage the periodicity of the function and reduce the quantity of interest to one that is within the range of values we have already computed. We have that:
\[
p_{10^6} = p_{10^6/2^8} \approx p_{3096} \approx 0.72134685.
\]Granted, the limiting amplitude of the sinusoid is so small that even if we use our original guess of 0.72135, we will be correct to 5 decimal places for any $n$ sufficiently large. So at the end of the day, we don’t gain much by doing all this extra work.

As with the Lonesome King puzzle, we are once again faced with this puzzling phenomenon of marginally stable oscillations. There are several questions that seem natural to ask:

Is there a closed-form expression for the limiting amplitude?
Is there a closed-form expression for the limiting mean?
Is the limiting waveform actually a true sinusoid?

To answer these questions, we have to go deeper…

Deeper down the rabbit hole

After more thought (and many great comments here and on Twitter), I decided to take another crack at the problem. I started with the approach used by Allen Gu, where he derived a non-recursive formula for $p_n$:
\[
p_n = \frac{n}{2}\sum_{t=0}^\infty 2^{-t} \left( 1-2^{-t} \right)^{n-1}
\]Define the inside function as
\[
g_n(t) = \frac{n}{2}2^{-t} \left( 1-2^{-t} \right)^{n-1}
\]Here is what $g_n(t)$ looks like for different values of $n$:

It appears that $g_n(t)$ gets shifted to the right as we increase $n$. The height and width of the peak appears to remain constant. To confirm this mathematically, we can solve $\frac{\mathrm{d}}{\mathrm{d}t}g_n(t) = 0$ and we find that the peak occurs at $t = \log_2(n)$. Since our plan is to sum $g_n(t)$ over all $t$ anyway, then we might as well shift it so that it is centered at its peak. If the offset from the peak is $\tau$, then the shifted function is:
\begin{align}
h_n(\tau) &= g_n(\log_2(n) + \tau) \\
&= 2^{-\tau-1} \left(1-\frac{2^{-\tau}}{n}\right)^{n-1}
\end{align}We can evaluate the limit as $n\to\infty$, and we see that:
\[
h(\tau) = \lim_{n\to\infty}h_n(\tau) = 2^{-\tau-1} e^{-2^{-\tau}}
\]It follows that in the limit of large $n$, the probability $p_n$ will take on all values in the range
\[
p(x) = \sum_{\tau=-\infty}^\infty 2^{-\tau-1+x} e^{-2^{-\tau+x}}\quad\text{for }x\in[0,1]
\]Of course, due to the translation invariance of the sum, this must be a periodic function. Here is what it looks like:

It looks very close to a sinusoid. However, it is not. If it were a true sinusoid, then $p(x)+p(x+\tfrac{1}{2})$ would be constant for all $x$. We can verify that this is not the case by doing an accurate numerical evaluation of the sums for two different values of $x$. We find:
\begin{align}
\tfrac{1}{2}(p(0.3)+p(0.8)) &\approx 0.72134752045108504298 \\
\tfrac{1}{2}(p(0.5)+p(1.0)) &\approx 0.72134752043912493759
\end{align}So they differ in the $11^\text{th}$ decimal digit.

Given that the function $p(x)$ is not a true sinusoid, we have to re-examine what we mean by “limiting value”. One way to define this would be the average value over one period. Thankfully, this integral can be evaluated exactly, and we obtain:
\begin{align}
\int_0^1 p(x)\,\mathrm{d}x &= \int_0^1 \sum_{\tau=-\infty}^\infty 2^{-\tau+x-1} e^{-2^{-\tau+x}}\,\mathrm{d}x \\
&= \sum_{\tau=-\infty}^\infty \int_0^1 2^{-\tau+x-1} e^{-2^{-\tau+x}}\,\mathrm{d}x \\
&= \sum_{\tau=-\infty}^\infty \frac{e^{2^{-\tau}}-e^{-2^{-(\tau-1)}}}{2\log (2)} \\
&= \lim_{M,N\to\infty} \sum_{\tau=-N}^M \frac{e^{2^{-\tau}}-e^{-2^{-(\tau-1)}}}{2\log (2)} \\
&= \frac{1}{2\log(2)} \left[ \lim_{M\to\infty} e^{-2^{-M}}-\lim_{N\to\infty}e^{-2^{N+1}} \right]\\
&= \frac{1}{2\log(2)}
\end{align}The bi-infinite sum is actually a telescoping series, which is pretty nice!

In conclusion, we showed that the limiting function is periodic (peak-to-peak is about $1.4\times 10^{-5}$), it is not a true sinusoid (but very close), and its mean value is exactly $\frac{1}{2\log(2)}$.

Rabbits all the way down

Our previous expression for the limiting distribution was:
\[
p(x) = \sum_{\tau=-\infty}^\infty 2^{-\tau+x-1} e^{-2^{-\tau+x}}\quad\text{for }x\in[0,1]
\]This is a periodic function, so why not find its Fourier series? To refresh your memory, since $p(x)$ is periodic with period $1$, we can write
\[
p(x) = \sum_{n=-\infty}^\infty c_n e^{2\pi i n x},\quad\text{where:}\quad
c_n = \int_0^1 p(x) e^{-2\pi i n x}\,\mathrm{d}x
\]Amazingly, we can evaluate $c_n$ in (somewhat) closed form:
\begin{align}
c_n &= \int_0^1 p(x) e^{-2\pi i n x} \\
&= \int_0^1 \sum_{\tau=-\infty}^\infty 2^{-\tau+x-1} e^{-2^{-\tau+x}} e^{-2\pi i n x}\,\mathrm{d}x \\
&= \sum_{\tau=-\infty}^\infty \int_0^1 2^{-\tau+x-1} e^{-2^{-\tau+x}} e^{-2\pi i n x}\,\mathrm{d}x
\end{align}Substitute $u = 2^{-\tau+x}$, so $\mathrm{d}u = u\log(2) \mathrm{d}x$, and $x=\tau+\log_2(u)$:
\begin{align}
c_n &= \sum_{\tau=-\infty}^\infty \frac{1}{2\log(2)}\int_{2^{-\tau}}^{2^{-\tau+1}} e^{-u}e^{-2\pi in(\tau+\log_2(u))}\,\mathrm{d}u \\
&= \frac{1}{2\log(2)} \sum_{\tau=-\infty}^\infty \int_{2^{-\tau}}^{2^{-\tau+1}} e^{-u}u^{-2\pi in/\log(2)}\,\mathrm{d}u \\
&= \frac{1}{2\log(2)} \int_{0}^{\infty} e^{-u}u^{-2\pi i n/\log(2)}\,\mathrm{d}u \\
&= \frac{\Gamma\bigl(1-\tfrac{2\pi i n}{\log(2)}\bigr)}{2\log(2)},
\end{align}where $\Gamma(z)$ is the gamma function. Therefore, the Fourier series of $p(x)$ is
\[
p(x) = \sum_{n=-\infty}^\infty \frac{\Gamma\bigl(1-\tfrac{2\pi i n}{\log(2)}\bigr)}{2\log(2)} e^{2\pi i n x}
\]Define $\alpha_n := \Gamma\bigl(1-\tfrac{2\pi i n}{\log(2)}\bigr) = |\alpha_n|e^{i \arg(\alpha_n)}$, and we can write the series in a more familiar way using standard trig functions:
\[
p(x) = \frac{1}{2\log(2)} + \frac{1}{\log(2)}\sum_{n=1}^\infty |\alpha_n| \cos\bigl( 2\pi n x+\arg(\alpha_n) \bigr)
\]There is actually a closed-form expression for the magnitude of the gamma function when the real part of the argument is an integer (see here). Specifically, we have $|\Gamma(1+ib)|^2 = \frac{\pi b}{\sinh(\pi b)}$ for all $b\in \mathbb{R}$. Therefore, we can write the Fourier series for $p(x)$ as:
\[
p(x) = \frac{1}{2\log(2)} + \sum_{n=1}^\infty \sqrt{ \frac{2 \pi^2 n \,\text{csch}\bigl(\tfrac{2 \pi ^2 n}{\log (2)}\bigr)}{\log^{3}(2)}
} \cos\bigl( 2\pi n x + \arg(\alpha_n) \bigr)
\]where $\text{csch}(z) := \frac{2}{e^z-e^{-z}}$ is the hyperbolic cosecant. The Fourier coefficients $|\alpha_n|$ decrease very rapidly with increasing $n$. In fact, we have:
\begin{align}
c_n &= \sqrt{ \frac{2 \pi^2 n \,\text{csch}\bigl(\tfrac{2 \pi ^2 n}{\log (2)}\bigr)}{\log^{3}(2)}} \\
c_1 &\approx 7.1301\times 10^{-6} \\
c_2 &\approx 6.60339\times 10^{-12} \\
c_3 &\approx 5.29624\times 10^{-18} \\
\end{align}As we increase $n$, we have $c_n \sim \exp\bigl(-\tfrac{n \pi^2}{\log(2)}\bigr) \approx 10^{-6.184 n}$. This is why the amplitude of $p(x)$ is so small, and why it looks so close to being a sinusoid when we plot it.

Using this information, we can return to the original problem and write out an exact expression for the steady-state function in terms of the original $n$ (number of coins). This yields the expression:
\[
p_n \to \frac{1}{2\log(2)} + \sum_{k=1}^\infty \sqrt{ \frac{2 \pi^2 k \,\text{csch}\bigl(\tfrac{2 \pi ^2 k}{\log (2)}\bigr)}{\log^{3}(2)}
} \cos\Bigl( 2\pi k \log_2(n) + \arg \Gamma\bigl(1-\tfrac{2\pi i k}{\log(2)}\bigr) \Bigr)
\]

The first few terms of the series are:
\begin{multline}
p_n \approx 0.72134752044448170368 \\
+ \left(7.130117\times 10^{-6}\right) \cos\bigl( 2\pi \log_2(n) +0.872711\bigr) \\
+ \left(6.603386\times 10^{-12}\right) \cos\bigl( 4\pi \log_2(n) +2.51702\bigr) + \dots
\end{multline}

Finally, here is the numerical data overlaid with the first harmonic approximation (only keeping the first term of the sum). Looks pretty good!

Equipped with this series, we can now evaluate the probability of finding a “luckiest coin” when there are $n=10^6$ coins. Here is a very accurate approximation (all digits significant)
\[
p_{1,000,000} \approx 0.7213539630726652117823954768
\]

Vehicular trouble

This week’s Riddler Classic is about steady-state mixing of fluids. Here is the paraphrased problem.

Your old van holds 12 quarts of transmission fluid. At the moment, all 12 quarts are “old.” But changing all 12 quarts at once carries a risk of transmission failure. Instead, you decide to replace the fluid a little bit at a time. Each month, you remove one quart of old fluid, add one quart of fresh fluid and then drive the van to thoroughly mix up the fluid. Unfortunately, after precisely one year of use, what was once fresh transmission fluid officially turns “old.” You keep up this process for many, many years. One day, immediately after replacing a quart of fluid, you decide to check your transmission. What percent of the fluid is old?

Here is my solution:
[Show Solution]

To completely describe the state of the fluid, we need to keep track of how much fluid is $k$ months old for $k=0,1,\dots.$. Once a parcel of fluid reaches 12 months old, we can simply consider it “old”, and we can stop counting months. So we need 13 numbers in total to do the job. Let $x_t \in \mathbb{R}^{13}$ be the state of the fluid at month $t$. Initially, all 12 quarts of fluid are old, so we can write:
\[
x_0 = \begin{bmatrix}0 \\ 0 \\ \vdots \\ 0 \\ 12\end{bmatrix}
\]During the course of the month, the fluid ages, so we shift all components of the vector down by one. We can represent this as matrix multiplication:
\[
x_{t}^{\text{shift}} = \underbrace{\begin{bmatrix}0 & 0 & 0 & \cdots & 0 & 0 & 0\\
1 & 0 & 0 & \cdots & 0 & 0 & 0\\
0 & 1 & 0 & \cdots & 0 & 0 & 0\\
\vdots & \ddots & \ddots & \ddots & \vdots & \vdots & \vdots\\[-3mm]
0 & 0 & \ddots & \ddots & 0 & 0 & 0\\[-3mm]
0 & 0 & 0 & \ddots & 1 & 0 & 0\\
0 & 0 & 0 & \cdots & 0 & 1 & 1\end{bmatrix}}_S x_t.
\]At the start of the next month, we remove one total quart of fluid, and add one quart of new fluid. The key here is that since the fluid is evenly mixed when we remove the quart, we remove the same proportion of each different type of fluid. In other words, we lose $\tfrac{11}{12}$ of all fluid types. Then, we add one quart of fluid, all of which is new. We can represent this as the linear combination:
\[
x_{t+1} = \frac{11}{12} x_t^{\text{shift}} + \underbrace{\begin{bmatrix}1 \\ 0 \\ \vdots \\ 0\end{bmatrix}}_{e_0}
\]Finally, the fraction of old fluid at time $t$ is found by extracting the last component, again representable by a matrix multiplication, and dividing it by 12 (the total amount of fluid):
\[
y_t = \underbrace{\begin{bmatrix}0 & 0 & \cdots & 0 & 1\end{bmatrix}}_{e_{12}^\textsf{T}} x_t
\]Putting all these facts together and eliminating $x_t^\text{shift}$, we get the following (state-space) representation for the dynamics:
\begin{align}
x_{t+1} &= \tfrac{11}{12}S x_t + e_0 \\
y_t &= \tfrac{1}{12}e_{12}^\textsf{T} x_t
\end{align}To find the steady-state distribution (as $t\to\infty$), we can solve for the fixed point $x_\star$ by setting $x_t=x_{t+1}$. This leads to the equations:
\begin{align}
x_\star &= \tfrac{11}{12}S x_\star + e_0 \\
y_\star &= \tfrac{1}{12}e_{12}^\textsf{T} x_\star
\end{align}Eliminating $x_\star$ leads to:

$\displaystyle
y_\star = \tfrac{1}{12}e_{12}^\textsf{T} \left( I-\tfrac{11}{12}S \right)^{-1} e_0
= \frac{3138428376721}{8916100448256} \approx 0.352
$

So after many months have passed, approximately 35.2% of the fluid will be old. We can also compute the steady-state distribution $x_\star$ directly via $x_\star = \left( I-\tfrac{11}{12}S \right)^{-1} e_0$, and we obtain the following distribution:

As anticipated, there should always be exactly one quart of new fluid. The amount in each subsequent bin decreases until we reach the old fluid, which accounts for 35.2% of the total amount of fluid (about 4.22 quarts).

Using the final value theorem

There are many ways to solve this problem (Markov chains come to mind). Since this is a problem concerning the steady-state response of a dynamical system, I thought I would present an interpretation of the problem using the language of control theory.

We can write the transfer function of this system, which turns out to be:
\[
Y(z) = \tfrac{1}{12}e_{12}^\textsf{T} \left( zI-\tfrac{11}{12}S \right)^{-1} e_0 =
\frac{3138428376721}{8916100448256 z^{12} (12 z-11)}
\]We are being asked what the steady-state output is to a constant input of 1. This is equivalent to asking for the steady-state response to a unit step. For this, we can apply the final value theorem, which tells us that the answer must be:
\[
\lim_{t\to\infty}\, y(t) = \lim_{z\to 1}\, (z-1) \cdot \frac{1}{z-1} \cdot Y(z) = Y(1) = \frac{3138428376721}{8916100448256}
\]which is the same as the answer found using the other approach.

But wait there’s more…

This problem exhibits a particularly interesting behavior. Unlike most linear systems that reach their steady-state asymptotically, this system does so in finite time! Specifically, we reach the steady state after exactly 12 months. There is no need to wait any longer than that!

To see why this is the case, return to the dynamics:
\begin{align}
x_\star &= \tfrac{11}{12}S x_\star + e_0 \\
y_\star &= \tfrac{1}{12}e_{12}^\textsf{T} x_\star
\end{align}We now make an observation about $S$:
\[
S^{12} = \begin{bmatrix}0 & 0 & \cdots & 0 & 0 \\
0 & 0 & \cdots & 0 & 0 \\
\vdots & \vdots & \cdots & \vdots & \vdots \\
0 & 0 & \cdots & 0 & 0 \\
1 & 1 & \cdots & 1 & 1 \end{bmatrix} = e_{12} 1^\mathsf{T}
\]This makes sense; after 12 months, all fluid has become old. In fact, $S^k = e_{12} 1^\mathsf{T}$ for all $k \geq 12$. Now consider the error, which is the deviation from the steady-state $x_\star$. The error dynamics are:
\[
(x_{t+1}-x_\star) = \tfrac{11}{12}S (x_0-x_\star)
\]After $k$ steps, where $k\geq 12$, the error is therefore:
\begin{align}
x_k-x_\star &= \left( \tfrac{11}{12}S \right)^k (x_0-x_\star) \\
&= \left( \tfrac{11}{12} \right)^k e_{12} \underbrace{1^{\textsf{T}} (x_0-x_\star)}_{=0}
\end{align}The reason this last term cancels is because the total sum of fluid of all ages is always 12. So the sum of the components of $x_0$ is equal to the sum of components of $x_{\star}$. So the fluid actually reaches its steady state after exactly 12 months! We can verify this using the following Matlab code:

S=diag(ones(12,1),-1); S(13,13)=1;
e0=[1;zeros(12,1)];
e12=[zeros(12,1);1];
Y=chgTimeUnit(ss(11/12*S,e0,e12',0,-1),'months');
lsim(Y,ones(16,1),0:15,12*e12)
ylabel('old fluid (quarts)')
xticks(0:2:15)
yticks(0:12)
grid

Which produces the following plot:

As we can see, after the 12th refill, we have exactly reached the steady state, and from that point forward, the amount of old fluid remains constant.

Can you eat all the chocolates?

This week’s Riddler Classic is a neat puzzle about eating chocolates.

I have 10 chocolates in a bag: Two are milk chocolate, while the other eight are dark chocolate. One at a time, I randomly pull chocolates from the bag and eat them — that is, until I pick a chocolate of the other kind. When I get to the other type of chocolate, I put it back in the bag and start drawing again with the remaining chocolates. I keep going until I have eaten all 10 chocolates.

For example, if I first pull out a dark chocolate, I will eat it. (I’ll always eat the first chocolate I pull out.) If I pull out a second dark chocolate, I will eat that as well. If the third one is milk chocolate, I will not eat it (yet), and instead place it back in the bag. Then I will start again, eating the first chocolate I pull out.

What are the chances that the last chocolate I eat is milk chocolate?

Here is my original solution:
[Show Solution]

Let’s solve the more general version of the problem, where we initially have $M$ milk chocolates and $N$ dark chocolates in the bag. When referring to a state of the game, I will write a pair $(m,n)$, where $m$ and $n$ denote the number of milk chocolates and dark chocolates remaining, respectively. Let’s define a few functions given that we are currently in state $(m,n)$:

$F(m,n)$: probability of winning given the state was just reset.
$H_1(m,n)$: probability of winning if we just ate a milk chocolate.
$H_2(m,n)$: probability of winning if we just ate a dark chocolate.

Note that $F(m,0)=H_1(m,0)=H_2(m,0) = 1$ for all $m\geq 1$, since we always win if there are only milk chocolates remaining. Also, we have that $F(0,n)=H_1(0,n)=H_2(0,n) = 0$ for all $n \geq 1$, since we always lose if there are only dark chocolates remaining.

We can write several recursive relationships relating these functions. First, once the game starts, we always eat the first chocolate we pick, and this determines our next state of the game:
\[
F(m,n) = \tfrac{m}{m+n} H_1(m-1,n) + \tfrac{n}{m+n} H_2(m,n-1)
\]Next, once we eat a chocolate, we will either continue to do so, or we will reset the state of the game, depending on which chocolate we pick next:
\begin{align}
H_1(m,n) &= \tfrac{m}{m+n}H_1(m-1,n) + \tfrac{n}{m+n}F(m,n) \\
H_2(m,n) &= \tfrac{m}{m+n}F(m,n) + \tfrac{n}{m+n}H_2(m,n-1)
\end{align}Let’s try to eliminate $H_1$ and $H_2$ in order to obtain a recursion in $F$ alone. Starting with the expression for $H_1$, we can recurse it backward to obtain:
\begin{align}
&H_1(m,n)\\
&= \tfrac{n}{m+n}F(m,n) + \tfrac{m}{m+n}H_1(m-1,n) \\
&= \tfrac{n}{m+n}F(m,n) + \tfrac{m}{m+1}\left( \tfrac{n}{m+n-1}F(m-1,n) + \tfrac{m-1}{m+n-1}H_1(m-2,n) \right) \\
&= \cdots
\end{align}Aggregating the terms and using the definition of the binomial coefficient, we can write the expression for $H_1$ compactly as
\[
H_1(m,n) = \frac{1}{\binom{m+n}{m}}\sum_{i=1}^m \binom{n+i-1}{n-1}F(i,n)
\]Similarly, we can recurse the expression for $H_2$ and obtain:
\[
H_2(m,n) = \frac{1}{\binom{m+n}{m}}\left[ 1 + \sum_{j=1}^n \binom{m+j-1}{m-1}F(m,j)\right]
\]The reason the expressions for $H_1$ and $H_2$ look slightly different in spite of all the symmetry in this problem is due to the different initial conditions; recall that $H_1(m,0)=H_2(m,0) = 1$, whereas $H_1(0,n)=H_2(0,n) = 0$.

Now substitute these expressions into our recursion for $F$ and after simplifying, we obtain:
\begin{align}
F(m,n) &= \tfrac{m}{m+n} H_1(m-1,n) + \tfrac{n}{m+n} H_2(m,n-1)\\
&= \frac{1}{\binom{m+n}{m}}\left[ 1 + \sum_{i=1}^{m-1} \binom{n+i-1}{n-1}F(i,n) + \sum_{j=1}^{n-1} \binom{m+j-1}{m-1}F(m,j)\right]
\end{align}Together with the initial conditions $F(i,0)=1$ and $F(0,j)=0$, this recursion gives us everything we need to find $F(m,n)$ for all $m$ and $n$.

Evaluating the first few terms, we discover something shocking… $F(m,n)$ appears to be exactly equal to $\tfrac{1}{2}$ for all values of $m$ and $n$! We can verify this by mathematical induction. Since our plan is to apply the above equation to generate all values of $F$, let’s assume that $F(i,j)=\tfrac{1}{2}$ for all terms on the right-hand side of the equation. Then we can compute:
\begin{align}
F(m,n) &= \frac{1}{\binom{m+n}{m}}\left[ 1 + \frac{1}{2}\sum_{i=1}^{m-1} \binom{n+i-1}{n-1} + \frac{1}{2}\sum_{j=1}^{n-1} \binom{m+j-1}{m-1}\right] \\
&= \frac{1}{\binom{m+n}{m}}\left[ 1 + \frac{1}{2}\left( \binom{m+n-1}{n}-1\right) + \frac{1}{2}\left( \binom{m+n-1}{m}-1 \right)\right]\\
&= \frac{1}{2\binom{m+n}{m}}\left[\binom{m+n-1}{n} + \binom{m+n-1}{m}\right] \\
&= \frac{1}{2\binom{m+n}{m}}\left[\binom{m+n-1}{n} + \binom{m+n-1}{n-1}\right] \\
&= \frac{1}{2\binom{m+n}{m}}\binom{m+n}{m} \\
&= \frac{1}{2}
\end{align}In the first step, we used the Hockey Stick Identity, in the third step we used the symmetry identity $\binom{n}{k}=\binom{n}{n-k}$, and in the fourth step we used the Pascal’s Triangle identity, namely $\binom{n}{k} + \binom{n}{k-1} = \binom{n+1}{k}$.

So amazingly, the probability of winning the game is independent of how many chocolates of each type are in the bag to begin with (as long as you have at least one of each type to begin with)! I have yet to come up with an elegant or intuitive way to understand why the answer should be independent of the initial number of chocolates, so if you have any ideas, leave them in the comments section below!

More types of chocolate?

Given the surprising solution to the previous problem, one might suspect that the pattern continues; if there are three different types of chocolates in the bag (say, milk chocolate, dark chocolate, and white chocolate), then perhaps the probability of eating a milk chocolate last is $\tfrac{1}{3}$? Unfortunately, this is not the case. Here are the recursions:
\begin{align}
F(m,n,p) &= \tfrac{m}{m+n+p} H_1(m-1,n,p) + \tfrac{n}{m+n+p} H_2(m,n-1,p) + \tfrac{p}{m+n+p} H_3(m,n,p-1) \\
H_1(m,n,p) &= \tfrac{m}{m+n+p}H_1(m-1,n,p) + \tfrac{n+p}{m+n+p}F(m,n,p) \\
H_2(m,n,p) &= \tfrac{n}{m+n+p}H_2(m,n-1,p) + \tfrac{m+p}{m+n+p}F(m,n,p) \\
H_3(m,n,p) &= \tfrac{p}{m+n+p}H_3(m,n,p-1) + \tfrac{m+n}{m+n+p}F(m,n,p)
\end{align}The initial conditions are a bit more complicated. For all $m\geq 1$:
\[
F(m,0,0)=H_1(m,0,0)=H_2(m,0,0)=H_3(m,0,0)=1
\]and for all $n\geq 0$ and $p \geq 0$ with $n+p \geq 1$, we have:
\[
F(0,n,p)=H_1(0,n,p)=H_2(0,n,p)=H_3(0,n,p)=0
\]A quick numerical check reveals that the solution is not $F(m,n,p)=\tfrac{1}{3}$. Here are a few values for $i=1,2,\dots$ and $j=1,2,\dots$:
\begin{align}
F(i,j,1) &= \begin{bmatrix}
\frac{1}{3} & \frac{13}{36} & \frac{3}{8} & \frac{23}{60} & \frac{7}{18} & \frac{11}{28} \\
\frac{5}{18} & \frac{7}{24} & \frac{71}{240} & \frac{107}{360} & \frac{25}{84} & \frac{25}{84} \\
\frac{1}{4} & \frac{31}{120} & \frac{31}{120} & \frac{115}{448} & \frac{685}{2688} & \frac{2041}{8064} \\
\frac{7}{30} & \frac{43}{180} & \frac{53}{224} & \frac{209}{896} & \frac{1237}{5376} & \frac{291}{1280} \\
\frac{2}{9} & \frac{19}{84} & \frac{299}{1344} & \frac{293}{1344} & \frac{493}{2304} & \frac{21343}{101376} \\
\frac{3}{14} & \frac{73}{336} & \frac{857}{4032} & \frac{199}{960} & \frac{10271}{50688} & \frac{40283}{202752} \\
\end{bmatrix} \\
F(i,j,2) &= \begin{bmatrix}
\frac{13}{36} & \frac{5}{12} & \frac{107}{240} & \frac{167}{360} & \frac{10}{21} & \frac{163}{336} \\
\frac{7}{24} & \frac{1}{3} & \frac{127}{360} & \frac{131}{360} & \frac{187}{504} & \frac{379}{1008} \\
\frac{31}{120} & \frac{53}{180} & \frac{4159}{13440} & \frac{203}{640} & \frac{31133}{96768} & \frac{52357}{161280} \\
\frac{43}{180} & \frac{49}{180} & \frac{2551}{8960} & \frac{1171}{4032} & \frac{2957}{10080} & \frac{4757}{16128} \\
\frac{19}{84} & \frac{65}{252} & \frac{929}{3456} & \frac{5507}{20160} & \frac{83627}{304128} & \frac{65203}{236544} \\
\frac{73}{336} & \frac{125}{504} & \frac{10393}{40320} & \frac{10529}{40320} & \frac{372011}{1419264} & \frac{13285}{50688} \\
\end{bmatrix} \\
F(i,j,3) &= \begin{bmatrix}
\frac{3}{8} & \frac{107}{240} & \frac{29}{60} & \frac{227}{448} & \frac{1405}{2688} & \frac{4309}{8064} \\
\frac{71}{240} & \frac{127}{360} & \frac{2561}{6720} & \frac{3567}{8960} & \frac{39623}{96768} & \frac{67351}{161280} \\
\frac{31}{120} & \frac{4159}{13440} & \frac{1}{3} & \frac{55961}{161280} & \frac{1275}{3584} & \frac{1283833}{3548160} \\
\frac{53}{224} & \frac{2551}{8960} & \frac{24679}{80640} & \frac{128069}{403200} & \frac{288059}{887040} & \frac{233803}{709632} \\
\frac{299}{1344} & \frac{929}{3456} & \frac{517}{1792} & \frac{105989}{354816} & \frac{1297237}{4257792} & \frac{3793873}{12300288} \\
\frac{857}{4032} & \frac{10393}{40320} & \frac{490247}{1774080} & \frac{1013213}{3548160} & \frac{3575315}{12300288} & \frac{1580549}{5381376} \\
\end{bmatrix}
\end{align}As anticipated, symmetry dictates that $F(k,k,k)=\tfrac{1}{3}$, and we do observe this in the numerical results above, but aside from that there is no immediately apparent pattern.

I didn’t put too much effort into trying to find a closed form for this case or the more general case with $k$ types of chocolates. If you think of a neat way to do it, leave a message in the comments section below!

And here is a far more elegant solution, courtesy of @rahmdphd on Twitter.
[Show Solution]

Delirious ducks

This week’s Riddler Classic is about random walks on a lattice:

Two delirious ducks are having a difficult time finding each other in their pond. The pond happens to contain a 3×3 grid of rocks.

Every minute, each duck randomly swims, independently of the other duck, from one rock to a neighboring rock in the 3×3 grid — up, down, left or right, but not diagonally. So if a duck is at the middle rock, it will next swim to one of the four side rocks with probability 1/4. From a side rock, it will swim to one of the two adjacent corner rocks or back to the middle rock, each with probability 1/3. And from a corner rock, it will swim to one of the two adjacent side rocks with probability 1/2.

If the ducks both start at the middle rock, then on average, how long will it take until they’re at the same rock again? (Of course, there’s a 1/4 chance that they’ll swim in the same direction after the first minute, in which case it would only take one minute for them to be at the same rock again. But it could take much longer, if they happen to keep missing each other.)

Extra credit: What if there are three or more ducks? If they all start in the middle rock, on average, how long will it take until they are all at the same rock again?

Here is my solution:
[Show Solution]

One duck

We will solve this problem by starting with a simpler one and building up to the full problem. Let’s consider a single duck on the grid of rocks. At each step, the duck moves randomly to one of the other nearby rocks. We can describe this process using a Markov chain. Each state is a different rock, and the directed edges are labeled with the transition probabilities. Here is a diagram of what the Markov chain looks like:

If we number the states in the same way as above, then the transition matrix for this Markov chain is:
\[
P = \begin{bmatrix}
0 & \tfrac12 & 0 & \tfrac12 & 0 & 0 & 0 & 0 & 0 \\
\tfrac13 & 0 & \tfrac13 & 0 & \tfrac13 & 0 & 0 & 0 & 0 \\
0 & \tfrac12 & 0 & 0 & 0 & \tfrac12 & 0 & 0 & 0 \\
\tfrac13 & 0 & 0 & 0 & \tfrac13 & 0 & \tfrac13 & 0 & 0 \\
0 & \tfrac14 & 0 & \tfrac14 & 0 & \tfrac14 & 0 & \tfrac14 & 0 \\
0 & 0 & \tfrac13 & 0 & \tfrac13 & 0 & 0 & 0 & \tfrac13 \\
0 & 0 & 0 & \tfrac12 & 0 & 0 & 0 & \tfrac12 & 0 \\
0 & 0 & 0 & 0 & \tfrac13 & 0 & \tfrac13 & 0 & \tfrac13 \\
0 & 0 & 0 & 0 & 0 & \tfrac12 & 0 & \tfrac12 & 0 \\
\end{bmatrix}
\]The way to interpret this matrix is that $P_{ij}$ is the probability that we will transition from $i$ to $j$. This explains why the rows sum to 1 (the matrix is right-stochastic). We can express this fact neatly using matrix multiplication as: $P \mathbf{1} = \mathbf{0}$, where $\mathbf{0}$ and $\mathbf{1}$ are column vectors of all-zeros and all-ones, repectively. Mathematically, we can propagate probability distributions through this Markov chain via matrix multiplication. If we have some initial distribution $\mathbf{a}$ across the states, for example
\[
\mathbf{a}^\mathsf{T} = \begin{bmatrix} \tfrac12 & \tfrac12 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}
\]i.e. we are equally likely to be in states 1 or 2. Then at the next step, the probability distribution will transition $\mathbf{a}\to\mathbf{b}$ with
\[
\mathbf{b}^\mathsf{T} = \mathbf{a}^\mathsf{T}P = \begin{bmatrix} \tfrac16 & \tfrac14 & \tfrac16 & \tfrac14 & \tfrac16 & 0 & 0 & 0 & 0 \end{bmatrix}
\]So in the next step, we might be in states 1, 2, 3, 4, 5 with the associated probabilities above.

It’s important to distinguish states e.g. $\{1,2,\dots,9\}$ from distributions over states, which are the vectors of probabilities such as the $\mathbf{a}$ and $\mathbf{b}$ used above. In the case where a distribution across states is degenerate, i.e. concentrated entirely on a particular state $s$, then I’ll use the following notation to denote the corresponding distribution over states:
\[
\left( \mathbf{e}_s\right)_i = \begin{cases} 1 & \text{if } i=s \\ 0 & \text{otherwise} \end{cases}
\]So, for example, $\mathbf{e}_2^\mathsf{T} = \begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}$.

Stopping time

We are interested in the expected hitting time, which is the number of steps it will take on average to travel from some initial state $s$ to some target state $t \in \mathcal{T}$. To keep this general, I’m assuming there can be multiple target states, indicated by the set $\mathcal{T}$. Let’s define $\mathbf{q}_s$ to be the expected hitting time to get from state $s$ to any state in our terminal set $\mathcal{T}$. It turns out that $\mathbf{q}$ satisfies the recursive relationship:
\[
\mathbf{q}_s = \begin{cases}
0 & \text{if }s \in \mathcal{T} \\
1 + \sum_{i} P_{si} \mathbf{q}_i & \text{otherwise}
\end{cases}
\]The first case is clear: if we start in the terminal set, then we have already arrived, so the hitting time is zero. If we are outside the terminal set, then the expected hitting time will be $1$ plus the weighted sum of the hitting times of wherever we end up after the next transition. The recursion above is essentially the Bellman equation from dynamic programming.

Define the vector $\mathbf{t}_i = \begin{cases} 0 & \text{if }i \in \mathcal{T}\\ 1 & \text{otherwise} \end{cases}$.
We can rewrite the equation above in concise vector form as:
\[
\mathbf{q} = \textrm{diag}(\mathbf{t}) \left( \mathbf{1} + P \mathbf{q} \right)
\]Using the fact that $\textrm{diag}(\mathbf{t})\mathbf{1} = \mathbf{t}$, we can further simplify and obtain:
\[
\left( I-\textrm{diag}(\mathbf{t}) P \right) \mathbf{q} = \mathbf{t}
\]For the case of one duck, if we set the terminal set to be $\mathcal{T} = \{5\}$, then we substitute $\mathbf{t} = \begin{bmatrix} 1 & 1 & 1 & 1 & 0 & 1 & 1 & 1 & 1 \end{bmatrix}^\mathsf{T}$ above and obtain:
\[
\begin{bmatrix}
1 & -\tfrac12 & 0 & -\tfrac12 & 0 & 0 & 0 & 0 & 0 \\
-\tfrac13 & 1 & -\tfrac13 & 0 & -\tfrac13 & 0 & 0 & 0 & 0 \\
0 & -\tfrac12 & 1 & 0 & 0 & -\tfrac12 & 0 & 0 & 0 \\
-\tfrac13 & 0 & 0 & 1 & -\tfrac13 & 0 & -\tfrac13 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & -\tfrac13 & 0 & -\tfrac13 & 1 & 0 & 0 & -\tfrac13 \\
0 & 0 & 0 & -\tfrac12 & 0 & 0 & 1 & -\tfrac12 & 0 \\
0 & 0 & 0 & 0 & -\tfrac13 & 0 & -\tfrac13 & 1 & -\tfrac13 \\
0 & 0 & 0 & 0 & 0 & -\tfrac12 & 0 & -\tfrac12 & 1 \\
\end{bmatrix}
\begin{bmatrix}\mathbf{q}_1\vphantom{\tfrac12} \\ \mathbf{q}_2\vphantom{\tfrac12} \\ \mathbf{q}_3\vphantom{\tfrac12} \\
\mathbf{q}_4\vphantom{\tfrac12} \\ \mathbf{q}_5\vphantom{\tfrac12} \\ \mathbf{q}_6\vphantom{\tfrac12} \\
\mathbf{q}_7\vphantom{\tfrac12} \\ \mathbf{q}_8\vphantom{\tfrac12} \\ \mathbf{q}_9\vphantom{\tfrac12} \end{bmatrix}
= \begin{bmatrix}
1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\
0\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12}
\end{bmatrix}
\]Inverting this matrix and solving for $\mathbf{q}$, here is the result in array form:
\[
\begin{array}{|c|c|c|} \hline
\mathbf{q}_1 = 6 & \mathbf{q}_2 = 5 & \mathbf{q}_3 = 6 \\ \hline
\mathbf{q}_4 = 5 & \mathbf{q}_5 = 0 & \mathbf{q}_6 = 5 \\ \hline
\mathbf{q}_7 = 6 & \mathbf{q}_8 = 5 & \mathbf{q}_9 = 6 \\ \hline
\end{array}
\]So if we start at node 1, it will take us on average 6 moves to reach node 5. It is by no means obvious that the expected hitting times should be integers, since they represent an average path length over an infinite number of possible paths. They just turn out to be integers in this case. The symmetry observed here also makes sense, and we could leverage it to simplify the problem to a Markov chain with only 3 states, but I’ll leave that as an exercise to the reader! (the reason I didn’t do this from the start is that I want to generalize easily to the multi-duck case, and the symmetries there are more complicated).

What about the variant where we start at the origin but we want to know how much time it will take us to return to the origin? In this case, we simply start counting from our second move, and we add 1 to the answer to account for the move we skipped. If we start at node 5, i.e. our initial distribution is $\mathbf{e}_5^\mathsf{T}$, then our distribution on the second turn is $\mathbf{e}_5^\mathsf{T} P$, and so the expected number of moves to start at 5 and return to 5 is: $\mathbf{e}_5^\mathsf{T} P \mathbf{q} + 1 = 6$.

Two ducks

At first glance, the two-duck version might seem much more challenging than the on-duck version. As we shall see, the problem is essentially the same once we look at it the right way. The key is to imagine a Markov chain where instead of the states being $\{1,2,\dots,9\}$, the states are $\{(1,1), (1,2), \dots, (9,9)\}$. In other words, there are 81 states, consisting of all possible pairs $(s_1,s_2)$, where $s_i \in \{1,2,\dots,9\}$ is the position of duck $i$.

What would be the transition matrix for such a Markov chain? The transition probability $(a_1, a_2) \to (b_1, b_2)$ is simply $P_{a_1,b_1} P_{a_2,b_2}$, i.e. the product of the respective transition probabilities for each duck. This means that if we order our states according to the first duck first (lexicographical ordering), then the transition matrix for this augmented Markov chain will be $P\otimes P$, where $\otimes$ is the Kronecker product.

So we can solve the problem in a very similar fashion. This time, our terminal set consists of all the nodes $(s_1, s_2)$ where $s_1=s_2$. There are 9 such nodes, and we can form the associated $\mathbf{t}$ by evaluating $\mathbf{t} = \mathbf{1}-\textrm{vec}(I)$, where $\mathrm{vec}$ is the vectorized version of the matrix, which you obtain by enumerating the columns. For example:
\[
\mathrm{vec} \begin{bmatrix}a_{11} & a_{12} & a_{13} \\
a_{21} & a_{22} & a_{23} \\
a_{31} & a_{32} & a_{33}
\end{bmatrix} =
\begin{bmatrix}
a_{11} \\ a_{21} \\ a_{31} \\ a_{12} \\ a_{22} \\ a_{32} \\ a_{13} \\ a_{23} \\ a_{33} \end{bmatrix}
\]the new equation for $\mathcal{q}$ is $\bigl( I-\textrm{diag}(\mathbf{t}) (P\otimes P) \bigr) \mathbf{q} = \mathbf{t}$ and the solution to the problem where you start at $(5,5)$ and want to know the expect time it takes for us to reach the terminal set (where the two ducks are on the same rock) can be computed following the same logic as in the one-duck case. That is, $\mathbf{e}_{(5,5)}^\mathsf{T} (P\otimes P) \mathbf{q} + 1$.

A cautionary note

While it may seem like this is straightforward algebra, there is a major issue I neglected to discuss: invertibility. When we solve the equation for $\mathbf{q}$ in the one-duck case, the equation looks like $A\mathbf{q} = \mathbf{t}$, and the solution is $\mathbf{q} = A^{-1}\mathbf{t}$. But what if $A$ isn’t invertible? This corresponds to cases where the Markov chain is not connected. So if you start on one island and the terminal set is on another island, you’ll never get there; the stopping time is infinite. That doesn’t happen with one duck, but it does happen with two ducks! This is because every time a duck moves to a new stone, the parity of the stone changes from odd to even or vice versa. So if one duck started on rock 1 and the other started on rock 2, they would never meet.

From a linear algebra standpoint, it just means that even though the matrix is not invertible, if we restrict our attention to the “island” (read: subspace) where the two ducks have the same parity, then the associated transition matrix for that island will be invertible and everything will be fine. From a practical standpoint, the equation above can still be solved; we just need to replace the inverse by a pseudoinverse, and ignore all the components of $\mathbf{q}$ for which the two components have different parity. More on this later…

Computation

I computed the result using the following Matlab code:

% transition matrix for one duck
P = [ 0 1/2 0 1/2 0  0  0  0  0
     1/3 0 1/3 0 1/3 0  0  0  0
      0 1/2 0  0  0 1/2 0  0  0
     1/3 0  0  0 1/3 0 1/3 0  0
      0 1/4 0 1/4 0 1/4 0 1/4 0
      0  0 1/3 0 1/3 0  0  0 1/3
      0  0  0 1/2 0  0  0 1/2 0
      0  0  0  0 1/3 0 1/3 0 1/3
      0  0  0  0  0 1/2 0 1/2 0 ];

T = kron(P,P);     % transition matrix
t = 1-vec(eye(9)); % terminal states

% starting distribution
s = kron([0 0 0 0 1 0 0 0 0],[0 0 0 0 1 0 0 0 0])';

% compute stopping time as a rational number
stop_time = s'*T*( (eye(81)-diag(t)*T)\t ) + 1
rats(stop_time)

The resulting expected stopping time was $\frac{363}{74}$, or approximately $4.905$ steps. We can also find an approximate answer via simulation. Here is what you get by simulating one million trials:

The staggered decrease of the probabilities is a real effect; it’s not a result of approximation error!

Many ducks

If we have $n$ ducks, we can clearly just generalize the approach used above. This time, the transition matrix will be $\underbrace{P\otimes \cdots \otimes P}_{n\text{ times}}$. The issue here is that our transition matrix will be quite large: $9^n \times 9^n$, to be precise. One way to shrink the number of states is to return to the idea of parity. The reason our transition matrix is large is because it describes all transitions between all possible duck configurations. Since we know all ducks start on the same node and that parity will be preserved, we only need to worry about the subset of states for which all ducks have the same parity. Namely, there will be $5^n$ odd states and $4^n$ even states. So we can restrict ourselves to a smaller $(5^n+4^n)\times(5^n+4^n)$ transition matrix. This is still large, but it’s an improvement. We will also leverage the fact that the transition matrix will be sparse, which will help with computation.

After all this work, we end up with the following result (for up to 6 ducks):

The values are as follows:

Number of ducks	Expected rendez-vous time
2	4.9054
3	18.4360
4	66.7420
5	237.3955
6	825.3364

It’s clear the expected rendez-vous time grows exponentially, but I haven’t found a practical way to approximate it or bound it. Unfortunately, it’s pretty difficult to extend my approach beyond $n=6$ ducks. To give you an idea of scale, our naive transition matrix would have had $9^6=531,441$ rows and columns. After our reduction procedure, we’re down to $5^6+4^6=19,721$. The transition matrix is sparse (about 98.5% of its entries are zeros), but that still leaves about 6 million nonzero entries, which makes computing $\mathbf{q}$ a challenge.

Further reductions are possible if we leverage symmetries. Let’s say we have 2 ducks. If we don’t do any reductions, there are $9^2=81$ states. If we reduce using even/odd parity, we drop to $5^2+4^2 = 41$.

From a stopping time standpoint, the only thing that matters is how many ducks are on each stone; it doesn’t matter which ducks are on which stones. This immediately allows us to reduce $5^n+4^n$ to ${n+4 \choose 4} + {n+3\choose 3}$, which is a dramatic improvement. In the two-duck case, this brings us to $25$.
Deeper symmetries occur as well. For example, a pair of ducks on 1 and 3 is equivalent to a pair of ducks on 7 and 9 by rotational symmetry. In the two-duck case, this brings us to $8$ states. Counting and bookkeeping these states might require some group theory…

Ultimately, these reductions would have a dramatic effect, reducing the scaling from exponential to polynomial. I didn’t implement any of these additional reductions because at this point (up to $n=6$) the exponential trend was clear. I didn’t think that extending the trend further would be particularly enlightening.

Beer pong

This interesting twist on the game of Beer Pong appeared on the Riddler blog. Here it goes:

The balls are numbered 1 through N. There is also a group of N cups, labeled 1 through N, each of which can hold an unlimited number of ping-pong balls. The game is played in rounds. A round is composed of two phases: throwing and pruning.

During the throwing phase, the player takes balls randomly, one at a time, from the infinite supply and tosses them at the cups. The throwing phase is over when every cup contains at least one ping-pong ball. Next comes the pruning phase. During this phase the player goes through all the balls in each cup and removes any ball whose number does not match the containing cup. Every ball drawn has a uniformly random number, every ball lands in a uniformly random cup, and every throw lands in some cup. The game is over when, after a round is completed, there are no empty cups.

How many rounds would you expect to need to play to finish this game? How many balls would you expect to need to draw and throw to finish this game?

Here is my solution:
[Show Solution]

We’ll address the second question first: How many balls would you expect to need to draw and throw to finish the game? For this question, the rounds don’t matter at all. The game will end once there is at least one correct ball in each cup, and the pruning phase at the end of each round only removes incorrect balls so we can effectively neglect it.

Each time a ball is thrown, there is a $1/N$ chance it will land in the correct cup. This is a classical coupon collector problem (see my write-up on the card collection completion problem for more background). The solution rests on the fact that if some trial is successful with probability $p$, then we must perform $1/p$ trials on average before we obtain our first success. The probability of landing our first correct ball is $1/N$ so it will take an average of $N$ tosses to do so. The probability of landing our next new correct ball (in a different cup) is $\tfrac{1}{N}\cdot \tfrac{N-1}{N}$. This takes on average $\tfrac{N^2}{N-1}$ more tosses. Continuing in this fashion and summing up all of the intervals until each cup is filled, the total expected number of tosses is:

$\displaystyle
T_\text{tot} = N^2 \sum_{k=1}^N \frac{1}{k} \approx N^2(\log N + \gamma)
$

where $\gamma \approx 0.5772$ is the Euler-Mascheroni constant.

Now let’s turn our attention to the first question: How many rounds would you expect to play to finish this game? This is a much more complicated question, so to illustrate the approach I’ll solve it for the case $N=3$. We’ll begin by breaking each round into a number of stages, which track the progress of cups as they are filled. In the diagram below, each cup is either empty (no ball), is red (incorrect ball), or is green (correct ball). We can also figure out all the transition probabilities and write them on the arrows.

In this diagram, self-loops have whatever probability is required so that the outgoing arrows sum to $1$. This is an example of a Markov chain. We start on the left (all cups empty), and we work our way to the right (all cups full). Here is what the transition matrix looks like for $N=3$ when we arrange the states from top to bottom and left to right:
\[
A =
\begin{bmatrix}
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
2/3 & 2/9 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
1/3 & 1/9 & 1/3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 4/9 & 0 & 4/9 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 2/9 & 4/9 & 2/9 & 5/9 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 2/9 & 0 & 1/9 & 2/3 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 2/9 & 0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1/9 & 2/9 & 0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 1/9 & 2/9 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 1/9 & 0 & 0 & 0 & 1
\end{bmatrix}
\begin{array}{c}
(k=0,j=0)\\
(k=1,j=0)\\
(k=1,j=1)\\
(k=2,j=0)\\
(k=2,j=1)\\
(k=2,j=2)\\
(k=3,j=0)\\
(k=3,j=1)\\
(k=3,j=2)\\
(k=3,j=3)\\
\end{array}
\]Note that the columns sum to $1$ since all outgoing probabilities must sum to $1$. The question is: what is the probability that we’ll end up in each of the possible ending point? (the four last states). This is called the limiting distribution. In our case, we can compute it easily. First, we can eliminate the self-loops by normalizing. This works because we don’t actually care how much time is spent at each stage. See below for an illustration of this normalization.

Here is what our transition matrix looks like after normalization
\[
A_\text{norm} =
\begin{bmatrix}
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
2/3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
1/3 & 1/7 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 4/7 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 2/7 & 2/3 & 2/5 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 1/3 & 0 & 1/4 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 2/5 & 0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1/5 & 1/2 & 0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 1/4 & 2/3 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 1/3 & 0 & 0 & 0 & 1
\end{bmatrix}
\begin{array}{c}
(k=0,j=0)\\
(k=1,j=0)\\
(k=1,j=1)\\
(k=2,j=0)\\
(k=2,j=1)\\
(k=2,j=2)\\
(k=3,j=0)\\
(k=3,j=1)\\
(k=3,j=2)\\
(k=3,j=3)\\
\end{array}
\]Once again, the columns sum to $1$, but this time only the absorbing states have self-loops. Finding the limiting distribution simply amounts to finding the limit $\lim_{t\to\infty} A_\text{norm}^t$. Since the upper-left block of $A_\text{norm}$ is nilpotent (the diagonal is zero), powers of the matrix will eventually be constant. It suffices to evaluate $A_{\text{norm}}^{2N}$. Extracting the last N+1 rows and the columns corresponding to the states $\{1,3,6,\dots,\tfrac{(N+1)(N+2)}{2}\}$ (the states where we have exactly $k$ correct balls in cups, for $k=0,1,\dots,N$), we obtain the transition matrix for one round of the game. Here is the result we obtain in the case $N=3$:
\[
P = \begin{bmatrix}
16/105 & 0 & 0 & 0 \\
41/105 & 1/3 & 0 & 0 \\
5/14 & 1/2 & 2/3 & 0 \\
1/10 & 1/6 & 1/3 & 1
\end{bmatrix}
\begin{array}{c}
(k=0)\\
(k=1)\\
(k=2)\\
(k=3)\\
\end{array}
\]In other words, $P_{ij}$ is the probability that if we currently have $j$ correct balls in cups, we will have $i$ after the next round is over. From here, our next task is to compute the expected number of rounds it will take to travel from $k=0$ to $k=N$. This can be computed as follows: if we call $E_k$ the expected number of rounds starting from $k$, then we have:
\begin{align}
E_0 &= 1 + P_{00} E_0 + P_{10}E_1 + \cdots + P_{N0}E_N \\
E_1 &= 1 + P_{01} E_0 + P_{11}E_1 + \cdots + P_{N1}E_N \\
&\vdots \\
E_{N-1} &= 1 + P_{0,N-1} E_0 + P_{1,N-1}E_1 + \cdots + P_{N,N-1} E_N \\
E_N &= 0
\end{align}Or, in matrix-vector notation: $E = \mathbb{1} – e_N + P^\mathsf{T} E$. This can be solved via simple matrix inversion. We can then extract $E_0$ from the solution, which is the expected number of rounds starting with all empty cups. As far as I can tell, there does not appear to be an easy way to extract an analytic formula, but all the steps above are very easy to compute numerically. Here are plots of the results for small values of $N$:

I also overlayed a result obtained via Monte-Carlo simulation using $1000$ simulated games. This shows that the analytical results agree well with simulation. Here are further numerical results, this time extended up to $N=100$:

It appears that the number of rounds required grows approximately linearly with $N$ while the number of balls required grows quadratically. This latter result is to be expected, since we know the formula in this case is approximately $N^2 \log N$.

Is this bathroom occupied?

After a brief hiatus from Riddling, I’m back! This Riddler problem is about probability and bathroom vacancy.

There is a bathroom in your office building that has only one toilet. There is a small sign stuck to the outside of the door that you can slide from “Vacant” to “Occupied” so that no one else will try the door handle (theoretically) when you are inside. Unfortunately, people often forget to slide the sign to “Occupied” when entering, and they often forget to slide it to “Vacant” when exiting.

Assume that 1/3 of bathroom users don’t notice the sign upon entering or exiting. Therefore, whatever the sign reads before their visit, it still reads the same thing during and after their visit. Another 1/3 of the users notice the sign upon entering and make sure that it says “Occupied” as they enter. However, they forget to slide it to “Vacant” when they exit. The remaining 1/3 of the users are very conscientious: They make sure the sign reads “Occupied” when they enter, and then they slide it to “Vacant” when they exit. Finally, assume that the bathroom is occupied exactly half of the time, all day, every day.

Two questions about this workplace situation:

1. If you go to the bathroom and see that the sign on the door reads “Occupied,” what is the probability that the bathroom is actually occupied?
2. If the sign reads “Vacant,” what is the probability that the bathroom actually is vacant?
Extra credit: What happens as the percentage of conscientious bathroom users changes?

Here is how I solved the problem:
[Show Solution]

This problem can be solved using Bayes’ Rule but you can also use Markov chains. I like Markov chains so that’s what we’ll do!

The first step is to define what are states and transition probabilities are. This is the tricky part; we might be tempted to think that because the bathroom can be either occupied or not, and the sign in front can either read “Vacant” or “Occupied”, that there should be four states (one for each possible pair).

However, this is not the case. The reason is because the transition probabilities are not independent of one another. Consider the state “bathroom is occupied and the sign says it’s occupied”. We must distinguish between the cases where the person occupying the bathroom is conscientious (they will definitely slide the sign to “Vacant” when they leave) or not (they will leave the sign as “Occupied” after they leave).

We must therefore add additional states to our Markov chain that correspond to the different ways in which the bathroom can be occupied. To keep things general, let’s assume the fraction of oblivious, forgetful, and conscientious users are $p$, $q$, and $r$, respectively. Therefore $0\le p,q,r \le 1$ and $p+q+r=1$. A diagram of the Markov chain:

If $\{x_k\}$ denotes the probability of being in state $k$ in the diagram, the stationary distribution satisfies the equation:

\[
\begin{bmatrix}
0 & 0 & 1 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 1 \\
p & 0 & 0 & 0 & 0 & 0 \\
r & r & 0 & 0 & 0 & 0 \\
q & q & 0 & 0 & 0 & 0 \\
0 & p & 0 & 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix}
x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \\ x_6
\end{bmatrix}
=
\begin{bmatrix}
x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \\ x_6
\end{bmatrix}
\]

Solving these equations, we obtain:

\begin{align}
x_1 &= \frac{r}{2 (q+r)}, & x_2 &= \frac{q}{2 (q+r)}, & x_3 &= \frac{p r}{2 (q+r)}, \\
x_4 &= \frac{r}{2}, & x_5 &= \frac{q}{2}, & x_6 &= \frac{pq}{2 (q+r)}
\end{align}

We can now compute the probabilities we’re after:

\begin{align}
\mathbb{P}(\text{Vacant}\mid \text{says “Vacant”})
&= \frac{x_1}{x_1+x_3} = \frac{1}{1+p}\\
\mathbb{P}(\text{Occupied}\mid \text{says “Occupied”})
&= \frac{x_4+x_5+x_6}{x_2+x_4+x_5+x_6} = \frac{r+q-pr}{r+2q-pr}
\end{align}

When $p=q=r=\frac{1}{3}$ (equal fraction of conscientious, forgetful, and oblivious users), we obtain:
\begin{align}
\mathbb{P}(V\mid V) &= \frac{3}{4}, &
\mathbb{P}(O\mid O) &= \frac{5}{8}
\end{align}

Bonus: different probabilities

We have the analytic expressions above for what happens as a function of $p,q,r$… so now it remains to visualize it somehow. One possibility is to set two of them equal and just make a plot. For example, if we only vary the percentage of conscientious users ($r$) and set the other two to be equal, i.e. $p=q=\frac{1-r}{2}$. We then obtain:
\begin{align}
\mathbb{P}(V\mid V) &= \frac{2}{3-r}, &
\mathbb{P}(O\mid O) &= \frac{r^2+1}{r^2-r+2}
\end{align}We can also plot these probabilities as a function of $r$:

If we want to see what happens as we vary the percentages of $p,q,r$ all at once, we can visualize the result using a ternary density plot. Here, we plot values in a triangle; the fraction of each type of bathroom user is determined by the length of the altitude to the opposite side of the triangle and the color indicates the probability. Here are the plots:

You can read off the appropriate percentages by looking at the tick marks. For example, the lower-left corner corresponds to everybody being conscientious, the bottom (horizontal) side of the triangle corresponds to nobody being oblivious, and as you move from left to right you increase the fraction of forgetful users while simultaneously decreasing the fraction of conscientious users. The center of the triangle corresponds to $p=q=r=\tfrac{1}{3}$, and so on.

One additional comment: you’ll notice the point at the top of the triangle (100% oblivious users) is missing. This is because this probability doesn’t make sense. When all users are oblivious, the sign will never change. So if it starts at “Vacant”, it will stay at “Vacant” forever.

Colorful balls puzzle

This Riddler puzzle about an interesting game involving picking colored balls out of a box. How long will the game last?

You play a game with four balls: One ball is red, one is blue, one is green and one is yellow. They are placed in a box. You draw a ball out of the box at random and note its color. Without replacing the first ball, you draw a second ball and then paint it to match the color of the first. Replace both balls, and repeat the process. The game ends when all four balls have become the same color. What is the expected number of turns to finish the game?

Extra credit: What if there are more balls and more colors?

Here is my solution to the first part (four balls):
[Show Solution]

Before I begin, I’d like to acknowledge Hector Pefo and Sawyer Tabony who also posted excellent solutions to this Riddler puzzle. We all arrived at the same answer (whew!) but our approaches differ slightly.

A natural way to think about this problem is that the game can exist in some number of states (the particular coloring of the balls in the box). At each step, at most one of the balls is recolored and we transition to another state. Such a collection of states and transition probabilities is known as a Markov Chain. Since we are only interested in how long it takes for all balls to be the same color, it’s not necessary to keep track of every possible state; we can aggregate states and simplify the problem.

Partition approach

One way to aggregate states is to count the number of balls of each different color without regard for the colors themselves. For example, 4 balls can be partitioned in 5 different ways:
\[
1+1+1+1,\qquad
2+1+1,\qquad
2+2,\qquad
3+1,\qquad
4
\]The partition “2+1+1”, for example, consists of cases where two of the balls are the same color and the other two balls are two other colors. Both the cases “red+red+green+blue” and “blue+blue+yellow+red” would fall into the category “2+1+1”. By using these five partitions as states in a Markov Chain, we can compute the transition probabilities to go from one state to the next. For example, the probability of transitioning from “2+1+1” to “3+1” is $\frac{1}{3}$ because in order for this transition to occur, we must first choose one of the two identically colored balls (probability $\frac{2}{4}$), then we must choose one of the other two balls out of the remaining three (probability $\frac{2}{3}$).

Here is what the complete Markov Chain looks like when all transition probabilities are filled in:

If we label the states in order, we can write the transition probabilities as a matrix that describes how states evolve.
\[
A = \begin{bmatrix}
0 & 0 & 0 & 0 & 0 \\
1 & \frac{1}{2} & 0 & 0 & 0 \\
0 & \frac{1}{6} & \frac{1}{3} & \frac{1}{4} & 0 \\
0 & \frac{1}{3} & \frac{2}{3} & \frac{1}{2} & 0 \\
0 & 0 & 0 & \frac{1}{4} & 0
\end{bmatrix}
\]For example, if $x$ is our current distribution, then $Ax$ is our distribution after one move. Note that the columns sum to $1$ because at each state, we must transition to some other state. The only exception is the last column, because the game stops once we reach the final state.

To compute the probability that the game ends after $k$ turns, we should find the probability that the game is in the final state after $k$ turns. Since we start in state $1$, the required quantity is $\left[ A^k \right]_{51}$. In other words, the $(5,1)$ component of $A^k$. We can compute this directly for each $k$ and we obtain the following distribution:

Note that the probability is zero for 1 or 2 turns because the game can’t end that quickly. It takes at least three turns (since there are four balls total) for all balls to be painted the same color. We can also compute the expected number of turns by evaluating the sum:
\[
A + 2A^2 + 3A^3 + \dots
= A(I-A)^{-2}
\]and then extracting the $(5,1)$ component. The result in this case turns out to be 9. When we compare the distribution to the mean in the figure above, we notice that the distribution has a heavy tail; although the mean is 9, the likeliest number of turns is actually 5.

Here is my solution to the general case with $N$ balls:
[Show Solution]

Although it’s possible to extend the partition used in the previous part to any other number of balls, in practice it is a challenge. One of the first obstacles is that the number of partitions as a function of $N$ is an intractable quantity. It is known as the partition function and is studied by number theorists for its deep connections to prime numbers. Partition functions were even mentioned in the recent movie about Ramanujan, “The Man Who Knew Infinity”. So we’re not taking that approach.

Counting Blue

We can still use Markov Chains, but we’ll have to come up with a different way of aggregating the states. One possibility is to define the states $k = 1,2,\dots,N$ which correspond to the number of Blue balls in the box. If Blue ends up winning, we will have started in state 1 and made our way eventually to state N without ever getting to state 0. The trick is to realize that no matter which color wins, the winning color has a Markov Chain like this one. All colors are indistinguishable so we might as well consider the case where Blue wins.

The next step is to compute transition probabilities. Here, we must be careful because we are conditioning on Blue winning. So we use Bayes’ Rule to compute the transitions. Given that there are $k$ Blue balls with $1 \le k \le N-1$, the probability that we transition $k\to k+1$ is:
\[
A_{k+1,k} = \frac{\mathbb{P}(k\to k+1 \text{, Blue wins})}{\mathbb{P}(\text{Blue wins})}
= \frac{\frac{k}{N}\cdot \frac{N-k}{N-1} \cdot \frac{k+1}{N}}{\frac{k}{N}} = \frac{(k+1)(N-k)}{N(N-1)}
\]Similarly, the probability that we transition $k\to k-1$ is:
\[
A_{k-1,k} = \frac{\mathbb{P}(k\to k-1 \text{, Blue wins})}{\mathbb{P}(\text{Blue wins})}
= \frac{\frac{N-k}{N}\cdot \frac{k}{N-1} \cdot \frac{k-1}{N}}{\frac{k}{N}} = \frac{(k-1)(N-k)}{N(N-1)}
\]The only other possibility is that we don’t transition at all, which is what happens the rest of the time. So $A_{k,k} = 1-A_{k+1,k}-A_{k-1,k}$. Simplifying the expressions, we obtain:
\[
A_{k+1,k} = \frac{(k+1)(N-k)}{N(N-1)},\quad
A_{k-1,k} = \frac{(k-1)(N-k)}{N(N-1)},\quad
A_{k,k} = 1-\frac{2k(N-k)}{N(N-1)}
\]As before, the $k^\text{th}$ column sums to 1. The only exception, as before, is the last column $A_{:,N}$, which is zero since the game ends once all colors are the same. The approach from here is identical to what we did when we used partitions. We can compute the probability of ending in $k$ turns by computing $\left[ A^k \right]_{N,1}$ and the expected number of turns is given by $\left[ A(I-A)^{-2} \right]_{N,1}$. Here are the distributions for $N=5$ and $N=8$.

As we can see, the distribution looks similar as we increase $N$. Curiously, the mean always appears to be $(N-1)^2$. In the next section, we will show that this is indeed the case in general!

Analytic expressions

Writing out the transition matrix explicitly, we have:
\[
A = \begin{bmatrix}
1-\frac{2(N-1)}{N(N-1)} & \frac{1(N-2)}{N(N-1)} & & & & 0 \\
\frac{2(N-1)}{N(N-1)} & 1-\frac{4(N-2)}{N(N-1)} & \frac{2(N-3)}{N(N-1)} & & & 0 \\
& \frac{3(N-2)}{N(N-1)} & 1-\frac{6(N-3)}{N(N-1)} & \ddots & & 0 \\
& & \frac{4(N-3)}{N(N-1)} & \ddots & \frac{(N-2)1}{N(N-1)} & 0 \\
& & & \ddots & 1-\frac{2(N-1)1}{N(N-1)} & 0 \\
& & & & \frac{N\cdot 1}{N(N-1)} & 0
\end{bmatrix}
\]If we define $a_k$ to be the expected number of turns until the game ends given that we are currently at state $k$, we can write the recursion:
\begin{align}
a_k &= 1 + A_{k+1,k} a_{k+1} + A_{k,k} a_k + A_{k-1,k} a_{k-1} \quad\text{for }k=1,\dots,N-1 \\
a_N &= 0
\end{align}Rearranging this expression, it looks roughly like $A^Ta+1=a$ (with the exception of the last row). Writing this as a single compact equation, we can factor things a bit and obtain a simpler equation:
\begin{multline}
\frac{1}{N(N-1)}
\begin{bmatrix}N-1 & & &\\ & \ddots & & \\ & & 2 & \\ & & & 1\end{bmatrix}
\begin{bmatrix}2 & -1 & & \\ -1 & 2 & \ddots & \\ & \ddots & \ddots & -1 \\ & & -1 & 2\end{bmatrix}\\
\times \begin{bmatrix}1 & & \\ & 2 & & \\ & & \ddots & \\ & & & N-1\end{bmatrix}
\begin{bmatrix}a_1 \\ a_2 \\ \vdots \\ a_{N-1}\end{bmatrix}
=
\begin{bmatrix}1 \\ 1 \\ \vdots \\ 1\end{bmatrix}
\end{multline}Or, simplified:
\[
\begin{bmatrix}2 & -1 & & \\ -1 & 2 & \ddots & \\ & \ddots & \ddots & -1 \\ & & -1 & 2\end{bmatrix}
\begin{bmatrix}a_1 \\ 2a_2 \\ \vdots \\ (N-1)a_{N-1}\end{bmatrix}
=
N(N-1)\begin{bmatrix}\frac{1}{N-1} \\ \frac{1}{N-2} \\ \vdots \\ 1\end{bmatrix}
\]By performing an LU decomposition, we may instead solve the following system of equations in sequence:
\begin{align}
\begin{bmatrix} 1 & & & & \\ -\frac{1}{2} & 1 & & & \\ & -\frac{2}{3} & 1 & & \\ & & \ddots & \ddots & \\ & & & -\frac{N-2}{N-1} & 1\end{bmatrix}
\begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_{N-2} \\ y_{N-1} \end{bmatrix}
&= N(N-1)\begin{bmatrix}\frac{1}{N-1} \\ \frac{1}{N-2} \\ \vdots \\ \frac{1}{2} \\ 1\end{bmatrix}\\
\begin{bmatrix} 2 & -1 & & & \\ & \frac{3}{2} & -1 & & \\ & & \frac{4}{3} & \ddots & \\ & & & \ddots & -1 \\ & & & & \frac{N}{N-1} \end{bmatrix}
\begin{bmatrix} a_1 \\ 2a_2 \\ \vdots \\ (N-2)a_{N-2} \\ (N-1)a_{N-1} \end{bmatrix}
&= \begin{bmatrix}y_1 \\ y_2 \\ \vdots \\ y_{N-2} \\ y_{N-1}\end{bmatrix}
\end{align}The first system is easily solved for $y$ by forward-substitution. We first have $y_1=N$. Then, $y_2 = \frac{N(N-1)}{N-2}+\frac{1}{2}y_1$, and so on. The result we find is:
\[
y_k = N(N-1)\sum_{i=1}^k \frac{i}{(N-i)k}
\]The second system can be solved for $a$ by back-subsitution, starting from the final equation and working our way back up. Doing this, we find:
\[
a_k = \sum_{j=1}^{N-k} \frac{1}{k+j} y_{k+j-1}
\]We may now combine both expressions and obtain:
\[
a_k = N(N-1) \sum_{j=1}^{N-k} \frac{1}{k+j} \sum_{i=1}^{k+j-1} \frac{i}{(N-i)(k+j-1)}
\]This is a bit messy, but we only care about $k=1$ (the expected number of turns given that we start with only one blue ball and blue ends up winning). This sum can be simplified with a bit of work:
\begin{align}
a_1 &= N(N-1) \sum_{j=1}^{N-1} \frac{1}{j(j+1)} \sum_{i=1}^{j} \frac{i}{(N-i)} \\
&= N(N-1) \sum_{j=1}^{N-1} \left( \frac{1}{j}-\frac{1}{j+1}\right) \sum_{i=1}^{j} \frac{i}{(N-i)}
\end{align}Define $Q_j := \frac{1}{j+1}\sum_{i=1}^j \frac{i}{N-i}$ and the sum telescopes:
\begin{align}
a_1 &= N(N-1)\sum_{j=1}^{N-1} \left( \frac{1}{N-j} + Q_{j-1}-Q_j\right)\\
&= N(N-1)\left( \sum_{j=1}^{N-1}\frac{1}{N-j}-Q_{N-1} \right) \\
&= N(N-1)\left( \sum_{j=1}^{N-1}\frac{1}{N-j}-\frac{1}{N}\sum_{j=1}^{N-1}\frac{j}{N-j} \right) \\
&= N(N-1) \sum_{j=1}^{N-1}\frac{1}{N} \\
&= (N-1)^2
\end{align}And this is precisely what we were trying to show! The expected number of turns for a game with $N$ balls is $(N-1)^2$. A similar expansion can be done to compute $a_k$ for $k\ne 1$ but in this case there is only partial telescoping and we are left with
\[
a_k = (N-1)^2-\frac{N(N-1)}{k}\sum_{i=1}^k \frac{k-i}{N-i}
\]As $N$ gets large, the sum of reciprocals is well-approximated by a logarithm. If we denote $\alpha = k/N$ the fraction of balls that are blue, then the expected number of turns until all balls are blue is:
\[
\frac{a_k}{(N-1)^2} \approx \left(\frac{1-\alpha}{\alpha}\right)\log\left(\frac{1}{1-\alpha}\right)
\]Here is a scaled plot of this distribution of expected number of turns:

A conjecture?

I was intrigued by the distribution of the number of turns rather than just the expected value. As seen in the plots above, the distribution converges to something, but it’s not clear what. And there doesn’t appear to be an easy way to compute or approximate powers of $A$ analytically. I suspect the distribution tends to a log-normal distribution as $N$ gets large because when plotted on a log scale, it sure looks close to normal:

If you have any thoughts, leave a comment below!

A supreme court puzzle

This timely Riddler puzzle is about filling supreme court vacancies…

Imagine that U.S. Supreme Court nominees are only confirmed if the same party holds the presidency and the Senate. What is the expected number of vacancies on the bench in the long run?

You can assume the following:

You start with an empty, nine-person bench.
There are two parties, and each has a 50 percent chance of winning the presidency and a 50 percent chance of winning the Senate in each election.
The outcomes of Senate elections and presidential elections are independent.
The length of time for which a justice serves is uniformly distributed between zero and 40 years.

Here is my solution:
[Show Solution]

The first thing to realize about this problem is that we can consider each of the nine supreme court separately. If $X_i$ is the vacancy of the $i^\text{th}$ seat ($1$ if it’s vacant, $0$ if it’s occupied), then the quantity we’re after is:
\[
\mathbb{E}\left[ \sum_{i=1}^9 X_i \right] = \sum_{i=1}^9 \mathbb{E}\left[X_i\right] = 9\, \mathbb{E}\left[ X_1 \right]
\]and this follows by linearity of expectation. So it suffices to consider the problem of a one-person bench and multiply our final result by $9$.

General solution strategy

Let’s say that the government is “aligned” if the same party holds the presidency and the senate. Imagine the distribution of $X_1$ (vacancy of the first seat) over time. When the seat becomes occupied, $X_1 = 0$, we must wait some number of years that is uniformly distributed between zero and $40$ years before it becomes vacant again. When the seat become vacant, $X_1 = 0$, two things can happen:

If the government is aligned when the seat becomes vacant, then it can be filled right away with no wait required. The seat will remain filled for an average duration of 20 years (the expected value of the uniform distribution on the interval $[0,40]$).
If the government is not aligned when the seat becomes vacant then the seat will remain vacant until the next election. Then, if the election aligns the government, the seat will fill. Otherwise, the wait continues.

Senate elections happen every two years while presidential elections happen every four years. Every two years, there is an independent probability of $1/2$ that the government becomes aligned. This occurs when the senate aligns itself with the presidency. It doesn’t matter whether the president changes or not since the outcome of the senate race is independent of the outcome of the presidential race. Given repeated independent events each with probability $p$, the expected number of events before we encounter our first success is $1/p$. This fact was demonstrated in the post about the Monsters’ Gems puzzle. In this case, $p=1/2$ so we can expect to wait $2$ elections before the government becomes aligned. Here is a tally of the possible cases when the seat becomes vacant:

With probability $1/2$, the government is aligned and the wait is zero.
With probability $1/2$, the government is not aligned. In this case, we must wait an average of $1$ year until the next election. Then, we must wait an average of $2$ more elections. This will only incur a wait of $2$ years on average because we must only wait the time between both elections. If the second election aligns the government, we can fill the seat right away! So the expected wait is $1+2=3$ years.

Combining these two facts, we are left with a net expected wait of $1.5$ years during which the seat is vacant.

On average, the empirical expected number of vacancies in the long run is:
\[
\frac{1.5}{20 + 1.5} = \frac{3}{43} \approx 0.069767
\]If we account for all nine seats, we multiply this number by nine:

$\displaystyle
(\text{Expected vacancies}) = \frac{27}{43} \approx 0.628
$

So the expected number of vacancies is less than one.

If we make the approximation that the variables $X_i$ are mutually independent, the expected number of vacancies is a Binomial distribution with parameters $n=9, p=\tfrac{3}{43}$. Using this fact, we can compute an approximate distribution of long-term vacancies:

We can see that the vacancies are $0$ about half the time, and are rarely greater than $2$.

Note: We made an approximation in the solution above. If a supreme court justice’s term expires while the same government is still in power (before any other elections take place), then the seat will always be filled immediately since the government is still aligned. This effect only matters for up to two years into the justice’s term. If we account for this special case, it will slightly increase the average duration of a term and therefore slightly decrease the expected rate of vacancy. Hector Pefo has a more detailed solution on his blog where he works out all the details.