counting - Book Proofs

Halloween Puzzle

This week’s Fiddler is about rounding!

You are presented with a bag of treats, which contains $n \geq 3$ peanut butter cups and some unknown quantity of candy corn kernels (with any amount being equally likely). You reach into the bag $k$ times, with $3 \leq k \leq n$, and pull out a candy at random. Each time, it’s a peanut butter cup! How many candy kernels do you expect to be in the bag?

My solution:
[Show Solution]

Define the following random variables:

$m$: number of candy corn kernels in the bag
$K$: whether we draw $k$ peanut butter cups in a row (true/false)

We are asked to calculate the expected number of corn kernels in the bag given that we draw $k$ peanut butter cups in a row. In other words, we want to find $\mathbb{E}(m\mid K)$.

An easier quantity to calculate is the probability that we draw $k$ peanut butter cups in a row given that there are $m$ candy corn kernels in the bag. This is given by:
\begin{align*}
\mathbb{P}(K\mid m)
&= \frac{n}{n+m}\cdot \frac{n-1}{n+m-1}\cdots\frac{n-k+1}{n+m-k+1} \\
&= \frac{n!(n+m-k)!}{(n-k)!(n+m)!} \\
&= \binom{n}{k}\binom{n+m}{k}^{-1}
\end{align*}Using Bayes’ rule, we can express the desired expectation in terms of the conditional probability above:
\begin{align*}
\mathbb{E}(m\mid K)
&= \sum_{m=0}^\infty m\cdot \mathbb{P}(m\mid K)
= \sum_{m=0}^\infty m\cdot \frac{\mathbb{P}(K \mid m)\mathbb{P}(m)}{\mathbb{P}(K)}
\end{align*} Now substitute: $\mathbb{P}(K) = \sum_{m=0}^\infty \mathbb{P}(K\mid m)\mathbb{P}(m)$ and obtain:
\begin{align*}
\mathbb{E}(m\mid K)
&= \frac{\sum_{m=0}^\infty m\cdot \mathbb{P}(K \mid m)\mathbb{P}(m)}{\sum_{m=0}^\infty \mathbb{P}(K \mid m)\mathbb{P}(m)} \\
&= \frac{\sum_{m=0}^\infty m\cdot \mathbb{P}(K \mid m)}{\sum_{m=0}^\infty \mathbb{P}(K \mid m)}
\end{align*}where in the last step, we used the fact that $\mathbb{P}(m)$ is the same for all $m$ so we canceled it from the numerator and denominator. Substituting our expression for $\mathbb{P}(K\mid m)$, we obtain:
\begin{align*}
\mathbb{E}(m\mid K)
&= \frac{\sum_{m=0}^\infty m\binom{n}{k}\binom{n+m}{k}^{-1}}{\sum_{m=0}^\infty \binom{n}{k}\binom{n+m}{k}^{-1}}
= \frac{\sum_{m=0}^\infty m\binom{n+m}{k}^{-1}}{\sum_{m=0}^\infty \binom{n+m}{k}^{-1}}
%= \frac{n-k+1}{k-2}
\end{align*}
To simplify this, define the following quantity:
\[
f(k,a,b) := \sum_{i=a}^b \binom{i}{k}^{-1}
\]and observe that we have:
\[
\mathbb{E}(m\mid K) = \frac{\sum_{j=1}^\infty f(k,n+j,\infty)}{f(k,n,\infty)}
\]It turns out we can evaluate $f$, and that it has the nice closed-form expression:
\[
f(k,a,b) = \frac{k}{k-1}\left( \binom{a-1}{k-1}^{-1}-\binom{b}{k-1}^{-1}\right)
\]This formula can be proved by induction and is related to something called the German tank problem.

With this formula in hand, we can take the limit $b\to\infty$ and see that:
\[
f(k,a,\infty) = \frac{k}{k-1}\binom{a-1}{k-1}^{-1}
\]Substituting into our expression for the desired expectation, we obtain:
\begin{align*}
\mathbb{E}(m\mid K)
&= \frac{\sum_{j=1}^\infty f(k,n+i,\infty)}{f(k,n,\infty)} \\
&= \frac{\sum_{j=1}^\infty \frac{k}{k-1}\binom{n+j-1}{k-1}^{-1}}{\frac{k}{k-1}\binom{n-1}{k-1}^{-1}} \\
&= \binom{n-1}{k-1}\sum_{j=1}^\infty \binom{n+j-1}{k-1}^{-1} \\
&= \binom{n-1}{k-1} f(k-1,n,\infty) \\
&= \binom{n-1}{k-1} \frac{k-1}{k-2} \binom{n-1}{k-2}^{-1} \\
&= \frac{(n-1)!}{(n-k)!(k-1)!} \frac{k-1}{k-2} \frac{(n-k+1)!(k-2)!}{(n-1)!} \\
&= \frac{n-k+1}{k-2}
\end{align*}

Therefore, our final answer is:

$\displaystyle
\mathbb{E}(m\mid K) = \frac{n-k+1}{k-2}
$

I undoubtedly found the most complicated possible solution for this problem… So if somebody can show me a more elegant solution (perhaps a counting argument?) then I would be much obliged!

Tiling a Tilted Square

This week’s Fiddler is a challenging counting problem.

Consider the following array of 25 squares:

You are filling the array with rectangles by repeating the following two steps:

Select one of the 12 squares along the outer perimeter that has not yet been selected as part of a rectangle.
Form the largest rectangle you can that includes the square you just selected and other squares that are not yet part of any such rectangle.

You repeat these steps until every square along the perimeter has been selected. Here are two final states you might encounter:

How many distinct final states are possible? (Note: States that are rotations or reflections of each other should be counted as distinct.)

My solution:
[Show Solution]

We will solve the problem in three steps, where each time we build up from the previous step and add layers of complexity.

Step 1: one quarter

First, consider the original problem, but with one quarter of the shape:

If we ask about tiling this shape using the procedure outlined in the problem statement, then we can write a recursion in the general case. To see why, start by picking one of the edge squares. This divides the shape into two similar but smaller shapes. For example, if there are $n=6$ edge squares and we pick the third one from the top, we obtain:

We are left with the cases $n=2$ and $n=3$, which we can further subdivide. If we call $c_n$ the number of tilings of this shape, we therefore have the recursion:
\[
c_0=1,\quad\text{and}\quad c_{n+1} = \sum_{k=0}^n c_k c_{n-k}\quad\text{for }n\geq 0
\]This is the well-known recurrence relation for the Catalan numbers. The first few Catalan numbers are (starting from $n=0$):
\[
\{c_n\} = \{1, 1, 2, 5, 14, 42, 132, 429, 1430,\dots\}
\]and a general formula is given by:
\[
c_n = \frac{1}{n+1}\binom{2n}{n}
\]

Step 2: one half

Now, consider the original problem, but with one half of the shape:

If we ask about tiling this shape, we can form a recurrence relation like we did in Step 1. This time, when we pick an edge square, we draw one large rectangle and we are left with a half-shape (up top) and identical quarter-shapes (on either side). For example, if there are $n=6$ edge squares and we pick the fourth from the top, we obtain:

If we call $d_n$ the number of tilings of this shape, we therefore have the recursion:
\[
d_0=1,\quad\text{and}\quad
d_{n+1} = \sum_{k=0}^n c_{n-k}^2 d_k\quad\text{for }n\geq 0
\]The first few numbers in this sequence are (starting from $n=0$):
\[
\{d_n\} = \{1, 1, 2, 7, 38, 274, 2350, 22531, 233292,\dots\}
\]

Step 3: the whole thing

Now consider the original problem. You can probably guess the pattern by now… By picking one of the edge squares, we subdivide the problem into four half-shapes in identical pairs (north-south and east-west). If we call $e_n$ the number of tilings of the whole shape, we therefore have the formula:
\[
e_0=1,\qquad\text{and}\quad
e_{n+1} = \sum_{k=0}^n d_{n-k}^2 d_k^2\quad\text{for }n\geq 0
\]The first few numbers in this sequence are (starting from $n=0$):
\[
\{e_n\} = \{1, 1, 2, 9, 106, 3002, 153432, 11209105, 1027079042\}
\]

Unfortunately, none of this is particularly satisfying since we do not have a closed-form solution for $e_n$. Let’s try to find one…

Attempt at a closed-form solution

A good place to start is to see how the closed-form solution for the Catalan numbers is derived. One way is to use generating functions. The idea is that we define an infinite polynomial (a power series) where the coefficients are the sequence we care about. For Catalan numbers, we have:
\[
C(x) = \sum_{n=0}^\infty c_n x^n
\]Now notice that the recurrence relation for Catalan numbers is a convolution, which we can obtain by squaring $C(x)$:
\begin{align*}
C(x)^2 &= \sum_{m=0}^\infty \sum_{k=0}^\infty c_m c_k x^{m+k} \\
&= \sum_{n=0}^\infty \left(\sum_{k=0}^n c_k c_{n-k} \right)x^n \\
&= \sum_{n=0}^\infty c_{n+1} x^n \\
&= \frac{1}{x}\left( C(x)-1\right)
\end{align*}Solving for $C(x)$, we obtain:
\[
C(x) = \frac{1-\sqrt{1-4x}}{2x}
\](The other root of the quadratic can be excluded since it does not satisfy $C(0)=c_0=1$) From here, we can perform a series expansion via the binomial theorem and extract the coefficient of $x^n$, which yields the formula $c_n = \frac{1}{n+1}\binom{2n}{n}$.

We can use a similar argument to obtain a generating function for $d_n$. To this effect, define:
\[
D(x) = \sum_{n=0}^\infty d_n x^n
\]Now, we can write:
\begin{align*}
D(x) &= 1 + \sum_{n=0}^\infty d_{n+1} x^{n+1} \\
&= 1 + \sum_{n=0}^\infty \sum_{k=0}^n c_{n-k}^2 d_k x^{n+1} \\
&= 1 + \sum_{k=0}^\infty \sum_{n=k}^\infty c_{n-k}^2 d_k x^{n+1} \\
&= 1 + \sum_{k=0}^\infty \sum_{n=0}^\infty c_{n}^2 d_k x^{n+k+1} \\
&= 1 + \sum_{k=0}^\infty d_k x^k \sum_{n=0}^\infty c_n^2 x^{n+1} \\
&= 1 + D(x)\sum_{n=0}^\infty c_n^2 x^{n+1}
\end{align*}Therefore, we can express the generating function for $d_n$ in terms of the generating function of squared Catalan numbers:
\[
\hat C(x) := \sum_{n=0}^\infty c_n^2 x^{n+1}\qquad\text{and}\qquad
D(x) = \frac{1}{1-\hat C(x)}
\]According to Mathematica, the series involving squared Catalan numbers can be evaluated in terms of a hypergeometric function. Namely:
\[
\hat C(x) =
\frac{1}{4} \bigl(\, _2F_1(-\tfrac{1}{2},-\tfrac{1}{2};1;16 x)-1\bigr)
\]Unfortunately, this isn’t particularly helpful as there does not appear to be any way to obtain a formula for the coefficient of $x^n$ in the series for $D(x)$.

Continuing in this fashion, we can also define $E(x)$ as the generating function for $e_n$, which we can express in terms of the square of the generating function for $d_n^2$. Namely:
\[
E(x) = \sum_{n=0}^\infty e_n x^n,\qquad
\hat D(x) = \sum_{n=0}^\infty d_n^2 x^{n+1},\qquad
E(x) = 1 + \tfrac{1}{x}\hat D(x)^2
\]But again, this doesn’t really seem helpful as we can’t evaluate $D(x)$ or much less $\hat D(x)$.

If anybody else can make progress on this problem I would love to hear your approach!

Asymptotics

It was pointed out by commenter MarkS that $c_n^2/d_n$ and $c_n^4/e_n$ appear to tend to finite limits as $n\to\infty$. Here is what we get when we plot these quantities up to $n=2000$:

Is there a way we can find the exact values of these limits? Maybe! We’ll make use of the following fact. If the sequences $a_n$ and $b_n$ have corresponding generating functions $A(x)$ and $B(x)$, and these have radius of convergence $R$, then we can write:
\[
\lim_{n\to\infty}\frac{a_n}{b_n} = \lim_{x\to R^-} \frac{A(x)}{B(x)}
\]This works so long as $A(x)$ and $B(x)$ go to $\infty$ as $x\to R^-$, because any finite truncation of the series will be dominated by its last term (largest power of $x$), so the ratio of truncated series just behaves like the ratio of its last terms, which is what we care about (apologies for the hand-waving; hopefully this makes sense!).

Applying this idea, we would like to write:
\[
\lim_{n\to\infty}\frac{c_n^2}{d_n} = \lim_{x\to R^-}\frac{\hat C(x)}{x D(x)} = \lim_{x\to R^-} x \hat C(x)\bigl(1-\hat C(x)\bigr)
\]But there is just one problem… We already established that
\[
\hat C(x) = \sum_{n=0}^\infty c_n^2 x^{n+1} =
\frac{1}{4} \bigl(\, _2F_1(-\tfrac{1}{2},-\tfrac{1}{2};1;16 x)-1\bigr)
\]and as it turns out, the series for this function has a radius of convergence of $R=\frac{1}{16}$. If you’re curious as to why, remember that this is a series where the coefficients are squared Catalan numbers. Catalan numbers have a well-known property that
\[
c_n \sim \frac{4^n}{n^{3/2}\pi},\qquad\text{i.e.,}\quad
\lim_{n\to\infty}\frac{c_n}{\left(\frac{4^n}{n^{3/2}\pi}\right)} = 1
\]Therefore, $c_n^2 \sim 16^n/n^3\pi$, so a necessary condition for $\hat C(x)$ to converge is that $|x|\leq \frac{1}{16}$ (the general term of the series must tend to zero).

So what’s the problem? Here is a 3D plot of $\hat C(z)$ for complex $z$ (vertical axis is magnitude, color-coded by argument) and plotted for $|z|\leq \frac{1}{16}$.

As we can see, nothing is going to infinity near the boundary. How is it possible? It turns out the radius of convergence is $\frac{1}{16}$ because there is a different kind of discontinuity at this point; the argument (rather than the magnitude) becomes discontinuous. You can see it more clearly when the plot is zoomed out to $|z|\leq 1$ (notice the color discontinuity along the positive real axis).

So how do we deal with this problem? We need to transform the series so that it goes to infinity as $x\to\frac{1}{16}$. One way to do this is to take derivatives. This will yield a series with terms like $n a_n x^{n-1}$ and $n b_n x^{n-1}$, and the ratio remains unchanged! Here is a plot of the first derivative, $\hat C'(z)$:

This is better, as now the function appears non-differentiable near $z=\frac{1}{16}$, but it still doesn’t go to infinity. Let’s differentiate again! Here is a plot of $\hat C'{}'(z)$:

Now we’re in business. So the limit we’re looking for is:
\[
\lim_{n\to\infty}\frac{c_n^2}{d_n}
= \lim_{x\to\frac{1}{16}} \frac{\frac{d^2}{dx^2}\Bigl(\frac{1}{x}\hat C(x)\Bigr)}{\frac{d^2}{dx^2}\Bigl(\frac{1}{1-\hat C(x)}\Bigr)}=\left(5-\frac{4}{\pi}\right)^2\approx 13.88874349
\]I used Mathematica in the last step, which can analytically differentiate and evaluate hypergeometric functions. In conclusion, we have the asymptotic formula:

$\displaystyle
d_n \sim \frac{c_n^2}{\left( 5-\tfrac{4}{\pi}\right)^2}
\sim \frac{16^n}{\pi\left( 5-\tfrac{4}{\pi}\right)^2 n^3}
$

Now as for $e_n$, I’m stuck again. In principle we could use the same technique, which would require evaluating:
\[
\lim_{n\to\infty}\frac{c_n^4}{e_n}
= \lim_{x\to R^-} \frac{\frac{d^k}{dx^k}\sum_{n=0}^\infty c_n^4 x^n}{\frac{d^k}{dx^k}E(x)}
= \lim_{x\to R^-} \frac{\frac{d^k}{dx^k}\sum_{n=0}^\infty c_n^4 x^n}{\frac{d^k}{dx^k}\left(1+\frac{1}{x}\hat D(x)^2\right)}
\]where we differentiate however many times necessary to ensure the series diverges as $x\to R^-$. Mathematica is able to evaluate the numerator (it’s hypergeometric functions again, but a different kind this time), so I was able to determine that $R = \frac{1}{256}$ and the smallest $k$ we can use is $5$ (yikes!). But the sticking point is the denominator. Although we found an asymptotic expansion for $d_n$, we don’t have an actual formula. Without this, we can’t evaluate the general term of $\hat D(x)^2$, which depends on all the $d_n$.

Again, if anybody has ideas, let me know in the comments!

Picking a speaker at random

This week’s Fiddler is about selecting a speaker of the house at random. How long will it take?

There are three candidates who want the job of Speaker. All 221 members of the party vote by picking randomly from among the candidates. If one candidate earns the majority of the votes, they become the next Speaker. Otherwise, the candidate with the fewest votes is eliminated and the process is repeated with one less candidate. If two or more candidates receive the same smallest number of votes, then exactly one of them is eliminated at random. What is the average number of rounds needed to select a new Speaker?

Extra Credit
What if there were 10 candidates running for Speaker?

My solution:
[Show Solution]

We will consider a more general case; suppose there are $n$ candidates and $2m+1$ voting members in the party. We write $2m+1$ since there should be an odd number of party members to break ties. Define the following probability:
\[
p(n,m) = \left\{
\begin{array}{l}
\text{Probability that one of the $n$ candidates} \\
\text{wins a majority of the $2m+1$ votes.}
\end{array}
\right\}
\]Let’s find a formula for $p(n,m)$. For a particular candidate to win $k$ votes, $k$ members must vote for the candidate and the remaining members must vote for one of the other $n-1$ candidates. This occurs with probability $\binom{2m+1}{k}\left(\frac{1}{n}\right)^k \left(1-\frac{1}{n}\right)^{2m+1-k}$. The candidate wins a majority if they garner at least $m+1$ of the $2m+1$ votes. We multiply the final result by $n$ since we are counting the probability that any of the $n$ candidates wins a majority. Therefore,
\[
p(n,m) = n\sum_{k=m+1}^{2m+1} \binom{2m+1}{k}\left(\frac{1}{n}\right)^k \left(1-\frac{1}{n}\right)^{2m+1-k}
\]Note: We can add the probabilities of the different candidates winning as we did because the events are disjoint; it is impossible for two candidates to simultaneously win a majority of the votes. If the criteria allowed for overlap, e.g., one candidate winning at least 40% of the vote, then two candidates could win simultaneously, and we would have to account for this case to avoid double-counting.

Making the change of variables $k \mapsto 2m+1-k$, we obtain

$\displaystyle
p(n,m) = n \sum_{k=0}^m \binom{2m+1}{k} \left(1-\frac{1}{n}\right)^{k}\left(\frac{1}{n}\right)^{2m+1-k}
$

We recognize this as $n$ times the CDF of a binomial distribution with $2m+1$ trials and probability $1-\frac{1}{n}$.

What we’re really after is the expected number of rounds required for a candidate to secure a majority of the vote.
\[
f(n,m) = \left\{
\begin{array}{l}
\text{Expected number of rounds until one of the } \\
\text{that one of the $n$ candidates is elected.}
\end{array}
\right\}
\]This will happen in one round with probability $p(n,m)$. It will happen in $k$ rounds if it does not happen in the first $k-1$ rounds, and happens in the $k^\text{th}$ round. Putting this together, we get:

$\displaystyle
f(n,m) = \sum_{k=1}^{n-1} k\cdot p(n-k+1,m) \prod_{i=0}^{k-1} \bigl( 1-p(n-k+1+i,m) \bigr)
$

Unfortunately, these formulas do not simplify in any meaningful way, so we must evaluate them numerically. The first question asks us to compute $f(3,110)$ and the extra credit asks for $f(10,110)$. Calculating these, we obtain:
\begin{align}
f(3,110) &\approx 1.9999995110943353682 \\
f(10,110) &\approx 8.9999995110943296464
\end{align}These are rational numbers, of course, but $f(3,110)$ and $f(10,110)$ have 105-digit and 1357-digit denominators, respectively, when expressed as irreducible fractions, so I will not write them out here! One thing for sure is that it is very unlikely that the voting will end early. It is a virtual certainty that it will go the full $n-1$ rounds until there are only two candidates remaining.

Bonus: Changing the number of voters

The reason it is so difficult to elect a Speaker early is that there are so many voting members, and therefore it is unlikely that one candidate will obtain a majority of the votes by random chance. We can see how this probability would scale if we changed the number of voting members. Here is a plot of $f(n,m)$ for different values of $n$ (number of candidates), as a function of $2m+1$ (number of voting members). In all cases, we see that the expected number of rounds rapidly approaches $n-1$.

The slow car chase

This week’s Fiddler is a problem about probability and a slow car chase!

You and your pursuer are stopped at traffic lights one block away from each other, as shown below. For both cars, it takes 1 minute to get to each subsequent light, and there is a 50-50 chance any light will be red upon arrival (independent of what happened previously). If a light is red, you must wait one minute before it turns green. What is the probability you will make it to the city limits without being caught by your pursuer?

Extra Credit
In the same scenario as before, imagine there are infinitely many lights. On average, how many minutes will it be until you are ultimately caught?

My solution:
[Show Solution]

Can we escape before getting caught?

Let $P(n)$ be the probability that our pursuer will catch us at the $n^\text{th}$ light. To compute this probability, consider the gap between both cars. Initially, this gap is 1, and every time we reach a new light, depending on whether we got a red light and had to wait, or whether our pursuer got a red light and had to wait, the gap can change. Specifically,

The gap does not change if we both get red lights or both get green lights, which occurs with probability $\frac{1}{4}+\frac{1}{4} = \frac{1}{2}$.
The gap decreases by one if we get a red light and our pursuer gets a green light, which occurs with probability $\tfrac{1}{4}$.
The gap increases by one if we get a green light and our pursuer gets a red light, which occurs with probability $\tfrac{1}{4}$.

Initially, both cars have green lights, and since we get caught on the $n^\text{th}$ light, our final light must be red and our pursuer’s final light must be green. This leaves $n-2$ intermediate lights, at which the distance can either increase, decrease, or stay the same. If the distance could only increase or decrease, the total number of ways would be given by $C_{(n-2)/2}$, where $C_n$ is the $n^\text{th}$ Catalan number (see the interpretation of Catalan numbers as “counting mountain ranges”).

In our case, the gap can also stay constant, which complicates matters slightly. If we assume there are $k$ “ups” and $k$ “downs”, this leaves $n-2-2k$ “flats”. It follows that the total number of mountain ranges that may include plateaus is given by:
\[
\sum_{k=0}^{\left\lfloor \tfrac{n-2}{2} \right\rfloor} C_k \binom{n-2}{2k}
\]Each such path has a probability of occurring that depends on the number of flats, namely, it is
\[
\left(\frac{1}{4}\right)^{2k+1} \left(\frac{1}{2}\right)^{n-2-2k}
\]The reason for the $2k+1$ in the first term is that we must account for the “down” at the very end of the path. Putting everything together and substituting the formula for the Catalan number in, we obtain:
\[
P(n) = \sum_{k=0}^{\left\lfloor \tfrac{n-2}{2} \right\rfloor}
\frac{1}{k+1}\binom{2k}{k} \binom{n-2}{2k} \left(\frac{1}{4}\right)^{2k+1} \left(\frac{1}{2}\right)^{n-2-2k}
\]Amazingly, this sum has a closed-form formula, which I found using Mathematica. The probability that our pursuer will catch at the $n^\text{th}$ light is given by:

$\displaystyle
P(n) = \begin{cases}
\displaystyle\frac{1}{4^n\bigl(n-\tfrac{1}{2}\bigr)} \binom{2n}{n} & \text{for }n=2,3,4,\dots \\[1mm]
0 & \text{otherwise}
\end{cases}
$

We can double check that this is a legitimate probability distribution by verifying that $\sum_{n=2}^\infty P(n) = 1$. In other words, we will eventually be caught, it’s just a matter of when. I used Mathematica to evaluate the infinite sum, and it is indeed equal to one!

If we have $n$ lights to go before we reach the city limit, the probability we will escape before getting caught is $1-\sum_{k=2}^n P(k)$, or in other words:
\[
\text{Prob}(\text{escape}) = 1-\sum_{k=2}^n \frac{1}{4^k\bigl(k-\tfrac{1}{2}\bigr)} \binom{2k}{k}
\]Amazingly, this also has a closed-form expression… and thanks again to Mathematica, it is given by:

$\displaystyle
\text{Prob}(\text{escape}) = \begin{cases}
\displaystyle\frac{2}{4^n} \binom{2n}{n} & \text{for }n=2,3,4,\dots \\[1mm]
0 & \text{otherwise}
\end{cases}
$

This is an embarrassingly simple-looking expression, so there is likely a more elegant solution waiting to be found. The solution includes a denominator $4^n$ which is the total number of ways of assigning $\{\text{red},\text{green}\}$ to each of the $2n$ lights. The numerator is $2\binom{2n}{n}$, which must be the number of light assignments that do not lead to capture. If you can figure this out, please let me know in the comments!

In the original problem, we had 5 lights ahead of us, so substituting $n=5$, we obtain a probability of escape of $\frac{63}{128}$, or approximately 49.2%.

We can also check whether our formula holds up numerically. In the plot below, I compare the empirical probability from running 1,000,000 Monte Carlo simulations, and the analytic formula produced above, and we observe that we have a perfect match.

Note: This problem is somewhat similar to an old Riddler problem from 2017 about a horse race. You can check out my write-up for that problem if you’re interested!

Extra credit: time until caught

It takes at least one minute per light to visit all the lights (possibly more, depending on how many times we stop for red lights). Therefore, whenever we are eventually caught, the number of minutes that elapsed must be at least as much as the number of lights we saw.

We are asked to calculate the expected time elapsed until capture. Based on the argument above, we must have:
\[
\mathbb{E}(\text{lights seen until capture}) \leq
\mathbb{E}(\text{time until capture})
\]Let’s calculate the left-hand side. We already calculated $P(n)$, which is the the probability that we will be captured after seeing exactly $n$ lights. Therefore, the left-hand side is equal to $\sum_{n=2}^\infty n P(n)$. It turns out that this sum does not converge!. To see why, we can use the following asymptotic formula (source), which holds as $n\to\infty$.
\[
\binom{2n}{n} \sim \frac{4^n}{\sqrt{\pi n}}
\]This leads us to the approximation
\[
P(n) \sim \frac{1}{\bigl(n-\tfrac{1}{2}\bigr)\sqrt{\pi n}}
\]And now we see why there is no convergence… The function $n P(n)$ decays like $\frac{1}{\sqrt{n}}$, and therefore is not summable (it is a $p$-series with $p=\tfrac{1}{2}$). Since the expected number of lights is unbounded, then the time until capture must also be unbounded.

But does this make any sense? A probability distribution that has no average? Yes! In general, distributions for which some of the moments are unbounded are called heavy-tailed distributions, and there are many well-known examples. In particular, examples of distributions with unbounded first moments (average does not exist) include the Cauchy distribution, the Pareto distribution, and $P(n)$ defined above.

Braille puzzle

This week’s Fiddler is a counting problem about the Braille system.

Braille characters are formed by raised dots arranged in a braille cell, a three-by-two array. With six locations for dots, each of which is raised or unraised, there are 26, or 64, potential braille characters.

Each of the 26 letters of the basic Latin alphabet (shown below) has its own distinct arrangement of dots in the braille cell. What’s more, while some arrangements of raised dots are rotations or reflections of each other (e.g., E and I, R and W, etc.), no two letters are translations of each other. For example, only one letter (A) consists of a single dot – if there were a second such letter, A and this letter would be translations of each other.

Of the 64 total potential braille characters, how many are in the largest set such that no two characters consist of raised dot patterns that are translations of each other?

Extra Credit In addition to six-dot braille, there’s also an eight-dot version. But what if there were even more dots? Let’s quadruple the challenge by doubling the size of the cell in each dimension. Consider a six-by-four array, where a raised dot could appear at each location in the array. Of the 224 total potential characters, how many are in the largest set such that no two characters are translations of each other?

My solution:
[Show Solution]

We will solve a general version of this problem, where the Braille characters are $m\times n$. That is, there are $m$ rows and $n$ columns.

When two characters are translations of one another, we will say they are equivalent. The equivalence class of a character is the set of all characters equivalent to that particular one. For example, in the $3\times 2$ case, a single dot has six equivalents, and two vertically stacked dots have 4 equivalents. In the figure below, each row shows equivalent characters.

The size of an equivalence class is determined by the footprint of its characters. That is, the smallest rectangle that contains all the dots. For example, in the $6\times 4$ case below, the character shown has a $4\times 4$ footprint (the red outline):

If we translate this character, we see that it can occupy $3$ different vertical locations and $2$ different horizontal locations, for a total of $3\times 2 = 6$ equivalent characters. In general, if the footprint of a character is $i\times j$, then its equivalence class contains $(m-i+1)(n-j+1)$ elements.
Let us define:

$f(m,n)$: The size of the largest set of $m\times n$ characters such that no two characters are translations of one another. In other words, the number of different equivalence classes.
$g(i,j)$: The number of different characters that have a footprint $i\times j$.

The problem is asking us to calculate $f(m,n)$, but based on the discussion above, it is clear that:
\[
f(m,n) = \sum_{i=1}^m\sum_{j=1}^n g(i,j)
\]Computing $g(i,j)$ is fundamentally a (tedious) counting problem: Given a $i\times j$ rectangle, how many dot patterns can we make that fit snugly in that rectangle? The size of the footprint is determined by the dots on the border, so we separate the corners, edges, and middle section of the footprint:

Depending on which corners are filled in, this determines what the edges can look like so that the character fit snugly in the $i\times j$ rectangle. In particular, we have the following five cases:

No corners filled in. There must therefore be at least one dot on each edge. The vertical and horizontal edges contain $i-2$ and $j-2$ slots, respectively. Since we want to exclude the “empty” edge, we can fill in the edges in $\bigl(2^{i-2}-1\bigr)^2\bigl(2^{j-2}-1\bigr)^2$ ways.
One corner filled in. The adjacent edges can be filled in arbitrarily, while the two opposite edges must each contain a dot. Thus we can fill the edges in $\bigl(2^{i-2}-1\bigr)2^{i-2}\bigl(2^{j-2}-1\bigr)2^{j-2}$ ways. This can happen in 4 different ways (depending on which corner is filled).
Two bottom (or top) corners filled in. The three adjacent edges can be filled in arbitrarily, while the opposite vertical edge must contain a dot. Thus we can fill the edges in $\bigl(2^{i-2}\bigr)^2\bigl(2^{j-2}-1\bigr)2^{j-2}$ ways. This can happen in 2 different ways (we can pick the bottom or top corners).
Two left (or right) corners filled in. The three adjacent edges can be filled in arbitrarily, while the opposite horizontal edge must contain a dot. Thus we can fill the edges in $\bigl(2^{i-2}-1\bigr)2^{i-2}\bigl(2^{j-2}\bigr)^2$ ways. This can happen in 2 different ways (we can pick the left or right corners).
Two opposite corners filled in (or more). In this case, all four edges can be filled in arbitrarily since opposite corners completely determine the footprint. We therefore count $\bigl(2^{i-2}\bigr)^2\bigl(2^{j-2}\bigr)^2$ ways. This can happen in 7 different ways, since we can have only the opposite corners (2 ways), any three corners (4 ways), or all four corners (1 way).

For each of these cases, we can fill in the central $(i-2)\times(j-2)$ portion of the character arbitrarily, which means we should add all the above cases together, and multiply the total by $2^{(i-2)(j-2)}$.

Putting all of this together, we get our final tally:
\begin{multline}
g(i,j) = 2^{(i-2)(j-2)}\biggl(\bigl(2^{i-2}-1\bigr)^2\bigl(2^{j-2}-1\bigr)^2 \\
+ 4 \cdot \bigl(2^{i-2}-1\bigr)2^{i-2}\bigl(2^{j-2}-1\bigr)2^{j-2}
+ 2 \cdot \bigl(2^{i-2}\bigr)^2\bigl(2^{j-2}-1\bigr)2^{j-2} \\
+ 2 \cdot \bigl(2^{i-2}-1\bigr)2^{i-2}\bigl(2^{j-2}\bigr)^2
+ 7 \cdot \bigl(2^{i-2}\bigr)^2\bigl(2^{j-2}\bigr)^2 \biggr)
\end{multline}One final wrinkle: The above formula only makes sense if $i\geq 2$ and $j\geq 2$. Otherwise, our footprint will be a degenerate rectangle (it will not have four corners). Accounting for this, we obtain the complete version of $g(i,j)$:
\[
g(i,j) = \begin{cases}
\text{(formula above)} & i\geq 2\text{ and }j\geq 2 \\
2^{i-2} & i\geq 2\text{ and }j=1\\
2^{j-2} & i=1\text{ and }j\geq 2\\
1 & i=j=1
\end{cases}
\]To verify that the formula is correct, we can count the total number of characters by separating by footprint:
\begin{align}
& \text{(total $m\times n$ characters)} \\
&= \sum_{i=1}^m\sum_{j=1}^n \text{(characters with an $i\times j$ footprint)}\cdot \text{(size of $i\times j$ equiv. class)}\\
&= \sum_{i=1}^m\sum_{j=1}^n g(i,j) \cdot (m-i+1)(n-j+1) \\
&= 2^{mn}-1
\end{align}The last line was evaluated with Mathematica, and it actually worked! The reason for the “$-1$” is that we are excluding the “all-empty” character in our count.

Now that we have verified our formula for $g(i,j)$, we can substitute it in to the formula for $f(m,n)$ we derived earlier:
\[
f(m,n) = \sum_{i=1}^m\sum_{j=1}^n g(i,j)
\]Again, I relied on Mathematica to derive the analytic expression for this sum, and I obtained the following formula:

$\displaystyle
f(m,n) = \bigl(2^{m-1}-1\bigr) \bigl(2^{n-1}-1\bigr) 2^{(m-1) (n-1)}+2^{m n-1}
$

Given the relative simplicity of this final answer, I’m expecting there to be a much more elegant way to solve this problem than the approach I used. Let me know in the comments!

Solutions to special cases

Now that we have a general formula, we need only substitute the appropriate values for $m$ and $n$ to solve the problem. For the first problem, we find:
\[
f(3,2) = 44
\]Here is a picture that shows all $44$ characters. For each equivalence class, I chose the representative in the lower-left corner of the grid. I also indicated in red the number of elements in each equivalence class. As expected, if you sum all the red numbers you obtain $63 = 2^{6}-1$, which is the total number of non-empty characters.

For the extra credit problem, we find:
\[
f(6,4) = 15,\!499,\!264
\]

Making something out of nothing

This week’s Fiddler is a problem about composing functions. Here it goes:

Consider $f(n) = 2n+1$ and $g(n) = 4n$. It’s possible to produce different whole numbers by applying combinations of $f$ and $g$ to $0$. How many whole numbers between $1$ and $1024$ (including $1$ and $1024$) can you produce by applying some combination of $f$’s and $g$’s to the number $0$?

Extra Credit: Now consider the functions $g(n) = 4n$ and $h(n) = 1−2n$. How many integers between $-1024$ and $1024$ (including $-1024$ and $1024$) can you produce by applying some combination of $g$’s and $h$’s to the number $0$?

My solution:
[Show Solution]

Initially, we must start by applying $f$, otherwise we would just stay at zero. So for all intents and purposes, we can assume we start at $1$. At this point, it helps to write numbers out in binary (base two).

Applying $f(n)=2n+1$ has the net effect of tacking on a “1” to the end of the number in binary form, For example, if we start at $12$ (which is $1100$ in binary), then $f(12)=25$, which is $11001$ in binary.
Applying $g(n)=4n$ has the net effect of tacking on a “00” to the end of the number in binary form, so for example, if we start at $3$ (which is $11$ in binary), then $g(3)=12$, which is $1100$ in binary.

Now back to the original problem. Our goal is to count how many numbers between $1$ and $1024$ we can reach using compositions of $f$ and $g$. Since $1024=2^{10}$, let’s consider the more general problem of asking how many numbers between $1$ and $2^n$ we can reach. Let’s call this number $a_n$. Let’s look at the reachable numbers that have $k$ digits, and organize them in a table (in base 2):
\[
\begin{array}{c|l}
\text{digits} & \text{reachable numbers} \\ \hline
1 & 1 \\
2 & 11 \\
3 & 111, 100\\
4 & 1111, 1001, 1100\\
5 & 11111, 10011, 11001, 11100, 10000\\
\end{array}
\]The number of entries in row $k$ is simply $F_k$, the $k^\text{th}$ Fibonacci number. To see why, observe that to get to row $k$, we can start with the entries from row $k-1$ and apply $f$ (add a “$1$” to the end), or we can start with the entries from row $k-2$ and apply $g$ (add a “$00$” to the end). Since the first two rows each have one entry, we get the standard Fibonacci sequence (1,1,2,3,5,8,13,…).

Counting the reachable numbers up to $2^n-1$ means including all numbers with up to $n$ digits, so summing $F_1+F_2+\cdots+F_{n}$. However, we also want to include $2^n$. It’s clear from the pattern that this number is reachable whenever $n$ is even. Therefore, the reachable numbers $a_n$ is given by:
\[
a_n = F_1+\cdots+F_{n} + \begin{cases}
1 & n\text{ is even} \\
0 & n\text{ is odd}
\end{cases}
\]We can simplify this sum further. To see how, write:
\begin{align}
F_n &= F_{n+2}-F_{n+1}\\
F_{n-1} &= F_{n+1}-F_n \\
F_{n-2} &= F_n-F_{n-1} \\
\vdots\,\, &= \,\,\vdots\\
F_2 &= F_4-F_3\\
F_1 &= F_3-F_2
\end{align}Summing both sides, we get cancelations on the right-hand side:
\[
F_1+\cdots+F_{n} = F_{n+2}-F_2 = F_{n+2}-1
\]Therefore, we can simplify our formula and write:
\[
a_n = F_{n+2}-1+\frac{1}{2}\left( (-1)^n+1 \right)
\]Further simplification gets us our final answer:

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

We can also find a closed-form expression if we use Binet’s formula:

$\displaystyle
a_n = \frac{\phi^{n+2}-\psi^{n+2}}{\sqrt{5}}+\frac{(-1)^n-1}{2}
$

where $\phi=\frac{1+\sqrt{5}}{2}$ and $\psi=\frac{1-\sqrt{5}}{2}$ are the golden ratio and its conjugate, respectively. The first several values of $a_n$ are:
\[
\begin{array}{c|c|c|c|c|c|c|c|c|c|c|c}
n & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 \\ \hline
a_n & 1 & 1 & 3 & 4 & 8 & 12 & 21 & 33 & 55 & 88 & 144 \\
\end{array}
\]

The original problem asked about numbers up to $1024 = 2^{10}$, so we have:
\[
a_{10} = F_{12} = 144
\]

Extra credit

This problem is similar to the first part, but we have to take a bit more care with edge cases. In particular, the first part had the nice feature that applying $f$ always added one digit, and applying $g$ always added two digits (in base 2). This is (almost) still true. Applying $g$ still adds two digits, but applying $h$ does the following:

If the number is negative, it adds a one and flips the sign. For example: $h(-10)=21$ i.e., $-1010\mapsto 10101$.
If the number is positive, it adds a zero, flips the sign, then subtracts $1$. For example: $h(6)=-11$, i.e., $110\mapsto -1011$.

It can happen that applying $h$ will not increase the number of digits. For example, $h(4)=-7$, i.e., $100\mapsto -111$. Nevertheless, we can still make a table similar to how we did before, except we will not arrange numbers by digits. Instead, arrange them in such a way that the Fibonacci pattern remains true, like this:
\[
\begin{array}{c|l}
\text{row} & \text{reachable numbers} \\ \hline
0 & 1 \\
1 & -1 \\
2 & 11, 100 \\
3 & -101, -111, -100 \\
4 & 1011, 1111, 1001, 1100, 10000 \\
5 & -10101, -11101, -10001, -11001, -11111, -10100, -11100, -10000 \\
\end{array}
\]Note that the rows alternate positive and negative numbers. Also, row $k$ for $k\geq 1$ contains numbers with $k$ binary digits, except the even rows, which also contain the number $2^{k}$, which has $k+1$ digits. The reason for this discrepancy is so that we maintain the Fibonacci pattern. Now, we can see that row $k$ is formed by applying $h$ to row $k-1$, or applying $g$ to row $k-2$. Therefore, row $k$ contains $F_{k+1}$ numbers.

Now, we can see that summing rows $0$ through $k$ will give us reachable numbers between $-(2^{k}-1)$ and $2^{k}$. To make sure we include $-2^{k}$, we should also add one whenever $k$ is even. So far, we counted positive and negative reachable numbers, but what about zero? Since $g(0)=0$, zero is also reachable, so we should add one more to include it. Therefore, we find that:
\[
b_n = F_1+\cdots+F_{n+1} + 1 + \begin{cases}
1 & n\text{ is even} \\
0 & n\text{ is odd}
\end{cases}
\]Similarly to the previous case, we obtain the solution:

$\displaystyle
b_n = F_{n+3}+\frac{1}{2}\left( (-1)^n+1 \right)
$

or, using Binet’s formula,

$\displaystyle
b_n = \frac{\phi^{n+3}-\psi^{n+3}}{\sqrt{5}}+\frac{(-1)^n+1}{2}
$

The first several values of $b_n$ are:
\[
\begin{array}{c|c|c|c|c|c|c|c|c|c|c|c}
n & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 \\ \hline
b_n & 3 & 3 & 6 & 8 & 14 & 21 & 35 & 55 & 90 & 144 & 234 \\
\end{array}
\]The original problem asked about numbers between $-1024$ and $1024$, so since $1024=2^{10}$, we have:
\[
b_{10} = F_{13}+1 = 234
\]

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.

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