Riddler - Book Proofs

How much can you pull out of a hat?

This week’s Riddler Classic is a strategy game about maximizing payout. What is the optimal strategy?

You start with just the number 1 written on a slip of paper in a hat. Initially, there are no other slips of paper in the hat. You will draw from the hat 100 times, and each time you draw, you have a choice: If the number on the slip of paper you draw is k, then you can either receive k dollars or add k higher numbers to the hat.

For example, if the hat were to contain slips with the numbers 1 through 6 and you drew a 4, you could either receive $4 or receive no money but add four more slips numbered 7, 8, 9 and 10 into the hat. In either case, the slip with the number 4 would then be returned to the hat.

If you play this game perfectly — that is, to maximize the total amount of money you’ll receive after all 100 rounds — how much money would you expect to receive on average?

My solution:
[Show Solution]

When we are playing the game, there are two relevant variables: (1) the number of slips of paper in the hat, and (2) the number of draws remaining. We will therefore define the following quantity:
\[
w(n,t) = \left\{
\begin{array}{l}
\text{Amount we can expect to win from this point} \\
\text{forward, if there are }n\text{ slips of paper in the hat} \\
\text{and we have }t\text{ draws remaining, assuming} \\
\text{optimal play for the rest of the game.}
\end{array}
\right\}
\]The problem is asking us to calculate $w(1,100)$. We will do this by calculating in order: $w(n,0)$, $w(n,1)$, and so on, using a process called dynamic programming. The idea is that the optimal strategy when there are $t+1$ draws remaining depends on the optimal strategy when there are $t$ draws remaining. So if we work our way backwards from $t=0$, we can derive the strategy for any $t$.

When $t=0$, the game is over, so we have $w(n,0)=0$ for all $n$.
When $t=1$, we are on our last turn, so there is no benefit to adding more slips of paper to the hat. In this case, we should always take the money. How much money will this be, on average? There are $n$ slips of paper, so we can expect to win $w(n,1) = \frac{1}{n}\left(1+2+\cdots+n\right) = \frac{n+1}{2}$. Therefore,
\[
w(n,1) = \frac{n+1}{2}.
\]
When $t=2$, we have to think more carefully. Suppose we draw the slip numbered $k$ from the hat. If we take the money, then we will pocket $k$ dollars and then earn on average $w(n,1)$ next turn. If instead we add more slips of paper to the hat, we will pocket nothing, but earn on average $w(n+k,1)$ next turn. Comparing these two quantities:
\[
k + w(n,1) = k + \frac{n+1}{2}\qquad\text{is larger than}\qquad
w(n+k,1) = \frac{n+k+1}{2}
\]Therefore, it is always in our best interest to take the money when there are two turns remaining. When we do this, we will, on average, earn
\[
w(n,2) = \frac{1}{n}\sum_{k=1}^n \left( k + \frac{n+1}{2} \right) = n+1.
\]
The step $t=3$ is similar to $t=2$. We need to decide whether $k+w(n,2)$ is larger or smaller than $w(n+k,2)$. In this case, we have
\[
k+w(n,2) = n+k+1\qquad\text{is equal to}\qquad w(n+k,2) = n+k+1
\]So both choices earn the same on average. In either case, we obtain
\[
w(n,3) = \frac{1}{n}\sum_{k=1}^n \left( n+k+1 \right) = \frac{3}{2}(n+1)
\]
When $t = 4$, it is always best to add more slips of paper. This is because
\[
k+w(n,3) = k+\frac{3}{2}(n+1)\qquad\text{is smaller than}\qquad w(n+k,3)=\frac{3}{2}(n+k+1)
\]And when choosing the larger one, we obtain
\[
w(n,4) = \frac{1}{n}\sum_{k=1}^n \frac{3}{2}(n+k+1) = \frac{9}{4}(n+1).
\]

It turns out we can continue in this fashion and prove that whenever $n\geq 4$, it is always better to add more slips of paper. We can show this using induction. Based on the pattern we’ve observed so far, we observe that $w(n,t) = \left(\frac{3}{2}\right)^{t-2}(n+1)$ for $t=2,3,4$. To show that the pattern continues for larger $t$, suppose it holds for a given $t\geq 3$. To compute $t+1$, we notice that
\[
k+w(n,t) = k+\left(\frac{3}{2}\right)^{t-2}(n+1)\quad\text{is smaller than}\quad w(n+k,t) = \left(\frac{3}{2}\right)^{t-2}(n+k+1)
\]Therefore, we should add more slips of paper, and we obtain
\begin{align}
w(n,t+1) &= \frac{1}{n}\sum_{k=1}^n \left(\frac{3}{2}\right)^{t-2}(n+k+1) \\
&= \frac{1}{n}\left(\frac{3}{2}\right)^{t-2}\sum_{k=1}^n (n+k+1) \\
&= \left(\frac{3}{2}\right)^{t-1}(n+1)
\end{align}which verifies that our pattern continues for $t+1$, and by induction, it continues for all larger $t$. We therefore have the formula:

$\displaystyle
w(n,t) = \begin{cases}
\left(\frac{3}{2}\right)^{t-2}(n+1) & \text{for }t\geq 2 \\
\frac{n+1}{2} & \text{for }t=1
\end{cases}
$

Returning to the original question, we can compute $w(1,100)$:
\[
w(1,100) = \left(\frac{3}{2}\right)^{98}(1+1) = \frac{3^{98}}{2^{97}}
\]So if we play optimally, we can expect to win approximately \$361,387,713,364,635,766.57, or about 361 quadrillion dollars!

Of course, this is only an average. In general, there will be a very large variability in how much we actually win! Here is a plot of the distribution of winnings (one million trials) along with the average shown as a vertical line. The distribution is heavily skewed, so I plotted the winnings on a log scale. This skewness is manifested in the fact that although the average is about $3.6\times 10^{17}$, the median is only $5.3 \times 10^{16}$. So on average we win a lot, but most of the time we win about 7 times less.

N Bottles of Beer

This week’s Riddler Classic is puzzle about the world’s most annoying song.

You and your friends are singing the traditional song, “99 Bottles of Beer.” With each verse, you count down the number of bottles. The first verse contains the lyrics “99 bottles of beer,” the second verse contains the lyrics “98 bottles of beer,” and so on. The last verse contains the lyrics “1 bottle of beer.” There’s just one problem. When completing any given verse, your group of friends has a tendency to forget which verse they’re on. When this happens, you finish the verse you are currently singing and then go back to the beginning of the song (with 99 bottles) on the next verse. For each verse, suppose you have a 1 percent chance of forgetting which verse you are currently singing. On average, how many total verses will you sing in the song?

Extra credit: Instead of “99 Bottles of Beer,” suppose you and your friends are singing “N Bottles of Beer,” where N is some very, very large number. And suppose your collective probability of forgetting where you are in the song is 1/N for each verse. If it takes you an average of K verses to finish the song, what value does the ratio of K/N approach?

My solution:
[Show Solution]

Let’s solve a more general version of the problem. Suppose there are $n$ bottles of beer, and there is a probability $p$ of forgetting at every verse. Define $E_k$ to be the expected number of additional verses we will sing assuming we are currently singing verse $k$. Our task is to solve for $E_1$. Since at each verse, we have a probability $p$ of returning to verse $1$, and a probability $1-p$ of moving to the next verse, we can write:

\begin{align}
E_1 &= 1 + p E_1 + (1-p) E_2 \\
E_2 &= 1 + p E_1 + (1-p) E_3 \\
&\;\vdots \\
E_{n-2} &= 1 + p E_1 + (1-p) E_{n-1} \\
E_{n-1} &= 1 + p E_1 + (1-p)
\end{align}In the last equation, we substituted $E_n=1$, since if we start on the last verse, we will sing that verse and the song will end. We have a system of $n-1$ equations in the unknowns $E_1,E_2,\dots,E_{n-1}$ and our task is to solve for $E_1$.

The obvious approach is to use back-substitution: solve for $E_{n-1}$ in the last equation, substitute into the previous equation to solve for $E_{n-2}$, and continue in this fashion until we reach the first equation. However this leads to a complicated formula. A more efficient approach is to multiply the $k^\text{th}$ equation by $(1-p)^{k-1}$, and obtain:

\begin{align}
E_1 &= 1 + p E_1 + (1-p) E_2 \\
(1-p)E_2 &= (1-p)(1 + p E_1) + (1-p)^2 E_3 \\
(1-p)^2E_3 &= (1-p)^2(1 + p E_1) + (1-p)^3 E_4 \\
&\;\vdots \\
(1-p)^{n-3}E_{n-2} &= (1-p)^{n-3}(1 + p E_1) + (1-p)^{n-2} E_{n-1} \\
(1-p)^{n-2}E_{n-1} &= (1-p)^{n-2}(1 + p E_1) + (1-p)^{n-1} E_n
\end{align}We can get all the terms involving $E_2,\dots,E_{n-1}$ to cancel if we add all these equations together. This yields:

\begin{align}
E_1 &= (1 + p E_1)\sum_{k=0}^{n-2} (1-p)^{k} + (1-p)^{n-1} \\
&= (1+ p E_1) \frac{1-(1-p)^{n-1}}{p} + (1-p)^{n-1} \\
&= \frac{1-(1-p)^{n-1}+p(1-p)^{n-1}}{p} + \bigl(1-(1-p)^{n-1}\bigr)E_1 \\
&= \frac{1-(1-p)^{n}}{p} + \bigl(1-(1-p)^{n-1}\bigr)E_1
\end{align}Now solve for $E_1$ and obtain

$\displaystyle
E_1 = \frac{1-(1-p)^{n}}{p(1-p)^{n-1}}
$

As a simple sanity check, if $n=1$, we obtain $E_1=1$ as we should. It’s also possible to check that if $p\to 0$, we get $E_1\to n$. The problem asks what happens to the ratio $E_1/n$ if we set $p=1/n$ and $n$ is very large. We can write this as

\begin{align}
\lim_{n\to\infty}\left.\frac{E_1}{n}\right|_{p=\frac{1}{n}} &= \lim_{n\to\infty} \frac{1}{n}\frac{1-\left(1-\tfrac{1}{n}\right)^{n}}{\tfrac{1}{n}\left(1-\tfrac{1}{n}\right)^{n-1}} \\
&= \lim_{n\to\infty} \frac{1-\left(1-\tfrac{1}{n}\right)^{n}}{\left(1-\tfrac{1}{n}\right)^{n-1}} \\
&= \frac{1-e^{-1}}{e^{-1}} \\
&= e-1 \\
&\approx 1.71828
\end{align}In the second-last step, we used the well-known fact that $\lim_{n\to\infty} \left(1+\frac{x}{n}\right)^{n} = e^x$, which implies that $\lim_{n\to\infty} \left(1-\frac{1}{n}\right)^{n} = e^{-1}$.

We can verify this limit by plotting $\frac{1}{n}E_1$ with $p=\frac{1}{n}$, which yields

Riddler Football Playoffs

This week’s Riddler Classic is a probability question inspired by the ongoing World Cup.

The Riddler Football Playoff (RFP) consists of four teams. Each team is assigned a random real number between 0 and 1, representing the “quality” of the team. If team $A$ has quality $a$ and team $B$ has quality $b$, then the probability that team $A$ will defeat team $B$ in a game is $\frac{a}{a+b}$.

In the semifinal games of the playoff, the team with the highest quality (the “1 seed”) plays the team with the lowest quality (the “4 seed”), while the other two teams play each other as well. The two teams that win their respective semifinal games then play each other in the final.

On average, what is the quality of the RFP champion?

My solution:
[Show Solution]

Suppose $a \gt b \gt c \gt d$. Consider the event that $A$ is the champion. In order for this to happen, $A$ must first defeat $D$ in the semifinals. Then, either $B$ wins their semifinal match and $A$ defeats $B$ in the finals, or $C$ wins their semifinal match and $A$ defeats $C$ in the finals. We can use a similar argument to compute the probability that $B$, $C$, or $D$ wins, and we obtain the following formulas:
\begin{align}
P_A &= \frac{a}{a+d}\left(\frac{b}{b+c}\cdot\frac{a}{a+b}+\frac{c}{b+c}\cdot\frac{a}{a+c}\right) \\
P_B &= \frac{b}{b+c}\left(\frac{a}{a+d}\cdot\frac{b}{a+b}+\frac{d}{a+d}\cdot\frac{b}{b+d}\right) \\
P_C &= \frac{c}{b+c}\left(\frac{a}{a+d}\cdot\frac{c}{a+c}+\frac{d}{a+d}\cdot\frac{c}{c+d}\right) \\
P_D &= \frac{d}{a+d}\left(\frac{b}{b+c}\cdot\frac{d}{b+d}+\frac{c}{b+c}\cdot\frac{d}{c+d}\right)
\end{align}One can check that $P_A+P_B+P_C+P_D=1$ for all values of $(a,b,c,d)$, since exactly one of these four outcomes must occur. The expected value of the quality of the champion is therefore
\[
x(a,b,c,d) = a P_A + b P_B + c P_C + d P_D.
\]Ultimately, we seek the average of $x(a,b,c,d)$ where each $a,b,c,d$ is chosen uniformly at random in $[0,1]$. The formulas above require $a\gt b \gt c \gt d$, but by symmetry each of the $4!=24$ permutations of $a,b,c,d$ is equally likely to occur. Therefore, we can restrict our averaging to the case where $a\gt b \gt c \gt d$ and multiply the result by $24$. This leads to the following formula for $\bar x$, the average value of $x$.
\[
\bar x = 24 \int_{a=0}^1 \int_{b=0}^a \int_{c=0}^b \int_{d=0}^c x(a,b,c,d)\,
\mathrm{d}d\,\mathrm{d}c\,\mathrm{d}b\,\mathrm{d}a
\]This integral is very difficult to evaluate by hand, so I used Mathematica, and found the result to be:
\begin{align}
\bar x
&= \frac{2}{5} \left(23-(29+2\pi^2)\log(2)+39\log^2(2)-8\log^3(2)-3\zeta(3)\right)\\
&\approx 0.6735417813867959221089
\end{align}Here, $\zeta$ is the Riemann zeta function. The constant $\zeta(3)=\sum_{n=1}^\infty\frac{1}{n^3}\approx 1.2020569$ is known as Apéry’s constant. The reason this constant shows up is because we are repeatedly integrating products of rational functions, and it turns out that
\[
\zeta(3) = \int_0^1\int_0^1\int_0^1\frac{1}{1-xyz}\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z
\]Fun fact: it is known that $\zeta(3)$ is irrational, but it is still not known whether it is transcendental!

Numerical verification

I also obtained an estimate for $\bar x$ by performing a Monte Carlo simulation of 10 million games of Riddler Football Playoff in Matlab. Here is my code:

N = 1e7;

quality = zeros(N,1);
for trial = 1:N
   
    x = sort(rand(4,1));
    a = x(4); b = x(3); c = x(2); d = x(1);
    
    % F = winner of Semifinal 1 (A vs D)
    if rand < a/(a+d)
        f = a;
    else
        f = d;
    end
    
    % G = winner of Semifinal 2 (B vs C)
    if rand < b/(b+c)
        g = b;
    else
        g = c;
    end
    
    % X = winner of Finals (F vs G)
    if rand < f/(f+g)
        quality(trial) = f;
    else
        quality(trial) = g;
    end
end
avg = mean(quality)
sem = std(quality)/sqrt(N);
conf_interval = [avg-1.96*sem,avg+1.96*sem]

Running the code above took about 10 seconds and produced an estimate of $0.67353$, with a 95% confidence interval $[0.67339, 0.67367]$. So it looks like there is good agreement between the numerical value derived from the exact formula and the empirical Monte Carlo estimate.

Randomly cutting a sandwich

This week’s Riddler Classic is geometry puzzle about randomly slicing a square sandwich.

I have made a square sandwich, and now it’s time to slice it. But rather than making a standard horizontal or diagonal cut, I instead pick two random points along the perimeter of the sandwich and make a straight cut from one point to the other. (These points can be on the same side.)

What is the probability that the smaller resulting piece has an area that is at least one-quarter of the whole area?

My solution:
[Show Solution]

Suppose the sandwich has side length 1. Let $t \in [0,1]$ be location of the first point (measured from a corner). We have three cases to consider:

If the second point is on the same side (probability $\frac{1}{4}$), the area of the resulting smallest piece will be zero.
If the second point is on an adjacent side (probability $\frac{1}{2}$), the area of the resulting smallest piece will be $\frac{1}{2}tr$, where $r \in [0,1]$ is the location of the second point on the adjacent side.
If the second point is on the opposite side (probability $\frac{1}{4}$), the area of the resulting piece will be the smaller of the two resulting quadrilaterals. Specifically, it will be $\min(\frac{t+r}{2},1-\frac{t+r}{2})$, where $r \in [0,1]$ is the location of the second point on the opposite side.

We will solve the general case, which is to find the probability that the smallest slice will have a fraction at least $p$ of the total area. the problem asks for $p=\frac{1}{4}$. Since the total area of the square is $1$, we are effectively asking: “what is the probability that the smallest area is at least $p$?”.

Case 1. If we are in the first case (points on the same side), the probability that the area is greater than $p$ is zero, since the area is always zero. So we conclude that:
\[
a_1 = 0
\]

Case 2. If we are in the second case (points on an adjacent sides), we want the probability that $\frac{1}{2}tr > p$, where $r$ and $t$ are chosen uniformly at random in $[0,1]$. This amounts to computing:
\begin{align}
a_2 &= \int_0^1 \int_0^1 \mathbf{1}\bigl( \tfrac{1}{2}tr \gt p \bigr)\,\mathrm{d}r\,\mathrm{d}t \\
&= \int_{2p}^1 \int_0^1 \mathbf{1}\bigl(r \gt \tfrac{2p}{t} \bigr)\,\mathrm{d}r\,\mathrm{d}t \\
&= \int_{2p}^1 \left( 1-\frac{2p}{t}\right)\,\mathrm{d}t \\
&= 1-2p+2p\log(2p)
\end{align}Note that we had to change the lower bound of the integral for $t$ because when $t \lt 2p$, we have $\tfrac{2p}{t}\gt 1$, so the integrand is zero. Geometrically, this is the fraction of the region $(t,r)\in[0,1]^2$ where $tr\gt p^2$, sketched below.

Case 3. If we are in the third case (points on opposite sides), we perform a similar computation with the area of the quadrilateral:
\begin{align}
a_3 &= \int_0^1 \int_0^1 \mathbf{1}\bigl( \min(\tfrac{t+r}{2},1-\tfrac{t+r}{2}) \gt p \bigr)\,\mathrm{d}r\,\mathrm{d}t \\
&= \int_0^1 \int_0^1 \mathbf{1}\bigl( \tfrac{t+r}{2}\gt p\text{ and } 1-\tfrac{t+r}{2} \gt p \bigr) \,\mathrm{d}r\,\mathrm{d}t \\
&= \int_0^1 \int_0^1 \mathbf{1}\bigl( p \lt \tfrac{t+r}{2} \lt 1-p\bigr) \,\mathrm{d}r\,\mathrm{d}t \\
&= \int_0^1 \int_0^1 \mathbf{1}\bigl( 2p \lt t+r \lt 2-2p\bigr) \,\mathrm{d}r\,\mathrm{d}t \\
&= 1-4p^2
\end{align}A trick to evaluating the above integral is to do it geometrically. Sketching the region for which $2p\lt t+r\lt 2-2p$, we see that it is the square of area $1$ with the two corners cut off, whose area is $(2p)^2$, so the result is $1-4p^2$. See below for an illustration.

Putting everything together, the probability that the area of the smallest piece is at least $p$ is given by sum of the probabilities we already calculated, each weighted by their likelihood of occurrence, so $f(p) = \frac{1}{4}a_1 + \frac{1}{2}a_2 + \frac{1}{4}a_3$. Plugging in the results from above and simplifying, we obtain our final answer:

$\displaystyle
f(p) = \frac{3}{4}-p-p^2+p \log (2 p)
$

Here is what the function looks like:

Naturally, the plot only extends to $p=\frac{1}{2}$ because the area of the smaller piece cannot be more than half of the total area. When $p=\frac{1}{4}$ (what the original problem asked about), we obtain $f(\frac{1}{4}) = \frac{7}{16}-\frac{\log (2)}{4}\approx 0.264213$. So the probability that the smaller piece is at least one quarter of the whole area is about 26.4%.

Fall color peak

This week’s Riddler Classic is a seasonal puzzle about leaves changing color.

The trees change color in a rather particular way. Each tree independently begins changing color at a random time between the autumnal equinox and the winter solstice. Then, at a random later time for each tree — between when that tree’s leaves began changing color and the winter solstice — the leaves of that tree will all fall off at once. At a certain time of year, the fraction of trees with changing leaves will peak. What is this maximal fraction?

My solution:
[Show Solution]

Another way to solve the problem, courtesy of Matthew Wallace:
[Show Solution]

Loteria

This week’s Riddler Classic is about Lotería, also known as Mexican bingo!

A thousand people are playing Lotería, also known as Mexican bingo. The game consists of a deck of 54 cards, each with a unique picture. Each player has a board with 16 of the 54 pictures, arranged in a 4-by-4 grid. The boards are randomly generated, such that each board has 16 distinct pictures that are equally likely to be any of the 54.

During the game, one card from the deck is drawn at a time, and anyone whose board includes that card’s picture marks it on their board. A player wins by marking four pictures that form one of four patterns, as exemplified below: any entire row, any entire column, the four corners of the grid and any 2-by-2 square.

After the fourth card has been drawn, there are no winners. What is the probability that there will be exactly one winner when the fifth card is drawn?

My solution:
[Show Solution]

Suppose $m$ players are playing the game ($m=1000$) and there are $n$ distinct cards $(n=54)$. Let’s consider the game from the point of view of one player. Define the following events:
\begin{align}
A &: \left\{\text{No winning pattern in the first four cards}\right\} \\
B &: \left\{\text{No winning pattern in the first five cards}\right\}
\end{align}If we count all the possible winning patterns, we see that there are 18. Meanwhile, there are $\binom{n}{4}$ ways of picking four cards. The probability of no winning pattern is 1 minus the probability of a winning pattern:
\[
\mathbf{P}(A) = 1-\frac{18}{\binom{n}{4}}
= \frac{35137}{35139} \approx 0.999943
\]The approach is similar when there are five cards. Each of the 18 winning patterns uses 4 cards. The fifth card can be any of the remaining $n-4$ cards, so we have:
\[
\mathbf{P}(B) = 1-\frac{18(n-4)}{\binom{n}{5}}
= \frac{35129}{35139} \approx 0.999715
\]

Now let’s ask the question: what is the probability that there is a winning pattern using all five cards given that there is no pattern in the first four cards. We can calculate this using Bayes’ rule and the law of total probability:
\begin{align}
\mathbf{P}(\bar B \mid A)
&= 1-\mathbf{P}(B\mid A) \\
&= 1-\frac{\mathbf{P}(A\cap B)}{\mathbf{P}(A)}\\
&= 1-\frac{\mathbf{P}(B)}{\mathbf{P}(A)}\\
&= \frac{8}{35137} \approx 0.00022768
\end{align}where in the third step, we used the fact that $B \subseteq A$ and therefore $A\cap B = B$ (if there is no winning pattern in the five cards, then there cannot be a winning pattern in the first four).

The problem asks us to find the probability that if $m$ players play the game and nobody wins on the fourth card, exactly one person wins on the fifth card. Let’s call this probability $q$. Since each player is playing independently, there are $m$ different people that can win, with probability $\mathbf{P}(\bar B\mid A)$, and the other $m-1$ players must lose, each with probability $\mathbf{P}(B \mid A)$. So the probability that exactly one person wins on the fifth card is:
\[
q\,=\, m \cdot \mathbf{P}(\bar B\mid A) \cdot \mathbf{P}(B\mid A)^{m-1}
= m \cdot \left( 1-\frac{\mathbf{P}(B)}{\mathbf{P}(A)} \right)\cdot \left( \frac{\mathbf{P}(B)}{\mathbf{P}(A)} \right)^{m-1}
\]Substituting the numerical values found earlier for $\mathbf{P}(A)$ and $\mathbf{P}(B)$ together with $m=1000$, we obtain:
\begin{align}
q &= 1000\left(\frac{8}{35137}\right) \left( \frac{35129}{35137} \right)^{999} \\
&= \frac{8000\cdot 35129^{999}}{35137^{1000}} \approx 0.181356
\end{align}So the probability that exactly one person out of 1000 wins on the fifth card given that nobody won on the first four cards is about 18.136%.

The general case

If we leave the expression for $q$ in terms of $n$ (number of distinct cards) and $m$ (number of players), we obtain:
\[
q = 1728m\cdot \frac{\left(n^4-6 n^3+11 n^2-6 n-2160\right)^{m-1}}{\left(n^4-6 n^3+11 n^2-6 n-432\right)^{m}}
\]If we fix $m=1000$ and vary $n$, we obtain the following graph:

When $n$ is small, it is very likely that many people will win. So the probability that only one person wins on the 5th card is low. But when $n$ is very large, it is very unlikely that anybody will ever win, so the probability that one person wins on the 5th card is low again. The peak occurs at $n=38$ (with probability 36.8%).

If instead we fix $n=54$ and vary $m$, we obtain the following graph:

The shape of the curve is similar. When $m$ is small, there are not many people, so it is not likely that anybody will win on the fifth card. But when $m$ is large, it is likely that many people will win, but unlikely that only one person will win. The peak occurs at $m=4392$ (with probability 36.8% again).

If we vary $m$ and $n$ together, we see that for each $n$, there is an ideal number of players $m$ that would maximize the probability (red dashed line in the plot below).

The red dashed curve has the equation $m = \frac{n^4-6 n^3+11 n^2-6 n-432}{1728}$.

Numerical validation

To test whether the above formulas are correct, I ran a Monte Carlo simulation to estimate $\mathbf{P}(\bar B\mid A)$, which is the probability that a single player wins on the fifth card given that they did not win on the fourth card. I simulated a large number of games and only kept those that didn’t win after four cards. Then I computed the empirical fraction of the remaining cases that won on the fifth card. Here is my Python code.

# winning sets of indices
win_idx = [
    [0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15],  # four in a row
    [0,4,8,12],[1,5,9,13],[2,6,10,14],[3,7,11,15],  # four in a column
    [0,1,4,5],[1,2,5,6],[2,3,6,7],                  # 2x2 in first row
    [4,5,8,9],[5,6,9,10],[6,7,10,11],               # 2x2 in second row
    [8,9,12,13],[9,10,13,14],[10,11,14,15],         # 2x2 in third row
    [0,3,12,15]                                     # four corners
]

# returns True if a bingo card is a winner
def iswinning(matches):
    for idx in win_idx:
        if all([matches[i] for i in idx]):
            return True
    return False

# compute the probability that 5th card wins
# given that the 4th card did not win
N = 10**7

count,total = 0,0

n = 54                   # number of cards
S = list(range(n))       # full set of numbers
C = random.sample(S,16)  # 16 chosen numbers

for iter in range(N):
    F = random.sample(S,5)   # five drawn numbers

    # bingo score card after 4 and 5 turns
    matches4 = [c in F[:-1] for c in C]
    matches5 = [c in F for c in C]
    
    if not iswinning(matches4):
        total += 1
        if iswinning(matches5):
            count += 1

# empirical probability and 95% confidence interval
phat = count/total
conf = 1.96*np.sqrt(phat*(1-phat)/total)

print(f'Analytical: {1728 / (-432 - 6*n + 11*n**2 - 6*n**3 + n**4):.5g}')
print(f'Monte Carlo: {count}/{total} = {phat:.4g}, or [{phat-conf:.4g},{phat+conf:.4g}] (95% CI)')

The result of running this code is:

Analytical: 0.00022768
Monte Carlo: 2292/9999478 = 0.0002292, or [0.0002198,0.0002386] (95% CI)

Unfortunately, since we are trying to estimate such a small probability, a large number of samples are required to reach a narrow confidence interval. I used 10 million samples in the simulation above, which took about 4 minutes on my laptop.

Shared birthdays

This week’s Riddler Classic is a challenging counting problem about shared birthdays.

Suppose people walk into a room, one at a time. Their birthdays happen to be randomly distributed throughout the 365 days of the year (and no one was born on a leap day). The moment two people in the room have the same birthday, no more people enter the room and everyone inside celebrates by eating cake, regardless of whether that common birthday happens to be today.

On average, what is the expected number of people in the room when they eat cake?

Extra credit: Suppose everyone eats cake the moment three people in the room have the same birthday. On average, what is this expected number of people?

My solution:
[Show Solution]

Suppose there are $n$ days in a year, so $n$ different possible birthdays. If there are $k$ people in the room when they eat cake, then this means the first $k-1$ people had different birthdays, and the $k^\text{th}$ person had a birthday that matched one of the first $k-1$ birthdays. The number of ways this can happen is:
\[
\underbrace{n (n-1) (n-2) \cdots (n-k+2)}_{k-1\text{ terms}} \cdot (k-1)
= \frac{n!(k-1)}{(n-k+1)!}
\]Note that we must have $1 \leq k \leq n+1$ because if $k$ is any larger than $n+1$, the first $k-1$ people cannot have distinct birthdays and the product above is zero. Meanwhile, the total number of birthday combinations is $n^k$. So the expected number of people in the room when they eat cake is given by:

$\displaystyle
E_2 = \sum_{k=1}^{n+1} \frac{n!\cdot k(k-1)}{(n-k+1)! \cdot n^k}
$

Unfortunately, this sum does not have a closed-form solution, but we can evaluate it numerically. When $n=365$, we find that the expected number of people in the room is $24.61659$.

This particular sum was studied by Ramanujan, and he found that the sum has an asymptotic expansion (as a function of $n$) that grows like
\[
\sum_{k=1}^{n+1} \frac{n!\cdot k(k-1)}{(n-k+1)! \cdot n^k}
\approx \sqrt{\frac{\pi n}{2}} + \frac{2}{3} + \frac{1}{12}\sqrt{\frac{\pi}{2n}} + \cdots
\]This approximation is quite good. In fact, when $n=365$, we have:
\[
\sqrt{\frac{\pi n}{2}} + \frac{2}{3} = 24.6112
\quad\text{and}\quad
\sqrt{\frac{\pi n}{2}} + \frac{2}{3} + \frac{1}{12}\sqrt{\frac{\pi}{2n}} = 24.6167
\]which are both close to the true value.

Extra credit

When we look for 3 matched birthdays, the counting exercise becomes a lot more complicated. Let’s suppose that there are $\ell$ pairs of matched birthdays before the $k^\text{th}$ person enters the room. To count the number of ways this can happen, we have:

$\binom{n}{\ell}$ ways to pick the $\ell$ days of the year that will have paired birthdays.
$\binom{k-1}{2\ell}$ ways to pick the $2\ell$ people that have a paired birthday.
$\binom{2\ell}{\underbrace{2,\dots,2}_{\ell\text{ times}}} = \frac{(2\ell)!}{2^\ell}$ ways to assign the $\ell$ pairs of birthdays to $2\ell$ people.
$\underbrace{(n-\ell)(n-\ell-1)\cdots(n+\ell-k+2)}_{k-1-2\ell\text{ terms}} = \frac{(n-\ell)!}{(n+\ell-k+1)!}$ ways to assign unique birthdays to the rest of the first $k-1$ people.
$\ell$ ways for the $k^\text{th}$ person to match one of the $\ell$ pairs.

In terms of sum ranges, we must have $1 \leq \ell \leq \frac{k-1}{2}$, since there must be at least one matched pair in the first $k-1$ people, and $2\ell \leq k-1$ because these paired birthdays cannot outnumber the total number of people. Therefore, the expected number of people in the room when there are 3 matched birthdays is:
\[
E_3 = \sum_{k=1}^{n+1}\sum_{\ell=1}^{\left\lfloor\frac{k-1}{2}\right\rfloor}
\binom{n}{\ell}
\binom{k-1}{2\ell}
\frac{k\ell\cdot (n-\ell)!(2\ell)!}{(n+\ell-k+1)!\cdot 2^\ell n^k}
\]With a bit of effort, we can expand the binomials and simplify some of the factorial terms, and we obtain

$\displaystyle
E_3 = \sum_{k=1}^{n+1}\sum_{\ell=1}^{\left\lfloor\frac{k-1}{2}\right\rfloor}
\frac{n!\,k!}{(\ell-1)!\,(k-1-2\ell)!\,(n+\ell-k+1)!\cdot 2^\ell n^k}
$

As you might expect, this expression does not have a closed form either, but again, we can evaluate it numerically, and when $n=365$, we obtain $88.7389$.

Numerical verification

I wrote a simple Python script to simulate the game and see if we can confirm our findings above.

from random import randint
from operator import countOf
import numpy as np

N = 100000 # number of trials
n = 365    # number of days
m = 2      # we celebrate when there are m repeated birthdays

lengths = []
for trials in range(N):
    x = []
    for k in range(1,n+2):
        z = randint(1,n)
        if countOf(x,z) == m-1:
            lengths.append(k)
            break
        else:
            x.append(z)
            
a = np.mean(lengths)                 # empirical mean
sem = np.std(lengths) / np.sqrt(N)   # standard error
print(f'The expected number of people is {a} ± {1.96 * sem:.3f} (95% confidence)')

Running the above with $m=2$, I obtained $24.55923 ± 0.075$.
Running the above with $m=3$, I obtained $88.69134 ± 0.204$. The plus/minus represents a 95% confidence interval. Both simulations agree with the analytical expressions derived above.

Perfect pizza sharing

This week’s Riddler Classic is about how to cut a pizza to achieve precise area ratios between the slices.

Dean made a pizza to share with his three friends. Among the four of them, they each wanted a different amount of pizza. In particular, the ratio of their appetites was 1:2:3:4. Therefore, Dean wants to make two complete, straight cuts (i.e., chords) across the pizza, resulting in four pieces whose areas have a 1:2:3:4 ratio.

Where should Dean make the two slices?

Extra credit: Suppose Dean splits the pizza with more friends. If six people are sharing the pizza and Dean cuts along three chords that intersect at a single point, how close to a 1:2:3:4:5:6 ratio among the areas can he achieve? What if there are eight people sharing the pizza?

My solution:
[Show Solution]

Deriving formulas for the areas

Let’s start by assuming the pizza is a circle of radius 1, and let’s use a coordinate system where the origin is at the center of the pizza. Since all cuts must intersect at a common point on the inside of the circle, we can rotate the circle such that the intersection point lies on the positive $x$-axis. Let the intersection point be $X$ and have coordinates $(x,0)$. We make $n$ cuts through this point, which will result in $n$ chords that intersect the circle at $n$ points above the $x$-axis and $n$ points below, so $2n$ total points.

Let $\theta_i$ denote the angle that point $i$ makes with the positive $x$-axis, measured with respect to the intersection point $X$.
Let $\phi_i$ denote the angle that point $i$ makes with the positive $x$-axis, measured with respect to the origin $O$.
Let $L_i$ denote the length of the line segment between point $i$ and the intersection point $X$.

The figure below illustrates one possible set of cuts for the case $n=3$.

Now we must compute the areas of the different slices. One such slice, between cuts $1$ and $2$, is shaded in the diagram above. We will split each slice into its triangular part $\triangle_{12}$ and its rounded part $\bigcirc_{12}$.

The triangular part is the triangle $XP_1P_2$. Based on the side lengths $XP_1=L_1$ and $XP_2=L_2$ and angle $\angle P_1XP_2=\theta_2-\theta_1$, the area is
\[
\triangle_{12} = \frac{1}{2}L_1L_2 \sin(\theta_2-\theta_1)
\]
The rounded part is the difference between the circular sector $O P_1P_2$ and the triangle $OP_1P_2$. Since the pizza has radius $1$, $OP_1=OP_2=1$. The circular sector has area $\frac{1}{2}(\phi_2-\phi_1)$ and the triangle has area $\frac{1}{2}\sin(\phi_2-\phi_1)$. Therefore, the slice between cuts $1$ and $2$ has total area:
\[
\bigcirc_{12} = \frac{1}{2}\left( (\phi_2-\phi_1)-\sin(\phi_2-\phi_1)\right)
\]

Our next task will be to express $\phi_i$ and $L_i$ in terms of $x$ and $\theta_i$. Consider the point $P_i$. We can write its coordinates in two different ways depending on whether we arrive there via $O\to P_i$ or $O\to X\to P_i$. This yields
\[
P_i = \begin{bmatrix}\cos\phi_i \\ \sin\phi_i\end{bmatrix}
= \begin{bmatrix}x + L_i \cos\theta_i \\ L_i\sin\theta_i \end{bmatrix}
\]squaring these equations and summing them to eliminate $\phi_i$, we obtain
\[
1 = (x+L_i\cos\theta_i)^2+L_i^2\sin^2\theta_i
\]Solving for $L_i$ and keeping the positive root (since $L_i \gt 0$), we find
\[
L_i = \sqrt{1-x^2\sin^2\theta_i}-x\cos\theta_i.
\]Putting everything together, we now have a formula for the area of a slice! In the case of the slice between cuts $1$ and $2$, the area would be:

$\displaystyle
\begin{aligned}
A_{12}(x,\theta_1,\theta_2) &= \frac{1}{2}L_1L_2 \sin(\theta_2-\theta_1) + \frac{1}{2}\bigl( (\phi_2-\phi_1)-\sin(\phi_2-\phi_1)\bigr) \\
\text{where:}\,\,& \left\{\begin{aligned}
L_i &= \sqrt{1-x^2\sin^2\theta_i}-x\cos\theta_i \\
\phi_i&= \arccos\left(\cos\theta_i\sqrt{1-x^2\sin^2\theta_i} + x \sin^2\theta_i \right)
\end{aligned}\right.
\end{aligned}
$

Similar formulas hold true for the remaining slices, so if we have $n$ cuts and we specify the location $x$ of the intersection and the angles $\theta_1,\dots,\theta_n$, we can find formulas for areas of each of the slices $A_{ij}$.

Finding an optimal set of cuts

If we have $n$ cuts, there will be $2n$ slices. To keep things simple, let’s normalize the areas so that the total area is $1+2+\dots+2n=2n^2+n$, and let’s label the normalized areas $A_1,\dots,A_{2n}$ (in counterclockwise order, starting from $A_{12}$). The reason for doing this is that for example in the case $n=2$, we want the $4$ slices to be in the ratio $1:2:3:4$. So we normalize the total area to be $1+2+3+4=10$. That way we can directly try to make $A_1=1$, $A_2=2$, etc.

Each area will be a function of the parameters $\{x,\theta_1,\dots,\theta_n\}$. The goal is to find the set of parameters that give areas that are as close to the numbers $1,2,\dots,2n$. There are many ways to define “closeness”. I will use the least squares criterion, which does a good job in this case because it ensures that all areas are “roughly equally close” to what they should be and there are no large deviations.

For example, in the case $n=2$, our task would be to pick $(x,\theta_1,\theta_2)$ such that the following function is minimized:
\[
J_1 = (A_1-1)^2 + (A_2-2)^2 + (A_3-3)^2 + (A_4-4)^2
\]But there is more… What if the slices do not occur in this particular order 1,2,3,4? We should also check other permutations! For example, we should also try the functions:
\begin{align}
J_2 &= (A_1-1)^2 + (A_2-2)^2 + (A_3-4)^2 + (A_4-3)^2 \\
J_3 &= (A_1-1)^2 + (A_2-3)^2 + (A_3-2)^2 + (A_4-4)^2 \\
&\cdots
\end{align}In all, there will be $(2n)!$ permutations, or $24$ in the case $n=2$. This number can be reduced by a factor of $2$ if we leverage symmetry in reflection.

So our process will be as follows:

pick a permutation $p$ of $\{1,2,\dots,2n\}$
Find $(x,\theta_1,\dots,\theta_n)$ that minimizes the least squares cost
\[
J_p = \sum_{i=1}^{2n} \bigl(A_i(x,\theta_1,\dots,\theta_n)-p_i\bigr)^2
\]subject to the constraints $0\leq x\leq 1$ and $0\leq \theta_1\lt\cdots\lt\theta_n\leq \pi$.
If we can obtain $J_p=0$, then we have a perfect match and we can stop. Otherwise, return to step 1 and pick a different permutation $p$. Repeat until all permutations have been tested.

The function $J(x,\theta_1,\dots,\theta_n)$ is a difficult (read: nonconvex) function, but still relatively well-behaved (read: continuously differentiable), so standard off-the-shelf solvers are adequate for the task. I used Matlab’s constrained nonlinear optimizer fmincon.

To jump straight to the results:
[Show Solution]

Spotting a rare creature

This week’s Riddler Classic is a question about large numbers of attempts at a very unlikely thing.

Graydon is about to depart on a boating expedition that seeks to catch footage of the rare aquatic creature, F. Riddlerius. Every day he is away, he will send a hand-written letter to his new best friend, David Hacker. But if Graydon still has not spotted the creature after $n$ days (where $n$ is some very, very large number), he will return home.

Knowing the value of $n$, Graydon confides to David there is only a 50 percent chance of the expedition ending in success before the $n$ days have passed. But as soon as any footage is collected, he will immediately return home (after sending a letter that day, of course).

On average, for what fraction of the $n$ days should David expect to receive a letter?

My solution:
[Show Solution]

We will assume the probability of spotting the creature on any given day is $p$ and is independent one day to the next. We are told the probability the expedition will fail is $\frac{1}{2}$. This occurs when the creature is not observed for $n$ days. This event has probability $(1-p)^n$. Therefore, we have
\[
(1-p)^n = \frac{1}{2}
\]We can solve this equation for $p$ and we find
\[
p = 1-2^{-1/n}
\]

If Graydon sends $k$ letters (with $1\leq k \leq n-1$), this means that the creature was spotted on the $k^\text{th}$ day, after not being spotted for the first $k-1$ days. The probability of this event is $p(1-p)^{k-1}$. If $n$ letters are sent, then either the creature was spotted on the last day (probability $p(1-p)^{n-1}$), or the creature was not spotted at all (probability $(1-p)^n$). Let the average fraction of days where a letter is received be $A$. This quantity is equal to the expected number of letters received divided by $n$. Therefore, we obtain:
\[
A = \frac{1}{n} \sum_{k=1}^n k p(1-p)^{k-1} + (1-p)^n
\]We can evaluate this sum in closed form, and we obtain
\[
A = \frac{1-(1-p)^n}{n p}
\]Substituting what we found for $p$, we obtain
\[
A = \frac{1}{2n(1-2^{-1/n})}
\]We are told that “$n$ is very very large”, so it makes sense to examine what happens as $n\to\infty$. Doing so, we observe that $A$ tends to a constant:

With a bit of work, we can actually evaluate the limit exactly:
\[
\lim_{n\to\infty} \frac{1}{2n(1-2^{-1/n})}
= \frac{1}{\log (4)} \approx 0.721348
\]So, on average, David should expect to receive a letter on approximately 72.1% of the days.

Alternative approach

A slightly faster way to get to the same answer is to realize that if $(1-p)^n=\frac{1}{2}$ and $n$ is “very very large”, then it must be true that $p$ is “very very small”. Specifically, if we recall the fact that
\[
\lim_{n\to\infty} \left(1+\frac{x}{n}\right)^n = e^x
\]Then we can rearrange:
\[
(1-p)^n = \left(1 + \frac{(-np)}{n}\right)^n
\]So if $n\to\infty$, we conclude that $e^{-np} \to \frac{1}{2}$, or that $np \to \log(2)$. Therefore, we can return to our original expression for $A$ and substitute $(1-p)^n \to \tfrac{1}{2}$ and $np\to \log(2)$ and we find:
\[
A = \frac{1-(1-p)^n}{n p} \to \frac{1}{\log(4)}
\]which is the same thing we found using the first approach.

Pill splitting

his week’s Riddler classic is about splitting pills to get the right dose.

I’ve been prescribed to take 1.5 pills of a certain medication every day for 10 days, so I have a bottle with 15 pills. Each morning, I take two pills out of the bottle at random.

On the first morning, these are guaranteed to be two full pills. I consume one of them, split the other in half using a precision blade, consume half of that second pill, and place the remaining half back into the bottle.

On subsequent mornings when I take out two pills, there are three possibilities:

I get two full pills. As on the first morning, I split one and place the unused half back into the bottle.
I get one full pill and one half-pill, both of which I consume.
I get two half-pills. In this case, I take out another pill at random. If it’s a half-pill, then I consume all three halves. But if it’s a full pill, I split it and place the unused half back in the bottle.

Assume that each pill — whether it is a full pill or a half-pill — is equally likely to be taken out of the bottle.

On the 10th day, I again take out two pills and consume them. In a rush, I immediately throw the bottle in the trash before bothering to check whether I had just consumed full pills or half-pills. What’s the probability that I took the full dosage, meaning I don’t have to dig through the trash for a remaining half-pill?

My solution:
[Show Solution]

Problems of this sort may seem tedious at first glance, but there is a simple and systematic way to approach them. We can encode the pills in the bottle using polynomials. Specifically, we’ll write $x^n y^m$ when there are $n$ full pills and $m$ half-pills. Since there is an element of randomness with how the pills are picked, we will use a weighted sum of such polynomials to describe a mixture or possible states. For example, we could write:
\[
\tfrac{2}{3} x y + \tfrac{1}{3} y^3
\]to describe a scenario where there is a $\frac{2}{3}$ probability that the bottle contains one full pill and one half-pill, and a $\frac{1}{3}$ probability that the bottle contains three half-pills. Note that the coefficients sum to $1$, so exactly one of the two scenarios must occur.

So what happens when we take our daily dose? We can consider each scenario in turn.

With probability $\frac{\binom{n}{2}}{\binom{n+m}{2}}$, we will pick two full pills. In this case, we cut one in half and return half to the bottle. This means the bottle now contains $x^{n-2}y^{m+1}$.
With probability $\frac{nm}{\binom{n+m}{2}}$, we will pick one full pill and one half-pill. In this case, we have our exact dose and the bottle now contains $x^{n-1}y^{m-1}$.
With probability $\frac{\binom{m}{2}}{\binom{n+m}{2}}$, we will pick two half-pills. In this case, we must go back into the bottle and pick one more pill. There are two sub-cases here:
- With probability $\frac{n}{n+m-2}$, we pick one of the full pills. So we split it in half, and return half to the bottle. The jar has now lost a total of one full pill and one half-pill, so it contains $x^{n-1}y^{m-1}$.
- With probability $\frac{m-2}{n+m-2}$, we pick one of the half-pills. So we have our exact dose and the bottle now contains $x^n y^{m-3}$.

Combining all the possible cases, we find that every day, the polynomial describing the contents of the bottle change according to the mapping
\begin{multline}
x^n y^m \mapsto \tfrac{n (n-1)}{(m+n) (m+n-1)}x^{n-2} y^{m+1} + \tfrac{m n (3 m+2 n-5)}{(m+n) (m+n-1) (m+n-2)} x^{n-1} y^{m-1} \\
+ \tfrac{m (m-1) (m-2)}{(m+n) (m+n-1) (m+n-2)} x^n y^{m-3}
\end{multline}
There are two nice things about this formulation. First, edge cases are taken care of automatically. For example, if the bottle contains $n=2$ full pills and $m=2$ half-pills, the third term ($x^n y^{m-3}$) in the expansion is impossible, since the bottle does not contain three half-pills. But the associated coefficient is zero, so this term cancels automatically. Similarly, if $m=0$ (the bottle only contains full pills), two of the coefficients vanish and the third becomes $1$, so we obtain $x^n \mapsto x^{n-2} y$.

Second, when the bottle contains a mixture of possible states (described by a polynomial whose coefficients sum to 1), we can apply the above transformation to each term of the polynomial separately, and the final result will again be a mixture of possible states. This follows because the sum of the coefficients of the mapping above sum to 1 as well.

To compute the mixture of states after 10 days, I wrote a Mathematica script that starts at $x^{15}$ (15 full pills) and applies the above transformation recursively. Here is the Mathematica script:

P = x^15; (* initial polynomial *)
Print[P];
For[iter = 0, iter < 9, iter++,
  Q = 0;
  ML = MonomialList[P];
  For[i = 1, i <= Length[ML], i++,
    term = ML[[i]];
    coef = term /. {x -> 1, y -> 1};
    n = Exponent[ term, x];
    m = Exponent[term, y];
    Q = Q + coef ((m(m-1)(m-2))/((n+m)(n+m-1)(n+m-2)) x^n y^(m-3)
          + (n m(2n+3m-5))/((n+m)(n+m-1)(n+m-2)) x^(n-1) y^(m-1)
          + (n(n-1))/((n+m)(n+m-1)) x^(n-2) y^(m+1));
   ];
   P = Q // Expand;
   Print[P];
]

and here is the output from the script:
\begin{gather}
x^{15} \\
x^{13} y\\
\tfrac{1}{7}x^{12}+\tfrac{6 }{7}x^{11} y^2\\
\tfrac{36 }{91}x^{10} y+\tfrac{55 }{91}x^9 y^3\\
\tfrac{23 }{308}x^9+\tfrac{2385}{4004} x^8 y^2+\tfrac{30}{91} x^7 y^4\\
\tfrac{4 }{13}x^7 y+\tfrac{81}{143} x^6 y^3+\tfrac{18}{143} x^5 y^5\\
\tfrac{335}{4004} x^6+\tfrac{87}{154} x^5 y^2+\tfrac{185}{572} x^4 y^4+\tfrac{4}{143} x^3 y^6\\
\tfrac{22573}{56056} x^4 y+\tfrac{42605}{84084} x^3 y^3+\tfrac{101}{1144} x^2 y^5+\tfrac{1}{429}x y^7\\
\tfrac{44783}{240240} x^3+\tfrac{362603}{560560} x^2 y^2+\tfrac{27185}{168168} x y^4 + \tfrac{61}{12012} y^6\\
\tfrac{80529}{101920} x y+\tfrac{21391}{101920} y^3\\
\end{gather}

Therefore, on the $10^\text{th}$ day, there is a $\tfrac{21391}{101920} \approx 20.988\%$ chance the bottle contains three half-pills, and a $\tfrac{80529}{101920} \approx 79.012\%$ chance the bottle contains one full pill and one half-pill. This is the probability that the correct dosage was taken, which is what is asked in the original problem statement.

Large-pill limit

We saw that on the final day, the probability that the bottle contains one full pill and one half-pill is about $79\%$. A natural question to ask is: what happens to this probability as we vary the number of starting pills? If we start with only 3 pills, then we always have one full pill and one half-pill on the next day. So the probability is $100\%$. Here is a plot of what happens as we vary the number of initial pills:

It appears that as we increase the number of initial pills, the probability of having one full pill and one half-pill on the last day decreases. I don’t have a great explanation for why this happens; if you made it this far and have any ideas, write a comment; I would love to hear it!

Generating functions?

I am well aware that the polynomials I used look a lot like generating functions. I looked for a nice way to express the pill-taking transformation as a differential operator acting on the generating function, but I was thwarted by the denominators. The numerators are easy, e.g.
\[
\tfrac{\mathrm{d}^3}{\mathrm{d}y^3}(x^n y^m) = m(m-1)(m-2) x^n y^{m-3}
\]but I couldn’t find a nice way to deal with the denominators. For the case above, the denominator should be $(n+m)(n+m-1)(n+m-2)$. If anybody found a way to solve this using the full generating function machinery, I would love to hear about it!

Tag: Riddler

How much can you pull out of a hat?

N Bottles of Beer

Riddler Football Playoffs

Numerical verification

Randomly cutting a sandwich

Fall color peak

Loteria

The general case

Numerical validation

Shared birthdays

Extra credit

Numerical verification

Perfect pizza sharing

Deriving formulas for the areas

Finding an optimal set of cuts

The case of 4 slices

The case of 6 slices

The case of 8 slices

Spotting a rare creature

Alternative approach

Pill splitting

Large-pill limit

Generating functions?