Markov chains - Book Proofs

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.

The lonesome king

This Riddler puzzle is about a random elimination game. Will someone remain at the end, or will everyone be eliminated?

In the first round, every subject simultaneously chooses a random other subject on the green. (It’s possible, of course, that some subjects will be chosen by more than one other subject.) Everybody chosen is eliminated. In each successive round, the subjects who are still in contention simultaneously choose a random remaining subject, and again everybody chosen is eliminated. If there is eventually exactly one subject remaining at the end of a round, he or she wins and heads straight to the castle for fêting. However, it’s also possible that everybody could be eliminated in the last round, in which case nobody wins and the king remains alone. If the kingdom has a population of 56,000 (not including the king), is it more likely that a prince or princess will be crowned or that nobody will win? How does the answer change for a kingdom of arbitrary size?

Here is my solution:
[Show Solution]

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!