Tag: linearity of expectation

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.

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!

The key feature of this sequential elimination game is that the likelihood of each possible outcome only depends on how many people are playing. So we can define the vector $x_n$ to be the probability somebody will end up winning the game, given that $n$ people are currently playing. We have the base cases $x_0=0$ and $x_1=1$, and our goal is to figure out $x_{56000}$.

At every round, some number of people are eliminated. This continues until we have $0$ or $1$ people left. Define $B(m,n)$ to be the number of ways of having $m$ people left for the next round, given that we currently have $n$ people. We will now explain how to count this.

• Everyone is eliminated: this is the quantity $B(0,n)$. In order for this to happen, each person must get picked exactly once. The number of ways this can happen is precisely the number of permutations of $\{1,\dots,n\}$. But it’s a bit more complicated, because people cannot eliminate themselves. So we must count the number of permutations with no fixed points. These are called derangements. The total number of derangements satisfies the recursion:
  \begin{align}
  B(0,0) &= 1, \quad B(0,1) = 0\\
  B(0,n+1) &= n B(0,n) + n B(0,n-1)\quad\text{for: }n=1,2,\dots
  \end{align}There is also a closed-form formula that holds for $n\ge 1$, given by $B(0,n) = \lfloor \tfrac{n!}{e} + \tfrac{1}{2} \rfloor$.

• Partial elimination: suppose we drop from $n$ people down to $m$ people with $m > 0$. This is the quantity $B(m,n)$. In this case, we must count the number of surjections from $\{1,\dots,n\}$ to itself such that (i) there are no fixed points and (ii) precisely $m$ of the elements are left untouched. We can count this recursively as well. Let $F(n,k)$ be the number of surjections from $\{1,\dots,n\}$ to $\{1,\dots,k\}$ with no fixed points. Let’s count how the number of surjections changes if we remove element $n$ from the pre-image. There are two cases: if we are left with a surjection from $\{1,\dots,n-1\}$ to $\{1,\dots,k\}$, i.e. the removed element was a redundant pairing, then it could have linked to any of the $k$ elements from the image. In other words, it contributed $k F(n-1,k)$ surjections. Alternatively, if the removed pairing was a critical link, then we will be left with a surjection from $\{1,\dots,n-1\}$ to a subset of $\{1,\dots,k\}$ with one element removed. The missing link must connect to the removed element, which can happen in $k$ ways ($k$ possible elements missing), so the net contribution is $k F(n-1,k-1)$. Consequently, we have the recursion: $F(n,k) = kF(n-1,k) + kF(n-1,k-1)$. Since we want to count all possible surjections, we can relate $F$ to $B$ via the formula: $B(m,n) = {n \choose m} F(n,n-m)$. Substituting and simplifying, we obtain a recursion only in terms of $B$’s:
  \begin{align}
  B(m,n) = \frac{n(n-m)}{m} B(m-1,n-1) + n B(m,n-1)
  \end{align}

Let $A(m,n)$ be the probability that we transition from $n$ people to $m$ people in one move. To convert $B(m,n)$ to $A(m,n)$, we divide by the total number of possible assignments. Each person can pick one of the $(n-1)$ other people, so there are $(n-1)^n$ total choices possible. Therefore, $A(m,n) = (n-1)^{-n} B(m,n)$. Substituting into the formulas above, we can obtain a recursion for the $A$’s:
\begin{align}
A(0,0) &= 1,\quad A(0,1) = 0,\quad A(1,1) = 1,\quad A(n,n) = 0,\quad n=2,3,\dots \\
A(0,n+1) &= \tfrac{(n-1)^n}{n^n} A(0,n) + \tfrac{(n-2)^{n-1}}{n^n} A(0,n-1),\quad n=1,2,\dots\\
A(m,n+1) &= \tfrac{(n-1)^n(n+1)}{n^{n+1}} \left( \tfrac{n-m+1}{m} A(m-1,n) + A(m,n) \right),\quad m=1,2,\dots,n
\end{align}Using these formulas, we can compute any of the transition probabilities we like. Here is what the $A(m,n)$ matrix looks like for $0\le m,n \le 6$:
\[
A = \small\begin{bmatrix}
1.0000& 0& 1.0000& 0.2500& 0.1111& 0.0430& 0.0170\\
0& 1.0000& 0& 0.7500& 0.5926& 0.4102& 0.2458\\
0& 0& 0& 0& 0.2963& 0.4688& 0.5069\\
0& 0& 0& 0& 0& 0.0781& 0.2150\\
0& 0& 0& 0& 0& 0& 0.0154\\
0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0
\end{bmatrix}
\]So for example, if we had 4 people playing the game, there would be a 29.63% chance 2 people get eliminated, a 59.26% chance 3 people get eliminated, and a 11.11% chance everyone gets eliminated. Note that the columns of this matrix necessarily sum to 1. We can think of the game as a Markov Chain where $A$ matrix is the transition matrix. If the current distribution of the population is $z$ (i.e. $z_k$ is the probability that we have a population of $k$), then the distribution after one round of the game is found via the matrix multiplication $Az$. We want to know what happens after infinitely many rounds have been played. In other words, we want to find the stationary distribution.

The standard way to find a stationary distribution is to solve the eigenvector equation $x^T A = x^T$. Another way is by simply computing $A^k$ for some sufficiently large $k$. Upon doing this, we can plot the limiting distribution as a function of the initial population:

There is a fascinating oscillatory behavior that occurs as the population gets large. Namely, the probability of having a winner vs. having no winner oscillates about 0.5 with a period that increases exponentially. When plotted on a log-scale as above, the stationary distribution appears to converge to a sinusoid. I couldn’t go all the way up to 56,000 but the pattern that emerged was clear, so I extrapolated out what I couldn’t compute. The best-fit curve I found (empirically) for the probability of having an eventual winner for large $n$ is:
\[
x_n \approx 0.503+0.0265 \sin(14.5 \log_{10}(n)-1.27)
\]Substituting for $n=56000$, we find there is a 47.65% chance somebody wins and a 52.35% chance nobody wins.

Note: The limiting period can be explained as follows. The probability that one person doesn’t get picked is $(1-\frac{1}{n})^{n-1}$, so by linearity of expectation, we can expect to have $n(1-\frac{1}{n})^{n-1}$ people remaining after one round. In the limit $n\to\infty$, this is equal to $n e^{-1}$. Since the population roughly shrinks by a factor of $e$ at each iteration, we can expect the probability distribution to repeat with a period of $\log_{10}(e)$ on our plot. This explains why the number inside the sinusoid in our approximate formula turned out to be $\frac{2\pi}{\log_{10}(e)} \approx 14.5$.

Aside from the observation above, I couldn’t find a complete analytic characterization of the stationary distribution, so if anybody has ideas I’d love to hear them!