probability – Book Proofs

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.

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!

Knocking down the last bowling pin

This week’s Riddler classic is a tough one! Here is a paraphrased version of the problem.

Imagine $n^2$ bowling pins arranged in a rhombus. The image to the right illustrates the case $n=3$. We knock down the topmost pin. When any pin gets knocked down, each of the (up to two) pins directly behind it has a probability $p$ of being knocked over (independently of each other). We are interested in the probability that the bottommost pin gets knocked down in the limit of large $n$.

The original problem specifically asked about $p=0.5$ and $p=0.7$.

My solution:
[Show Solution]

Let $F(n,p)$ be the probability that the bottommost pin gets knocked over. The problem asks about $\lim_{n\to\infty} F(n,0.5)$ and $\lim_{n\to\infty} F(n,0.7)$.

Percolation theory

This problem is closely connected to a topic in statistical physics called percolation theory. Statistical physics is concerned with how microscopic properties (e.g., how atoms interact with their neighbors) can lead to macroscopic properties, such as temperature, phase (liquid, solid, gas), and more. We can reframe the problem in statistical physics terms by thinking of the bowling pins as atoms in a porous material. At the microscopic scale, water can squeeze through the space between molecules with probability $p$. Beneath each space, there are two spaces in the next layer. Then, we can ask whether water can percolate through the material or not. In other words, we are asking the macroscopic question: Is the material waterproof?

To get a feel for this, let’s run some simulations. We’ll start with $p=0.5$ and $n=20$. Here, the topmost pin is in the top left, and the bottommost pin is in the bottom right. We can see that the pins never really get rolling because the probability $p$ is too small. In this case, there is no percolation.

Now let’s increase the probability to $p=0.7$. For these plots, I zoomed out and used $n=100$. We can see that in most of the cases, the floodgates open and a large number of pins get knocked down (percolation!). But it doesn’t happen every time. Some cases still run into bad luck and get stopped before the chain reaction can develop.

In percolation theory, the scenario we’re studying is called “directed bond percolation on a 2D square lattice”. It is directed because the pins can only knock other pins over in one direction (downwards), and it is bond percolation because each link (bond) in the lattice is active with probability $p$, as opposed to site percolation where all links are active, but you make a particular pin active with probability $p$.

A fundamental result in percolation theory is that there is a critical probability $p_c$ at which percolation begins to happen. In other words, the probability of percolation $\theta(p)$ satisfies:
\[
\theta(p) \begin{cases}
=0 & \text{for }0 \leq p \lt p_c \\
\gt 0 & \text{for }p_c \lt p \leq 1
\end{cases}
\]The exact value of $p_c$ is only known for certain special cases. For example, it was proved in the following seminal paper that if you have an infinite square lattice and allow pins to knock each other down in all 4 directions, the critical probability is exactly $1/2$:

Harry Kesten. “The Critical Probability of Bond Percolation
on the Square Lattice Equals 1/2”. Commun. Math. Phys. 74, 41-59 (1980). A pdf of this manuscript is freely available here.

The case of directed bond percolation on a square lattice (which is our case of interest) is still an open problem. We know some bounds, however. For example, it was derived in various papers that
\[
0.5176 \leq p_c \leq 0.6298
\]Why is this relevant? Well the problem asks about the cases $p=0.5$ and $p=0.7$, which lie conveniently outside these bounds! Therefore,

When $p=0.5$, the bound above tells us that $p < p_c$. Therefore, $\theta(0.5)=0$ and we have no percolation. So, the probability of knocking over the bottommost pin when $n$ is large is zero.
When $p=0.7$, the bound above tells us that $p > p_c$. Therefore, $\theta(0.7) > 0$ and we have percolation. So, there is a nonzero probability that we will knock over the bottommost pin when $n$ is large.

An analytic lower bound

Given that finding the exact value of $p_c$ is still an open problem, it seems hopeless to look for an analytic solution to the original question involving the bottommost pin. However, we can use an analytic approach to get a bound that allows us to confirm the result for $p=0.5$.

Let’s count the total number of paths that connect the topmost pin to the bottommost pin. Each path consists of $2(n-1)$ steps, and exactly half of the steps are “left” and the other half are “right”. They can be in any order, so the total number of distinct paths is $\binom{2n-2}{n-1}$. Each path also has a probability $p^{2n-2}$ of occurring, since each link has a probability $p$ of being active. Therefore,
\begin{align}
&\mathrm{Prob}(\text{at least one active path})\\
&= 1-\mathrm{Prob}(\text{no active path}) \\
&= 1-\mathrm{Prob}\left(\text{path }k\text{ is inactive for }k=1,\dots,\textstyle\binom{2n-2}{n-1}\right) \\
&\leq 1-\prod_{k=1}^{\binom{2n-2}{n-1}} \mathrm{Prob}(\text{path }k\text{ is inactive}) \\
&=1-\left(1-p^{2n-2}\right)^{\binom{2n-2}{n-1}}
\end{align}Here, the inequality follows from the fact that the paths are not mutually independent (they may have edges in common), and the lower bound comes from assuming independent paths. As an example, suppose there are two paths, and $A$ and $B$ are the events that the paths are inactive, respectively. Then, the probability that $A$ is inactive increases if we know that $B$ is inactive (because the paths may have edges in common). This means that:
\[
P(A\cap B) = P(A\mid B)P(B) \geq P(A)P(B)
\]Note that we have equality if $A$ and $B$ are independent (the paths share no edge in common). Consequently, we now have an analytic upper bound on the probability we care about:

$\displaystyle
F(n,p) \leq 1-\left(1-p^{2n-2}\right)^{\binom{2n-2}{n-1}}
$

Let $L$ be the limit as $n\to\infty$. We can evaluate it as follows:
\begin{align}
L &= 1-\lim_{n\to\infty} \left(1-p^{2n-2}\right)^{\binom{2n-2}{n-1}} \\
&= 1-\lim_{m\to\infty} \left(1-p^{2m}\right)^{\binom{2m}{m}} \\
&= 1-\exp\,\lim_{m\to\infty} \binom{2m}{m} \log \left(1-p^{2m}\right) \\
&= 1-\exp\,\lim_{m\to\infty} \frac{2^{2m}}{\sqrt{m\pi}} \log \left(1-p^{2m}\right) \\
&= 1-\exp\,\lim_{m\to\infty} \frac{4^{m+1} \sqrt{m}\, p^{2 m} \log (p)}{\sqrt{\pi } \left(1-p^{2 m}\right) (m \log (16)-1)} \\
&= 1-\exp\,\lim_{m\to\infty} \frac{4 \cdot (2p)^{2m} \log (p)}{ \log (16)\sqrt{m\pi} } \\
&= 1-\exp\begin{cases}
0 & \text{for }0 \leq p \leq \frac{1}{2}\\
-\infty & \text{for }\frac{1}{2} < p \leq 1 \end{cases}\\ &= \begin{cases} 0 & \text{for }0 \leq p \leq \frac{1}{2}\\ 1 & \text{for }\frac{1}{2} < p \leq 1 \end{cases} \end{align} In the fourth step, we used an asymptotic approximation for the binomial coefficient. In the fifth step, we used L’Hôpital’s rule. Putting everything together, the upper bound on our probability of interest in the limit of large $n$ allows us to conclude that

$\displaystyle
\lim_{n\to\infty} F(n,p) = 0
\quad\text{for } 0 \leq p \leq \frac{1}{2}.
$

So we can answer part of the problem analytically: when $p\leq 0.5$, we can never knock over the bottommost pin in the limit of very large $n$.

Simulation results

To precisely pinpoint the probability of interest in the case $p>0.5$ requires more work, and here we turn to simulation. The first simulation I produced is intended to illustrate the behavior of $F(n,p)$ as $n$ grows.

Here, we can see that for $p \leq 0.5$, the probability of knocking over the bottommost pin clearly tends to zero, as we predicted. Moreover, for $p \geq 0.7$, the probability appears to tend to a constant value, again as predicted. The intermediate cases ($p=0.55, 0.60, 0.65$) are ambiguous, and it isn’t clear whether the limit is zero or if it tends to a positive constant. Note: The slight wiggles in the plot lines are due to the inherent randomness in the simulations. Each point is an average of $50{,}000$ simulations.

The second simulation I ran looked at the limiting probability only, by computing $F(n,p)$ for a large value of $n$. Here, I used $n=300$ and $10{,}000$ simulations per point. Here, we can clearly see the phase transition, which occurs near $p=0.6$. It is difficult to estimate exactly where the transition occurs because this is similar to estimating $p_c$. The closer we get to the critical probability, the longer our typical paths become. So it is not clear which paths are percolating because we’re not using a sufficiently large $n$ to see the true behavior.

Based on the plot, we can see that $\lim_{n\to\infty} F(n,0.7) \approx 0.58$. In order to get more accuracy, I computed how many simulations I would need to run in order to get more accuracy. Using the confidence interval for a Bernoulli process, we can show that an error of $\pm 0.001$ with $95\%$ confidence requires about $1{,}000{,}000$ simulations. This led me to the following slightly more accurate result.

$\displaystyle
\lim_{n\to\infty} F(n,0.7) \approx F(300,0.7) \approx 0.5782 \pm 0.001.
$

Tower of goats

This week’s Riddler classic is a counting problem. Can the goats fit in the tower?

A tower has 10 floors, each of which can accommodate a single goat. Ten goats approach the tower, and each goat has its own (random) preference of floor. Multiple goats can prefer the same floor. One by one, each goat walks up the tower to its preferred room. If the floor is empty, the goat will make itself at home. But if the floor is already occupied by another goat, then it will keep going up until it finds the next empty floor, which it will occupy. But if it does not find any empty floors, the goat will be stuck on the roof of the tower. What is the probability that all 10 goats will have their own floor, meaning no goat is left stranded on the roof of the tower?

My solution:
[Show Solution]

Suppose there are $n$ goats and the building has $n$ floors. Let $a_i$ be the floor preference of goat $i$, with $i=1,\dots,n$. We will say a vector of floor preferences $v=(a_1,a_2,\dots,a_n)$ is “crowded” if it leads to exactly one goat on each floor.

Let $f(n)$ be the number of crowded floor preferences. Since there are $n^n$ total possible preference vectors (each of the $n$ goats can prefer $n$ different floors), the probability that each goat finds a floor is $f(n)/n^n.$ Therefore, solving the problem amounts to finding $f(n)$.

The following elegant counting argument is due to Pollak (1974). We start by considering a modified version of the problem:

Add one extra floor (the roof). The goats are allowed to pick this $(n+1)^\text{st}$ floor as their preferred floor.
Just like the other floors, the roof can accommodate at most one goat.
If a goat ends up on the roof and there is no space for it, the goat “wraps around”; it takes the elevator back down to the first floor and continues its search for an open floor.

In our modified version of the problem, each goat will eventually occupy a floor (possibly the roof), and there will be exactly one empty floor. If the empty floor is the roof, then no goat had the roof as its preferred floor, no goat ever visited the roof, and no goat ever had to take the elevator back down. Therefore, the vector of floor preferences was crowded according to our original definition. Conversely, if one of the other floors ends up empty, then a goat ended up on the roof, which means the original vector of floor preferences was either invalid (one of the goats picked the roof), or the preferences were valid but a goat nonetheless overflowed onto the roof.

So not only is $f(n)$ equal to the number of crowded preference vectors in the original problem, it is also equal to the number of preference vectors in the modified problem that lead to the roof being empty.

The final piece to this argument is to realize that there must be an equal number of preference vectors that lead to the $k^\text{th}$ floor being empty, for $k=1,\dots,n+1$. This follows because of the symmetry in the modified version of the problem. Specifically, if a particular preference vector $(a_1,\dots,a_n)$ leads to floor $k$ being empty, then $(a_1+j,\dots,a_n+j)$ leads to floor $k+j$ being empty (where we wrap all numbers around so they are in the range between $1$ and $n+1$).

Since there are $(n+1)^n$ total possible preference vectors for the modified problem (each of the $n$ goats can prefer $n+1$ different floors), and there are $n+1$ possible empty floors, we obtain
\[
f(n) = \frac{(n+1)^n}{n+1} = (n+1)^{n-1}
\]Converting this number into a probability as explained at the beginning of the solution, we find that the probability that $n$ goats have a crowded set of preferences is

$\displaystyle
p(n) = \frac{(n+1)^{n-1}}{n^n} = \frac{1}{n+1}\left(1+\frac{1}{n}\right)^n
$

For the case $n=10$, we can calcluate:
\[
p(10) = \frac{11^9}{10^{10}} = \frac{2{,}357{,}947{,}691}{10{,}000{,}000{,}000} \approx 23.58\%
\]

Since $\lim_{n\to\infty} \left(1+\frac{1}{n}\right)^n = e$, it follows that in the asymptotic limit of large $n$, we have $p(n) \sim \frac{e}{n+1}$. Here is a plot comparing the true probability to this asymptotic approximation.

About this problem’s history

This sort of counting problem, where you start with a preference and move to the next available open slot, is also known as a parking problem, because its original formulation was about cars trying to park on a one-way street and taking the next available open spot. There are many other equivalent counting problems, which you can find more about in the entry A000272 in the Online Encyclopedia of Integer Sequences. For a wonderful explanation of the parking problem and other related counting problems, take a look at Richard Stanley’s slides on the topic. The solution in my post above is based on the argument in Richard’s slides, which is due to Pollak (1974).