Book Proofs - A blog for mathematical riddles, puzzles, and elegant proofs

Thanos snaps

This week’s Riddler Express problem is a nod to the recently released finale to the Avengers saga. No spoilers!

Thanos, the all-powerful supervillain, can snap his fingers and destroy half of all the beings in the universe.

But what if there were 63 Thanoses, each snapping his fingers one after the other? Out of 7.5 billion people on Earth, how many can we expect would survive?

If there were N Thanoses, what would the survival fraction be?

Here is my solution:
[Show Solution]

First things first: If there are $S$ Thanos snaps and the probability of survival after one snap is $p$, then the expected fraction of survivors is $p^S$. In the case where the population is $7.5$ billion and the probability of survival is $p=0.5$, if there are $S=63$ snaps, we can calculate the expected number of survivors to be:
\[
\frac{7.5\times 10^9}{2^{63}} = 8.13\times 10^{-10}
\]This number is so small that there will almost certainly be no survivors.

A more interesting interpretation: But what if Thanos were not immune to the snap? After all, when Thanos snaps his fingers, half of all the beings in the universe are destroyed! This should include Thanos! So when the first Thanos snaps his fingers, on average half of the remaining Thanoses will be destroyed before having the opportunity to snap their fingers.

This means that $S$, the total number of snaps, is a random variable. We can compute its probability mass function recursively. Specifically, let’s define $P_N(s)$ to be the probability that there will be $S=s$ snaps given that there are initially $N$ Thanoses. Let’s start with some easy cases.

If there are no Thanoses, there will be no snaps. So $P_0(0) = 1$.
If there is at least one Thanos, there will be at least one snap. So $P_N(0) = 0$ for $N=1,2,\dots$.

The recursion for larger values of $s$ is given by:
\begin{align}
P_N(s) &= \sum_{k=0}^{N-1} \mathrm{Prob}(k\text{ Thanoses survive first snap}) P_k(s-1) \\
&= \sum_{k=0}^{N-1} \binom{N-1}{k} p^{k}(1-p)^{N-1-k} P_{k}(s-1)
\end{align}This follows because if there are $s$ snaps remaining, then during the first one, the remaining $N-1$ Thanoses are at risk of destruction. The number of surviving Thanoses follows a binomial distribution with $N-1$ trials and probability $p$, which explains why the probability that $k$ survive the first snap is equal to $\binom{N-1}{k} p^{k}(1-p)^{N-1-k}$. This recursion can be evaluated numerically. Here is what the distribution of snaps looks like for $p=0.5$ and $N=63$:

So there will most likely be $5$ or $6$ snaps. Once we know how many snaps there are, the effect on the population of the universe is the same as in the first interpretation. So the expected fraction of survivors in the universe will be $\mathbb{E}(p^S)$, or:
\[
\mathbb{E}(\text{fraction of survivors if we have }N\text{ Thanoses}) = \sum_{s=0}^N p^s P_{N}(s)
\]Here is a plot of the surviving fraction as a function of $N$ when $p=0.5$:

When there are 63 Thanoses, the expected survival rate is about 2.22%. So with an initial population of 7.5 billion, we will be left with, on average, about 166.5 million.

If you’re curious about the shape of this graph, we can approximate it like this: With $N$ Thanoses and $p=0.5$, there are roughly $\log_2(N+1)$ flips that occur. If exactly that many flips occur, the fraction of survivors will be roughly $(0.5)^{\log_2(N+1)} = \tfrac{1}{N+1}$. This tells us that as $N$ gets large, the fraction of survivors will scale roughly like $\frac{1}{N}$, where $N$ is the number of Thanoses. This is a crude approximation, but it captures the correct asymptotic behavior.

If you’re interested in the (Python) code I wrote to generate this solution (and the figures), you can find it on github here.

Circular conundrum

This week’s Riddler problem is short and sweet:

If N points are generated at random places on the perimeter of a circle, what is the probability that you can pick a diameter such that all of those points are on only one side of the newly halved circle?

Here is my solution:
[Show Solution]

Label the points $\{1,2,\dots,N\}$. Given one of the points $i$ on the perimeter, we’ll call this point a clockwise split if all remaining points lie on the semi-circle starting from point $i$ in the clockwise direction. For example, in the figure below, Point 1 is not a clockwise split because Point 4 is outside the semicircle that starts at Point 1. We can also observe that Points 2 and 3 are not clockwise splits, but Point 4 is.

If all $N$ points are located on the same half-circle, then exactly one of these points will be a clockwise split (almost surely). Meanwhile, if there exists no diameter that puts all the points on the same half-circle, then no point will be a clockwise split. Define the random variable $X_i$ as follows.
\[
X_i = \begin{cases}1 & \text{Point }i\text{ is a clockwise split} \\ 0 & \text{otherwise}\end{cases}
\]Based on the observations above, the probability that there exists a diameter such that all the points are one side of the newly halved circle (let’s call this probability $P$) is given by:
\begin{align}
P &= \text{Prob}(\text{one of the points is a clockwise split}) \\
&= \text{Expected total number of clockwise splits} \\
&= \mathbb{E}\left( \sum_{i=1}^N X_i \right) \\
&= \sum_{i=1}^N \mathbb{E}(X_i)
\end{align}By symmetry, $\mathbb{E}(X_i)$ is the same for each $i$, and it’s equal to $\frac{1}{2^{N-1}}$, since once we’ve fixed point $i$, there are $N-1$ remaining points and each one has a probability of $\tfrac{1}{2}$ of being in the correct half of the circle. Therefore, the probability we seek is:

$\displaystyle
P = \frac{N}{2^{N-1}}
$

The last step in the derivation relies on linearity of expectation. Even though the $X_i$’s are not independent (far from it, since at most one of them can be $1$), linearity of expectation still holds, and it allows us to make short work of what would otherwise be a complicated calculation.

Note: There are many ways to solve this problem. In fact, Derek commented that a previous Riddler problem is quite similar to this one, and in his comments suggested two different ways of solving the present problem, one year ago!

Gift card puzzle

Here is a puzzle from the Riddler about gift cards:

You’ve won two gift cards, each loaded with 50 free drinks from your favorite coffee shop. The cards look identical, and because you’re not one for record-keeping, you randomly pick one of the cards to pay with each time you get a drink. One day, the clerk tells you that he can’t accept the card you presented to him because it doesn’t have any drink credits left on it.

What is the probability that the other card still has free drinks on it? How many free drinks can you expect are still available?

Here is my solution:
[Show Solution]

Let’s suppose each card starts with $n$ drinks. There are three ways the sequence of buying drinks can end:

The first card gets maxed out and the second card has $k$ drinks remaining with $k\ge 1$. This means we purchased a total of $2n-k$ drinks and $n$ of them ended up on the first card. Then, we tried to purchase one additional drink on the first card and the buying stopped. The probability of this occurring is: $(\tfrac{1}{2})^{2n-k+1}\binom{2n-k}{n}$.
both cards get maxed out and then we try to buy one more drink on either card. This means we purchased a total of $2n$ drinks and $n$ of them ended up on the first card. The probability of this occurring is $(\tfrac{1}{2})^{2n}\binom{2n}{n}$.
The second card gets maxed out and the first card has $k$ drinks remaining with $k\ge 1$. This is exactly analogous to the first case, and the probability of this occurring is again: $(\tfrac{1}{2})^{2n-k+1}\binom{2n-k}{n}$.

Putting all of this together, we end up with a tidy formula for $p(n,k)$, the probability that when the game ends, the other card has exactly $k$ drinks remaining on it:

$\displaystyle
p(n,k) = \frac{1}{2^{2n-k}} \binom{2n-k}{n}
$

To double-check that this is a valid probability mass function, we should have $\sum_{k=0}^n p(n,k) = 1$. This is indeed the case, but it’s rather challenging to verify. Here is a link to the WolframAlpha computation. Here is plot of the probability mass function for $n=50$.

We can now answer the first question: the probability that the other card still has drinks on it is one minus the probability that there are no drinks left, i.e. $1-p(n,0)$. This evaluates to:

$\displaystyle
\mathrm{Prob}\Bigl(\begin{smallmatrix}\text{there are drinks left}\\\text{on the other card}\end{smallmatrix}\Bigr)
= 1-\frac{1}{4^n} \binom{2n}{n} \approx 1-\frac{1}{\sqrt{\pi n}}
$

This makes sense; the only way there can be no drinks left is if we purchased $2n$ drinks, and we were lucky to use each card exactly $n$ times. We would expect this probability to get smaller as $n$ gets larger. The approximation I used above comes from Stirling-like approximations for binomial coefficients (see here). When $n=50$ as in the problem statement, the probability evaluates to $92.04\%$. Here is a plot showing how the probability of there being drinks on the other card slowly tends to $1$ as $n$ gets larger:

Now let’s address the second question: what is the expected number of drinks on the other card? We will compute the expected value directly from the probability mass function: $E_n = \sum_{k=0}^n k\, p(n,k)$. We find:

$\displaystyle
\mathbb{E}\Bigl(\begin{smallmatrix}\text{number of drinks left}\\\text{on the other card}\end{smallmatrix}\Bigr)
= \frac{2n+1}{4^n}\binom{2n}{n}-1 \approx \frac{2n+1}{\sqrt{\pi n}}-1
$

Again, the derivation is challenging so I omitted it, but here is the WolframAlpha link. Here is a plot of this function as it changes with $n$ along with the approximation (which is quite good!)

When $n=50$ as in the problem statement, the expected number of drinks remaining on the other card is $7.0385$. This corresponds to the expected value of the probability distribution plotted in the first image.

Elf music

This holiday-themed Riddler problem is about probability:

In Santa’s workshop, elves make toys during a shift each day. On the overhead radio, Christmas music plays, with a program randomly selecting songs from a large playlist.

During any given shift, the elves hear 100 songs. A cranky elf named Cranky has taken to throwing snowballs at everyone if he hears the same song twice. This has happened during about half of the shifts. One day, a mathematically inclined elf named Mathy tires of Cranky’s sodden outbursts. So Mathy decides to use what he knows to figure out how large Santa’s playlist actually is.

Help Mathy out: How large is Santa’s playlist?

Here is my solution:
[Show Solution]

Settlers in a circle

In this Riddler problem, the goal is to spread out settlements in a circle so that they are as far apart as possible:

Antisocial settlers are building houses on a prairie that’s a perfect circle with a radius of 1 mile. Each settler wants to live as far apart from his or her nearest neighbor as possible. To accomplish that, the settlers will overcome their antisocial behavior and work together so that the average distance between each settler and his or her nearest neighbor is as large as possible.

At first, there were slated to be seven settlers. Arranging that was easy enough: One will build his house in the center of the circle, while the other six will form a regular hexagon along its circumference. Every settler will be exactly 1 mile from his nearest neighbor, so the average distance is 1 mile.

However, at the last minute, one settler cancels his move to the prairie altogether (he’s really antisocial). That leaves six settlers. Does that mean the settlers can live further away from each other than they would have if there were seven settlers? Where will the six settlers ultimately build their houses, and what’s the maximum average distance between nearest neighbors?

Here is my solution:
[Show Solution]

I approached this problem from a modeling and optimization perspective. If there are $N$ settlers and we imagine the coordinates of the settlers are $(x_i,y_i)$ for $i=1,\dots,N$ and $d_i$ is the distance between the $i^\text{th}$ settler and its nearest neighbor, then we can model the problem as follows:
\[
\begin{aligned}\underset{x,y,d\,\in\,\mathbb{R}^N}{\text{maximize}} \quad & \frac{1}{N}\sum_{i=1}^N d_i \\
\text{such that} \quad & x_i^2 + y_i^2 \le 1 &&\text{for }i=1,\dots,N\\
& (x_i-x_j)^2 + (y_i-y_j)^2 \ge d_i^2 &&\text{for }i,j=1,\dots,N\text{ and }j\ne i
\end{aligned}
\]The objective is to minimize the average distance between each settler and its nearest neighbor (the average of the $d_i$). The first constraint says that each settler must be within the circle (a distance of at most $1$ from the origin). The second constraint says that the $i^\text{th}$ settler is a distance at least $d_i$ from each of the other settlers.

The game is then to find $(x_i,y_i,d_i)$ that satisfy these constraints and maximize the objective. Broadly speaking, this is an example of a mathematical optimization problem. Unfortunately, the problem as stated is nonconvex problem because of the form of the second constraint. This means that the problem may have local maximizers that are not globally optimal. There are two ways forward:

First, we could use local search: start with randomly placed settlers, wiggle their positions in ways such that the objective continues to increase, and then stop once we can no longer improve the objective. Examples of such approaches include hill-climbing and interior-point methods. In general, such approaches only find local maxima, so we should try many random starting positions to ensure we find the best possible answer.

Second, we could use global search: these are methods that attempt to find a globally optimal solution by considering many configurations simultaneously and adding random perturbations that can “kick” a solution out of a local maximum. Examples include simulated annealing and particle swarms. These approaches tend to be much slower than local search, but they have a much better chance of finding a global optimum.

Ultimately, I used local search with several random initializations. Here is how the maximum average distance scales with the number of settlements.

The bar plot also shows how well we would do if we distributed the settlements evenly on the circumference (“all on the border”) and if we put one in the middle and the rest evenly distributed on the circumference (“one in the center”). We can also examine what the settlement distributions actually look like. Here they are:

The case $N=6$ is particularly interesting because it has a settlement that is almost in the center, but not quite. Indeed, the average minimum distance to neighbors for this case is $1.002292$, so we benefit ever so slightly by placing one settlement off-center.

We can solve for the $N=6$ case analytically as well, but the solution isn’t pretty (fair warning!) It turns out the settlements are distributed throughout the circle in the following way:

Here, $x$ is the distance between the origin and the point $A$. The other distances satisfy $z \lt y \lt x+1$, and the average distance to the nearest neighbor is therefore $\tfrac{1}{6}(1+x+2y+3z)$. Note that once we set $x$, this uniquely determines the positions of all the other points. After some messy trigonometry, we obtain the following equations relating $x$, $y$, $z$ and two auxiliary variables $p$ and $q$:
\begin{align}
y^2&=1+x-x^2-x^3\\
z^2&=1-x^2+2x\left(pq-\sqrt{1-p^2}\sqrt{1-q^2}\right) \\
p&= \tfrac{1}{2}(1-2x-x^2)\\
q&= \tfrac{1}{2}(1-x+x^2+x^3)
\end{align}We can eliminate $\{y,z,p,q\}$ and solve for the average distance $\bar d = \frac{1+x+2y+3z}{6}$ as a function of $x$ alone:
\[
\bar d = \frac{x+1}{12} \left(2+4 \sqrt{1-x}+3 \sqrt{2}\sqrt{1-x} \sqrt{\left(x^3+2 x^2-x+2\right)-x \sqrt{x+3} \sqrt{x^3+x^2-x+3}} \right)
\]We can then plot this function of $x$, and we get the following curve:

Interestingly, this function reaches its maximum just a bit after $x=0$. Zooming in, we can see more clearly:

The maximum occurs at about $x=0.0555108$ and the maximum value is $1.00229$, just as we found before. If we try to solve for this analytically by taking the derivative of this function, setting it equal to zero, and eliminating rationals, we obtain the following result… That the optimal $x$ a root of the following $30^\text{th}$ order polynomial:
\begin{align}
x^{30}&+\tfrac{152}{9} x^{29}+\tfrac{9259}{81} x^{28}+\tfrac{6851}{18} x^{27}+\tfrac{383363 }{648}x^{26}+\tfrac{1999495} {5832}x^{25}\\
&+\tfrac{43478263}{52488} x^{24}+\tfrac{75726373}{23328} x^{23}+\tfrac{3282369305}{1679616} x^{22}-\tfrac{31642768891}{7558272} x^{21}\\
&+\tfrac{367775655235}{136048896} x^{20}+\tfrac{392027905801}{34012224} x^{19}-\tfrac{2016793647847}{136048896} x^{18}-\tfrac{1258157703385}{68024448} x^{17}\\
&+\tfrac{4545130566235}{136048896} x^{16}-\tfrac{22059424877}{17006112} x^{15}-\tfrac{8448216227365}{136048896} x^{14}+\tfrac{2878062707167}{68024448} x^{13}\\
&+\tfrac{7398046956073}{136048896} x^{12}-\tfrac{312655708903}{3779136} x^{11}+\tfrac{17266856539}{1679616} x^{10}+\tfrac{63670102441}{839808} x^9\\
&-\tfrac{29827845749}{559872} x^8-\tfrac{541401565}{23328} x^7+\tfrac{619158287}{23328} x^6+\tfrac{6168305}{2592} x^5\\
&-\tfrac{156285}{32} x^4-\tfrac{5153}{36} x^3+\tfrac{3923}{12} x^2+\tfrac{45}{2} x-\tfrac{35}{16}
\end{align}Yuck! Here is what you get if you plot this polynomial; as expected, there is a root at $x=0.0555108$.

Alice and Bob fall in love

In this interesting Riddler problem, we’re dealing with a possibly unbounded sequence of… children? Here it goes:

As you may know, having one child, let alone many, is a lot of work. But Alice and Bob realized children require less of their parents’ time as they grow older. They figured out that the work involved in having a child equals one divided by the age of the child in years. (Yes, that means the work is infinite for a child right after they are born. That may be true.)

Anyhow, since having a new child is a lot of work, Alice and Bob don’t want to have another child until the total work required by all their other children is 1 or less. Suppose they have their first child at time T=0. When T=1, their only child is turns 1, so the work involved is 1, and so they have their second child. After roughly another 1.61 years, their children are roughly 1.61 and 2.61, the work required has dropped back down to 1, and so they have their third child. And so on.

(Feel free to ignore twins, deaths, the real-world inability to decide exactly when you have a child, and so on.)

Five questions: Does it make sense for Alice and Bob to have an infinite number of children? Does the time between kids increase as they have more and more kids? What can we say about when they have their Nth child — can we predict it with a formula? Does the size of their brood over time show asymptotic behavior? If so, what are its bounds?

Here is an explanation of my derivation:
[Show Solution]

Let’s call $t_k$ the time at which Alice and Bob have their $k^\text{th}$ child. The first child is born at time zero, so $t_1 = 0$. Since the work due to a child is one divided by their age, the work due to child $k$ at time $t$ is given by:
\[
\text{work at time }t = \begin{cases}
\frac{1}{t-t_k} & \text{if }t > t_k \\
0 & \text{otherwise}
\end{cases}
\]We also know that each new child is born once the total work due to all previous children equals $1$. Therefore, we have the following equations:
\begin{align}
t_2\text{ satsifies:} && \frac{1}{t_2-t_1} &= 1 \\
t_3\text{ satisfies:} && \frac{1}{t_3-t_2} + \frac{1}{t_3-t_1} &= 1 \\
t_4\text{ satisfies:} && \frac{1}{t_4-t_3} + \frac{1}{t_4-t_2} + \frac{1}{t_4-t_1} &= 1 \\
\cdots\qquad \\
t_k\text{ satisfies:} && \sum_{i=1}^{k-1} \frac{1}{t_k-t_i} &= 1
\end{align}and so on. We also know that $t_1 \lt t_2 \lt t_3 \lt \ldots$ since the children are born sequentially. The nice thing about these equation is that they can be solved sequentially. We know that $t_1=0$. Substituting this into the first equation yields $t_2=1$. Substituting these values into the second equation yields:
\[
\frac{1}{t_3-1}+\frac{1}{t_3} = 1
\]Rearranging yields $t_3^2-3t_3+1=0$. This equation has two solutions, but we take the one that satisfies $t_3 \gt t_2$, and we obtain $t_3=\tfrac{3+\sqrt{5}}{2} \approx 2.618$. We can continue in this fashion, solving each subsequent equation by substituting the previous values… but things get complicated in a hurry. The next equation becomes:
\[
t_4^3-\left(\tfrac{11+\sqrt{5}}{2}\right)t_4^2+\left(\tfrac{13+3\sqrt{5}}{2}\right)t_4-\left(\tfrac{3+\sqrt{5}}{2}\right)=0
\]Once again, we find that there is a unique solution that satisfies $t_4 \gt t_3$. This solution is approximately $t_4\approx 4.5993$. There doesn’t appear to be any automatic way to continue this process. Every time we try to compute the next $t_k$, we have to find the root of an even more complicated polynomial!

We’ll take a numerical approach, but let me first address one issue: how can we be sure that these ever-more-complicated polynomials always have a root that satisfies our constraints? For example, what if the polynomial for $t_4$ had no roots larger than $t_3$? What if it had many such roots? Thankfully, there is a simple way we can guarantee that there is always exactly one admissible root to each polynomial by leveraging the intermediate value theorem. Specifically, the function
\[
f_k(t) = \frac{1}{t-t_1} + \frac{1}{t-t_2} + \cdots + \frac{1}{t-t_{k-1}}
\]with $t \gt t_{k-1} \gt \ldots \gt t_1$ is a continuous and strictly decreasing function. Moreover, $\lim_{t\to {t_{k-1}}^+}f_k(t) = +\infty$ and $\lim_{t\to\infty} f_k(t) = 0$. Therefore, by the intermediate value theorem, there must be some point $t_k \in (t_{k-1},\infty)$ such that $f_k(t_k)=1$. The strict monoticity of $f_k$ ($f_k$ is strictly decreasing) implies that this $t_k$ must be unique.

The root may be unique, but how can we find it efficiently? Here, we can leverage a further property of $f_k$; namely convexity. The fact that $f_k$ is convex (and smooth) means that we can leverage some efficient methods from convex analysis that make use of derivative information. There are many ways to get the answer from here… What I ended up doing was to use the L-BFGS algorithm on the objective $(f_k(t)-1)^2$. My implementation in Julia using the Optim.jl package took about 40 seconds to evaluate birth dates of the first 1,000 children.

Note: there are many ways to formulate this as an optimization problem. One way is to note that solving $f_k(t)=1$ is equivalent to extremizing the antiderivative of $f_k(t)-1$. In other words, we could solve:
\[
t_k = \underset{t \gt t_{k-1}}{\text{arg min}} \left( t-\sum_{i=1}^{k-1}\log(t-t_i) \right)
= \underset{t \gt t_{k-1}}{\text{arg max}}\,\, e^{-t}\prod_{i=1}^{k-1}(t-t_i)
\]Incidentally, this last expression for $t_k$ makes it clear that there should be a unique $t_k$. When $t\gt t_{k-1}$, the product is a polynomial that is strictly increasing, and it’s multiplied by $e^{-t}$, which will drive it to zero eventually. There must therefore be a unique maximizer. I tried implementing both of these approaches and both turned out to be slower than directly minimizing $(f_k(t)-1)^2$.

If you’re just interested in the answers to the questions, here they are:
[Show Solution]

Does it make sense for Alice and Bob to have an infinite number of children? Biological limitations aside, yes! We proved that for each $t_k$, there exists a unique $t_{k+1} \gt t_k$ that satisfies the constraints of the problem. So by induction, Alice and Bob can have as many children as they like.
Does the time between kids increase as they have more and more kids? The answer is yes, and I will provide a “proof by picture”. After finding $t_k$ for $k=1,2,\dots$, we know the birth dates of each child. Then it remains to compute $(t_{k+1}-t_k)$ for each $k$ (the time between births). Here is the plot

The first point on the left says that we wait $1$ year after child $1$ is born, we wait about $1.6$ years after child $2$ is born, and so on. Plotted on a log scale as above, we see the curve is almost a perfect straight line. Here is a formula we can use to approximate the number of years to wait:

$\displaystyle
\text{years to wait after child }k\text{ is born} \,\approx\, 1 + 2.19\cdot \log_{10} k
$

So yes, the time between children increases (without bound), but it does so slowly. By the time Alice and Bob have their 1,000th child, they will be waiting approximately 7.5 years between kids.
What can we say about when they have their $N^\text{th}$ child? Is there a formula we can write down? I couldn’t find an analytic expression, but I did find a pretty good approximation. First, let’s make a plot to see what we’re working with. Here is what the birth years is like for the first 19 children

We can see here that the growth rate is slightly faster than linear. Since the time between children fit so well to a curve of the form $a+b\log x$, it makes sense to fit the birth dates themselves to the integral of this form. In an effort to fit the simplest possible curve, I did a bit of experimenting and found that the form $a x + x \log x$ does a really good job. I found the equation of best fit using linear regression and put the equation in the plot above. If we extend the time horizon, we get a similar result:

This is a pretty good fit; the plot actually depicts the true data with the fit overlaid on top. The fit is so good that you can’t even see the data underneath! So for large $N$, a good approximation is:

$\displaystyle
\text{year when child }N\text{ is born} \,\approx\,-0.315N + N\log N
$
The final question is whether the size of the brood over time shows asymptotic behavior. We already plotted the time of birth vs child number, and found the relationship to be approximately: $y \approx -0.315 N + N\log N$. So here, we need to look at the inverse relationship and ask what $N$ does as a function of $y$. Since $y$ grows a little faster than linearly, this means $N$ will grow a little more slowly than linearly.

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$.

Hand sort

A card-rearranging problem on the Riddler blog. Here it goes:

You play so many card games that you’ve developed a very specific organizational obsession. When you’re dealt your hand, you want to organize it such that the cards of a given suit are grouped together and, if possible, such that no suited groups of the same color are adjacent. (Numbers don’t matter to you.) Moreover, when you receive your randomly ordered hand, you want to achieve this organization with a single motion, moving only one adjacent block of cards to some other position in your hand, maintaining the original order of that block and other cards, except for that one move.

Suppose you’re playing pitch, in which a hand has six cards. What are the odds that you can accomplish your obsessive goal? What about for another game, where a hand has N cards, somewhere between 1 and 13?

Here is my solution:
[Show Solution]

I was unable to find an analytic or closed-form expression for the solution to this problem, so I’ll present an old fashioned brute-force solution. The idea is simple:

Enumerate all possible hands
For each hand, enumerate all possible moves that consist of moving some group of adjacent cards to somewhere else in the hand.
If any such move succeeds in sorting the hand, mark the hand as “good”.
Once we’re done, count all “good” hands and divide that by the total number of hands

Simple, right? I’d like to point out is that this is not an approximate solution; there is no simulation involved. I’m counting all the possible scenarios, so this approach produces the exact answer. The only problem is that there are a lot of hands to check, and the bigger the hands get, the longer it takes to check them, since there are more possible moves. Successfully solving this problem requires coding it up in a way that a computer can find the solution in a reasonable amount of time. In order to make this work, I used several tricks to reduce the computation required:

We don’t have to enumerate all possible hands because only the suits matter, but we do have to be careful about how we count. We begin with a full deck of cards (52 distinct cards). Let’s suppose we want to count hands of size 2. There are $52\cdot 51 = 2652$ possible hands. But since numbers don’t matter, we can remove the number and just look at suits. There are then 16 possible hands:
\[
\begin{array}{cccc}
♠♠ & ♠♣ & ♠\color{red}{♥} & ♠\color{red}{♦} \\
♣♠ & ♣♣ & ♣\color{red}{♥} & ♣\color{red}{♦} \\
\color{red}{♥}♠ & \color{red}{♥}♣ & \color{red}{♥}\color{red}{♥} & \color{red}{♥}\color{red}{♦} \\
\color{red}{♦}♠ & \color{red}{♦}♣ & \color{red}{♦}\color{red}{♥} & \color{red}{♦}\color{red}{♦}
\end{array}
\]suppose for example that the hands with different suits but same color are not sortable. There are four such hands: $\{ ♠♣, ♣♠, \color{red}{♥}\color{red}{♦}, \color{red}{♦}\color{red}{♥} \}$. Then we might conclude that $\tfrac{4}{16} = 0.25$ of hands are not sortable. But this would be wrong! (thanks to Guy Moore for pointing this out). The reason is that these 16 hands don’t all occur with equal probability. There are $13\cdot 13 = 169$ ways of picking two cards of different suits, but only $13\cdot 12 = 156$ ways of picking two cards of the same suit! So each type of card should be weighted by its likelihood of occurrence. This means the probability of picking an unsortable hand should really be:
\[
\frac{4\cdot 169}{4\cdot 156 + 12\cdot 169} = \frac{13}{51} \approx 0.2549
\]So we can collapse the large number of possible hands to this more manageable size as long as we compensate by appropriately weighting the different items. In combinatorics, this collapsed group of cards (where we care about the order and the suit, but not the number of the card) is an example of a multiset permutation.

Note: If we don’t weigh the hands according to likelihood (i.e. assume all hands are equally likely), this is equivalent to assuming the deck has infinitely many cards (but the same number of cards in each suit).
Adjacent cards of the same suit can be collapsed to a single card. This is because there is never any benefit to moving a block of cards if it splits an existing contiguous block. So for example:
\[
\{♠,♠,\color{red}{♥},\color{red}{♥},\color{red}{♦}\} \implies \{♠,\color{red}{♥},\color{red}{♦}\}.
\]This greatly reduces the computations because if two different hands collapse to the same reduced hand, we only need to do the work for one of them! In computer programming, this technique of storing the values of computations so that you don’t end up computing the same thing many times is called memoization.
If a hand consists of $N$ cards (after we collapse adjacent suit repetitions as described above), then there are $\binom{N+1}{3}$ possible moves we can make. This is because every move is characterized by three locations. I’ll illustrate this with an example. Say our hand looks like this: $\{♠,\color{red}{♥},\color{red}{♦},♠,\color{red}{♥},♣\}$ and we want to sort it by moving cards 4 and 5 and inserting them between cards 1 and 2. In other words, we want to perform:
\[
\{♠,\color{red}{♥},\color{red}{♦},(♠,\color{red}{♥}),♣\} \implies \{♠,(♠,\color{red}{♥}),\color{red}{♥},\color{red}{♦},♣\}
\implies \{♠,\color{red}{♥},\color{red}{♦},♣\}
\]where we did one final collapse at the end. This transformation can be represented by inserting three “bars” between cards as follows. Then the swap simply consists of exchanging the cards between the two pairs of bars! Here is the diagram:
\[
\{\,♠\,\,|\underbrace{\color{red}{♥}\,\,\color{red}{♦}}_\text{swap this}|\underbrace{♠\,\,\color{red}{♥}}_\text{with this}|\,\,♣\,\}
\implies \{♠, ♠, \color{red}{♥},\color{red}{♥},\color{red}{♦},♣\}
\]Each triplet of bars corresponds to a possible swap, and we can insert bars in $N+1$ possible spots (in between each card and on either end).
If the hand has 8 or more cards in it after we perform the collapse, then it’s impossible to sort it in one move. This is because our final hand can have at most 4 groups of cards in it, and a single move can cause at most 3 collapses. With 8 cards, we can get down to 5 but not 4. Here is an example of a 7-card hand that can be sorted in one move:
\[
\{\,♠\,|\,\color{red}{♥}\,♣\,|\,♠\,\color{red}{♥}\,|\,♣\,\color{red}{♦}\,\}
\implies \{\,♠\,♠\,\color{red}{♥}\,\color{red}{♥}\,♣\,♣\,\color{red}{♦}\,\}
\implies \{\,♠\,\color{red}{♥}\,♣\,\color{red}{♦}\,\}
\]

Numerical results

There was a bit of confusion as to how to interpret the meaning of the requirement “the cards of a given suit are grouped together and, if possible, such that no suited groups of the same color are adjacent”. One way to interpret this statement is that we only require cards of a different suit but same color to be separated if it’s actually possible to separate them. For example, the two-card hand $\{\color{red}{♥}\,\color{red}{♦}\}$ counts as being sorted because there are no spades or clubs that can be placed between the heart and the diamond. In this scenario, here are the results:

Note that hands of up to size $N=4$ are always sortable in one move. The exact solution for the case $N=6$ using this interpretation is $\tfrac{51083}{83895} \approx 60.8892\%$. Another way to interpret the requirement is that cards of a different suit but same color can never be adjacent. This means that hands like $\{\color{red}{♥}\,\color{red}{♦}\}$ are never sortable no matter what you do. In this case, we get the result:

It’s interesting to note that the probability of a “sortable” hand in this interpretation actually reaches a maximum when $N=4$ and then decreases afterwards. The exact solution for the case $N=6$ using this interpretation is $\tfrac{1735996}{2936325} \approx 59.1214\%$.

I wrote my code in Julia, and you can view my notebook here. The computation gets slower in an exponential fashion as we increase $N$. The cases $N \leq 8$ take on the order of seconds, but computation time quadruples every time $N$ increases by 1. The last case, $N=13$, took about 40 minutes! While there is a significant amount of computation to do, there is also significant computational overhead due to using exact arithmetic instead of making approximations. For this, I used Julia’s BigInt datatype. So the solutions I obtained were exact. For example, the probability of a sortable hand for $N=13$ (using the first interpretation) turns out to be:
\[
\tfrac{10954472929065768960}{3954242643911239680000} = \tfrac{30785713171}{11112737293000} \approx 0.27703\%
\]Using approximations the whole way leads to a computation time of less than 10 minutes for the most complicated case.

Archipelago

A graph theory problem from the Riddler blog. Here it goes:

You live on the volcanic archipelago of Riddleria. Your archipelago is connected via a network of bridges, forming one unified community. In an effort to conserve resources, the ancient Riddlerians who built this network opted not to build bridges between any two islands that were already connected to the community otherwise. Hence, there is exactly one path from any one island to any other island.

Each island contains exactly one volcano. You know that if a volcano erupts, the subterranean pressure change will be so great that the volcano will collapse in on itself, causing its island — and any connected bridges — to crumble into the ocean. Remarkably, other islands will be spared unless their own volcanoes erupt. But if enough bridges go down, your once-unified archipelagic community could split into several smaller, disjointed communities.

If there were N islands in the archipelago originally and each volcano erupts independently with probability p, how many disjointed communities can you expect to find when you return? What value of p maximizes this number?

Here is my solution:
[Show Solution]

Here is an example of a possible archipelago with $N=9$ islands:

Of course, there are many possible ways to connect the bridges between the islands, but it turns out, this does not actually change the answer to the question! A useful way of thinking about this problem is using the language of graph theory, where we imagine that the archipelago is a graph where the islands are the vertices and the bridges are the edges. Here is the problem, rewritten using graph terminology:

We are given a tree. Each vertex is removed with probability $p$ and when a vertex is removed, so are all incident edges. What is the expected number of connected components in the resulting forest?

We can now leverage a key property of forests:
\[
(\text{# of edges}) + (\text{# of connected components}) = (\text{# of vertices})
\]In particular, if a forest has $N$ vertices and $N-1$ edges, it must have a single connected component, i.e. it’s a tree! You can check this with the tree shown above; it has $9$ vertices and $8$ edges.

This property is very powerful because it tells us exactly what happens when we delete edges and/or vertices from a tree: deleting an edge increases the number of connected components by one, and if we delete an isolated vertex (a vertex with no incident edges), the number of connected components decreases by one.

We now have everything we need to solve the problem. We can write:
\begin{align}
&\mathbb{E}(\text{connected components}) \\
&= \mathbb{E}(\text{vertices})-\mathbb{E}(\text{edges}) \\
&= (N-\mathbb{E}(\text{vertices removed}))-(N-1-\mathbb{E}(\text{edges removed}))\\
&= 1+\mathbb{E}(\text{edges removed})-\mathbb{E}(\text{vertices removed})\\
&= 1+(N-1)\mathrm{Prob}(\text{one edge removed})-N\,\mathrm{Prob}(\text{one vertex removed})\\
&= 1+(N-1)(2p-p^2)-Np \\
&= (1-p)(1-p+Np)
\end{align}In the third-to-last line, we used linearity of expectation. Each vertex has the same probability of being removed ($p$), and each edge also has the same probability of being removed, which is the probability that either of its adjacent vertices of removed, or one minus the probability neither of its vertices are removed: $1-(1-p)^2 = 2p-p^2$.

You might be thinking: but the probability that any two edges disappear is correlated! And it also depends on the connectivity pattern! Correct! but this does not affect the expectation. Linearity of expectation holds for any sum of random variables, even when those variables are correlated.

So here is the final result:

$\displaystyle
\mathbb{E}(\text{# of communities}) = (1-p)(1-p+Np)
$

Maximizing the number of communities

First of all, does this question even make sense? Let’s do a thought experiment. When $p$ is low, there are very few volcanoes, so the archipelago will not be disturbed much (will remain highly connected), so we should obtain a small number of communities. When $p$ is high, most islands disappear, so again we are left with a small number of communities. One might expect that with an intermediate value of $p$, we can maximize the number of communities. The formula we derived before is a concave quadratic function of $p$, and we can maximize it with respect to $p$ by setting the first derivative equal to zero (or completing the square). We obtain:

$\displaystyle
\underset{p}{\mathrm{max}}\,\, \mathbb{E}(\text{# of communities}) = \frac{N^2}{4(N-1)},
\qquad p_{\text{max}} = \frac{N-2}{2(N-1)}
$

As $N$ gets very large, the limiting behavior looks like:

$\displaystyle
\underset{p}{\mathrm{max}}\,\, \mathbb{E}(\text{# of communities}) \approx \frac{N}{4},
\qquad p_{\text{max}} \approx \frac{1}{2}
$

So for a very large archipelago, the probability of volcanic eruption that maximizes the number of communities is $\tfrac{1}{2}$, and the expected number of communities will grow linearly with the number of islands we had initially.

What about distributions rather than expectation?

The natural question to ask after we’ve computed the expected number of communities is to ask about the distribution of the communities. In other words, what are the probabilities associated with obtaining communities of different sizes? It turns out this question doesn’t make any sense unless we know how the islands are connected. For example, consider the simple case of four islands, and two possible ways of connecting them:

In the graph on the right, if we remove vertex $A$, we will have $3$ communities remaining. But in the graph on the left, no matter which combinations of vertices we remove, it’s impossible to obtain $3$ communities. This simple example shows that the distribution of communities must depend on the particular way in which the vertices are connected. Amazingly, though, the expected number of communities is the same no matter how the graph is connected.

Tether your goat!

A geometry problem from the Riddler blog. Here it goes:

A farmer owns a circular field with radius R. If he ties up his goat to the fence that runs along the edge of the field, how long does the goat’s tether need to be so that the goat can graze on exactly half of the field, by area?

Here is my solution:
[Show Solution]

There are many ways to solve this problem, so I figured I would pick an unorthodox way: using calculus! Consider the diagram below.

The circular field is centered at $O$. The goat is tethered at $C$ and the tether has length $r$. Let $B$ be the point diametrically opposed to $C$ and let $A$ be the farthest point along the perimeter that the goat can reach. Let $\theta$ be the angle $BOA$, as shown. If $r$ changes a little bit, by some amount $dr$, then $\theta$ will change by some amount $d\theta$. Likewise, the fraction of the total area reachable by the goat will change by $dA$. Let’s calculate how these quantities are related.

First, note that $dA$ is simply the area between the two circular arcs centered at $C$ and passing through $A$ and $A’$, respectively. We must also divide by the total area of the field in order to get a ratio. As $dr\to 0$, we have:
\[
dA = \frac{\tfrac{1}{2} \theta\, d(r^2)}{\pi R^2}
\]Our next task is to find out how $d\theta$ is related to $d(r^2)$. Look at the triangle $OAC$. By the law of cosines, we have:
$r^2 = 2R^2(1 + \cos(\theta))$. Taking differentials, we obtain:
\[
d(r^2) = -2R^2 \sin(\theta)\,d\theta
\]Substituting back into our expression for $dA$, we obtain:
\[
dA = -\tfrac{1}{\pi} \theta\,\sin(\theta)\,d\theta
\]We can now integrate! Note that $A=1$ when $\theta=0$. Therefore, we have:
\begin{align}
A(\theta) &= 1 + \int_{0}^\theta \tfrac{1}{\pi} t\, \sin(t)\,dt \\
&= \frac{\pi + \theta\,\cos(\theta)-\sin(\theta)}{\pi}
\end{align}This tells us the ratio of areas as a function of $\theta$. If we want this as a function of $r$ instead, we can use the law of cosines again to eliminate $\theta$ and write $A$ as a function of $r$. This yields:
\[
A(\rho) = 1+\frac{\left(\frac{\rho ^2}{2}-1\right) \cos ^{-1}\left(\frac{\rho ^2}{2}-1\right)}{\pi }-\frac{\rho \sqrt{4-\rho ^2}}{2\pi}
\]where I made the substitution $\rho := r/R$. We can plot this quantity as a function of $\rho$, and we get:

As we might expect, when $r=0$, $A=0$. And the area ratio grows monotonically with $r$ until we reach $r=2R$, at which point the goat can reach the farthest point (which is point $B$) and therefore can reach the entire field. There is no closed-form expression for the exact tether length that leads to a particular area ratio, since that would require solving the above $A(\rho)$ for $\rho$. However, we can easily solve the equation numerically. Doing so, we obtain:

$\displaystyle
\frac{\text{length of tether}}{\text{radius of field}} \approx 1.15873
$

Finally, here is a diagram of what the field and tether look like when the area ratio is exactly 1/2.