random walk – Book Proofs

Flipping your way to victory

This week’s Riddler Classic concerns a paradoxical coin-flipping game:

You have two fair coins, labeled A and B. When you flip coin A, you get 1 point if it comes up heads, but you lose 1 point if it comes up tails. Coin B is worth twice as much — when you flip coin B, you get 2 points if it comes up heads, but you lose 2 points if it comes up tails.

To play the game, you make a total of 100 flips. For each flip, you can choose either coin, and you know the outcomes of all the previous flips. In order to win, you must finish with a positive total score. In your eyes, finishing with 2 points is just as good as finishing with 200 points — any positive score is a win. (By the same token, finishing with 0 or −2 points is just as bad as finishing with −200 points.)

If you optimize your strategy, what percentage of games will you win? (Remember, one game consists of 100 coin flips.)

Extra credit: What if coin A isn’t fair (but coin B is still fair)? That is, if coin A comes up heads with probability p and you optimize your strategy, what percentage of games will you win?

Here is my solution:
[Show Solution]

At first, the problem appears to be asking the impossible. How can flipping fair coins over and over lead to anything other than winning half the time? The key here is that “winning” only requires getting a positive score. If we looked at the “expected total score” after 100 flips, this would be 0 no matter what strategy was used (you can see this by applying linearity of expectation). But since we’re counting number of positive outcomes rather than total score, it’s possible to win more than half of the time on average.

As a simpler example of how this can happen, consider a game where you roll a die. If you get a 1, you lose \$5. But if you get anything else, you win \$1. The expected score in the game is $\tfrac{1}{6}(-5) + \tfrac{5}{6}(+1) = 0$. However, if “winning” simply means obtaining a positive score, then we win $\tfrac{5}{6}$ of the time. The flipping game is similar; our total score can average out to zero, even though it ends up being positive more than half of the time. Practically speaking, this means that when we do lose, we lose big.

Probability of winning assuming optimal play

Let’s suppose we are doing $N$ flips, and we are choosing between:

Coin $A$: wins $+1$ with prob $a$ and wins $-1$ with prob $1-a$.
Coin $B$: wins $+2$ with prob $b$ and wins $-2$ with prob $1-b$.

Our task is to determine the percentage of games that we will win if we optimize our strategy. One way to approach this is using dynamic programming. The basic idea is that our best move at any point in the game should only depend on (1) our current score and (2) how many moves we have left. Therefore, we define:
\[
V_t(s) = \left\{\begin{aligned}
&\text{Probability that we win if we play optimally}\\
&\text{from here on out, given that our current score}\\
&\text{is }s\text{ and we are currently at turn }t\text{ of }N
\end{aligned}\right.
\]We can view each $V_t$ as a vector indexed by $s \in \{\dots,-1,0,1,\dots\}$, which is our current score. The idea is to solve for each $V_t$ by working our way backwards from the final step until we reach $V_0$, which is our probability of winning for each different starting score. Ultimately, the question is asking us to compute $V_0(0)$.

The terminal distribution is simply given by our probability of winning given that we have achieved a given final score. The result in this case is simple, since at that point there are no moves left to make:
\[
V_N(s) = \begin{cases}
0 & \text{if }s \leq 0\\
1 & \text{if }s > 0
\end{cases}
\]At each previous step, we pick the coin that will give us our highest chance of winning. This leads to the recursive formula for $t=n-1,\dots,2,1$:
\[
V_{t-1}(s) = \max\left\{ a V_t(s+1)+(1-a) V_t(s-1),\, b V_t(s+2)+(1-b) V_t(s-2) \right\}
\]I couldn’t find a way to evaluate this expression in closed form, but it is rather straightforward to evaluate it numerically. Here is efficient Python code that evaluates $V_0(0)$ for any choice of $(a,b,N)$.

import numpy as np

# probability of winning if coin A has prob a of +1 and (1-a) of -1,
# coin B has prob b of +2 and (1-b) of -2, and we do N flips in total
def compute_win_prob( a, b, N ):

    # shift the indices so that a score of "s" is at index v[s+2*N]    
    v = np.zeros(4*N+1)
    
    # initialize (time=N)
    v[:2*N] = 0
    v[2*N+1:] = 1
    
    # overwrite results as we recurse backwards to time=0
    # this code is efficient: uses vectorization and overwrites past results
    for t in range(N):
        coin_A = a*v[3:4*N] + (1-a)*v[1:4*N-2]
        coin_B = b*v[4:4*N+1] + (1-b)*v[0:4*N-3]
        v[2:4*N-1] = np.maximum(coin_A, coin_B)
    
    # return the probability of winning assuming we start with a score of 0
    return v[2*N]

When we run this script using $a=b=0.5$ and $N=100$, we obtain a probability of winning of about $0.6403$. We can get the exact answer as well using the following Mathematica code:

ComputeWinProb[a_,b_,n_]:=(
  v = ConstantArray[0,4n+1];
  v[[2n+2;;4n+1]] = 1;
  For[t=0, t<n, t++,
    CoinA = a v[[4;;4n]] + (1-a) v[[2;;4n-2]];
    CoinB = b v[[5;;4n+1]] + (1-b) v[[1;;4n-3]];
    v[[3;;4n-1]] = MapThread[Max,{CoinA,CoinB}];
  ];
  v[[2n+1]])
ComputeWinProb[1/2,1/2,100]

This gives us the exact result:
\[
\text{Prob}(\text{win})=\tfrac{811698796376000066208208781649}{1267650600228229401496703205376} \approx 0.640317
\]

Does doing more flips help?

We can win 64% of the time if we do 100 flips, but What if we did 1000 flips, or more? Could we keep increasing our chances of winning? It doesn’t appear to be the case. As we do more and more flips (increase $N$), the probability of winning increases slightly, but only up to a limit. To see this, I solved the problem for $N$ ranging up to 100,000 flips. Here is what I found:

Given how slowly the function is growing, it appears that there should be a finite limit, and that this limit is around $\frac{2}{3}$ (the last point on the plot is at $0.6658$). But is there another way we can compute this limit?

Limiting behavior for a large number of flips

The optimal strategy for the case $a=b=\tfrac{1}{2}$ is to flip coin A when our score is positive and to flip coin B when our score is negative. How often will we win the game as $N\to\infty$ if we use this strategy?

We can think of the coin-flipping problem in the context of random walks. It is well-known that random walks on the 1D line (with equal probability of moving left or right at each step) will visit every integer an infinite number of times, given enough time. It is also known that the expected number of steps before a walk first returns to its starting point is infinite. Therefore, as $N$ gets large, we can expect that the the score will spend only an infinitesimal fraction of the time near 0.

We can model this process as a Markov chain with states:
\[
\dots,-8,-6,-4,-2,0,1,2,3,4,\dots
\]When our score is $\leq 0$, we will flip coin B and move $\pm 2$. When our score is $\gt 0$, we will flip coin A and move $\pm 1$. Note that the states $-3,-5,-7,\dots$ can never be reached because we only move by $\pm 2$ when our score is negative. Based on the random walk argument above, we can expect the Markov chain to spend most of its time with a large-magnitude score (either on the positive or negative side).

As an approximation, we can imagine that once we hit $+2$, we will spend a large amount of time on the positive side (with a winning score) before we return to $+1$. Likewise, once we hit $0$, we will spend a large amount of time on the negative side (with a losing score) before we return. Since the positive and negative halves of this infinite Markov chain have identical probability distributions, we can expect the time spent on both halves to be equally large. We model this situation by truncating the Markov chain and assigning a transition probability of $\gamma \gg 0$ of remaining in a winning or losing state once we get there. Here is the resulting truncated Markov chain:

The transition matrix for this Markov chain is:
\[
A = \begin{bmatrix}
\gamma & \frac{1}{2} & 0 \\
0 & 0 & 1-\gamma \\
1-\gamma & \frac{1}{2} & \gamma
\end{bmatrix}
\]The steady-state distribution (right eigenvector corresponding to the eigenvalue of $1$) is given by: $\begin{bmatrix} \tfrac{1}{2} & 1-\gamma & 1\end{bmatrix}$. As anticipated, as $\gamma \to 1$, less and less probability is associated with the intermediate state 1. In the limit, the fraction of time spent in the winning state is
\[
\text{Prob}(\text{win as }N\to\infty) = \frac{1}{\tfrac{1}{2}+1} = \frac{2}{3}
\]So this justifies why we found the probability of winning tending to $\tfrac{2}{3}$ in the numerical simulation above.

As an aside, @electricvishnu with some help from @DarkAdonisSA found a formula that computes the probability of winning the game as a function of $N$. This formula is:

$\displaystyle
\text{Prob}(\text{win in }N\text{ steps}) = \frac{2}{3} + \frac{(-1)^N}{2^N}
\left[ \frac{1}{3} + \sum_{k=1}^{N-1} (-1)^k \binom{k}{\lfloor k/2 \rfloor} \right]
$

We don’t have a proof of formula as of yet, but it does produce the exact answer when $N=100$, and it also clearly reveals the limit of $\tfrac{2}{3}$ since the second term tends to zero as $N\to\infty$. There does not appear to be a meaningful way to simplify this formula either.

Changing the probabilities

We can also examine what happens when we try varying the probabilities away from $a=b=\tfrac{1}{2}$. Using the Python function shown above, this is a straightforward exercise. Here is what we obtain when we assume coin B is fair but coin A is not fair:

As expected, when coin A has a low probability, the optimal strategy is to flip coin B a lot more, which just leads to a probability of 0.5 of winning. As coin A’s probability increases, as do our chances of winning. We can also see what happens as we increase either $a$ or $b$ (neither coin is a fair coin). Here is what we obtain:

Once again, I struggled to find closed-form expressions for any of the above cases ($a$ and/or $b$ varying). Even for simpler instances, the solutions get complicated in a hurry. For example, the probability of winning as a function of $a$, with $b=\tfrac{1}{2}$, when the game only lasts 4 flips is:
\[
\text{Prob}(\text{win}) = \begin{cases}
\frac{1}{8}(a+4) & 0 \leq a \leq \tfrac{1}{2} \\
-\frac{1}{8} \left(16 a^4-32 a^3+20 a^2-11 a-1\right) & \tfrac{1}{2} \leq a \leq \gamma \\
-\frac{1}{2} a \left(6 a^3-11 a^2+6 a-3\right) & \gamma \leq a \leq 1
\end{cases}
\]where $\gamma \approx 0.7328$ is one of the roots of $8\gamma^3-4\gamma^2-1=0$. (this was found using the Mathematica function I included above). Each time we increase the number of flips, the number of parts in the piecewise polynomial can potentially double. So it seems highly unlikely that we’ll have a simple closed-form expression by the time we arrive at $N=100$!

Delirious ducks

This week’s Riddler Classic is about random walks on a lattice:

Two delirious ducks are having a difficult time finding each other in their pond. The pond happens to contain a 3×3 grid of rocks.

Every minute, each duck randomly swims, independently of the other duck, from one rock to a neighboring rock in the 3×3 grid — up, down, left or right, but not diagonally. So if a duck is at the middle rock, it will next swim to one of the four side rocks with probability 1/4. From a side rock, it will swim to one of the two adjacent corner rocks or back to the middle rock, each with probability 1/3. And from a corner rock, it will swim to one of the two adjacent side rocks with probability 1/2.

If the ducks both start at the middle rock, then on average, how long will it take until they’re at the same rock again? (Of course, there’s a 1/4 chance that they’ll swim in the same direction after the first minute, in which case it would only take one minute for them to be at the same rock again. But it could take much longer, if they happen to keep missing each other.)

Extra credit: What if there are three or more ducks? If they all start in the middle rock, on average, how long will it take until they are all at the same rock again?

Here is my solution:
[Show Solution]

One duck

We will solve this problem by starting with a simpler one and building up to the full problem. Let’s consider a single duck on the grid of rocks. At each step, the duck moves randomly to one of the other nearby rocks. We can describe this process using a Markov chain. Each state is a different rock, and the directed edges are labeled with the transition probabilities. Here is a diagram of what the Markov chain looks like:

If we number the states in the same way as above, then the transition matrix for this Markov chain is:
\[
P = \begin{bmatrix}
0 & \tfrac12 & 0 & \tfrac12 & 0 & 0 & 0 & 0 & 0 \\
\tfrac13 & 0 & \tfrac13 & 0 & \tfrac13 & 0 & 0 & 0 & 0 \\
0 & \tfrac12 & 0 & 0 & 0 & \tfrac12 & 0 & 0 & 0 \\
\tfrac13 & 0 & 0 & 0 & \tfrac13 & 0 & \tfrac13 & 0 & 0 \\
0 & \tfrac14 & 0 & \tfrac14 & 0 & \tfrac14 & 0 & \tfrac14 & 0 \\
0 & 0 & \tfrac13 & 0 & \tfrac13 & 0 & 0 & 0 & \tfrac13 \\
0 & 0 & 0 & \tfrac12 & 0 & 0 & 0 & \tfrac12 & 0 \\
0 & 0 & 0 & 0 & \tfrac13 & 0 & \tfrac13 & 0 & \tfrac13 \\
0 & 0 & 0 & 0 & 0 & \tfrac12 & 0 & \tfrac12 & 0 \\
\end{bmatrix}
\]The way to interpret this matrix is that $P_{ij}$ is the probability that we will transition from $i$ to $j$. This explains why the rows sum to 1 (the matrix is right-stochastic). We can express this fact neatly using matrix multiplication as: $P \mathbf{1} = \mathbf{0}$, where $\mathbf{0}$ and $\mathbf{1}$ are column vectors of all-zeros and all-ones, repectively. Mathematically, we can propagate probability distributions through this Markov chain via matrix multiplication. If we have some initial distribution $\mathbf{a}$ across the states, for example
\[
\mathbf{a}^\mathsf{T} = \begin{bmatrix} \tfrac12 & \tfrac12 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}
\]i.e. we are equally likely to be in states 1 or 2. Then at the next step, the probability distribution will transition $\mathbf{a}\to\mathbf{b}$ with
\[
\mathbf{b}^\mathsf{T} = \mathbf{a}^\mathsf{T}P = \begin{bmatrix} \tfrac16 & \tfrac14 & \tfrac16 & \tfrac14 & \tfrac16 & 0 & 0 & 0 & 0 \end{bmatrix}
\]So in the next step, we might be in states 1, 2, 3, 4, 5 with the associated probabilities above.

It’s important to distinguish states e.g. $\{1,2,\dots,9\}$ from distributions over states, which are the vectors of probabilities such as the $\mathbf{a}$ and $\mathbf{b}$ used above. In the case where a distribution across states is degenerate, i.e. concentrated entirely on a particular state $s$, then I’ll use the following notation to denote the corresponding distribution over states:
\[
\left( \mathbf{e}_s\right)_i = \begin{cases} 1 & \text{if } i=s \\ 0 & \text{otherwise} \end{cases}
\]So, for example, $\mathbf{e}_2^\mathsf{T} = \begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}$.

Stopping time

We are interested in the expected hitting time, which is the number of steps it will take on average to travel from some initial state $s$ to some target state $t \in \mathcal{T}$. To keep this general, I’m assuming there can be multiple target states, indicated by the set $\mathcal{T}$. Let’s define $\mathbf{q}_s$ to be the expected hitting time to get from state $s$ to any state in our terminal set $\mathcal{T}$. It turns out that $\mathbf{q}$ satisfies the recursive relationship:
\[
\mathbf{q}_s = \begin{cases}
0 & \text{if }s \in \mathcal{T} \\
1 + \sum_{i} P_{si} \mathbf{q}_i & \text{otherwise}
\end{cases}
\]The first case is clear: if we start in the terminal set, then we have already arrived, so the hitting time is zero. If we are outside the terminal set, then the expected hitting time will be $1$ plus the weighted sum of the hitting times of wherever we end up after the next transition. The recursion above is essentially the Bellman equation from dynamic programming.

Define the vector $\mathbf{t}_i = \begin{cases} 0 & \text{if }i \in \mathcal{T}\\ 1 & \text{otherwise} \end{cases}$.
We can rewrite the equation above in concise vector form as:
\[
\mathbf{q} = \textrm{diag}(\mathbf{t}) \left( \mathbf{1} + P \mathbf{q} \right)
\]Using the fact that $\textrm{diag}(\mathbf{t})\mathbf{1} = \mathbf{t}$, we can further simplify and obtain:
\[
\left( I-\textrm{diag}(\mathbf{t}) P \right) \mathbf{q} = \mathbf{t}
\]For the case of one duck, if we set the terminal set to be $\mathcal{T} = \{5\}$, then we substitute $\mathbf{t} = \begin{bmatrix} 1 & 1 & 1 & 1 & 0 & 1 & 1 & 1 & 1 \end{bmatrix}^\mathsf{T}$ above and obtain:
\[
\begin{bmatrix}
1 & -\tfrac12 & 0 & -\tfrac12 & 0 & 0 & 0 & 0 & 0 \\
-\tfrac13 & 1 & -\tfrac13 & 0 & -\tfrac13 & 0 & 0 & 0 & 0 \\
0 & -\tfrac12 & 1 & 0 & 0 & -\tfrac12 & 0 & 0 & 0 \\
-\tfrac13 & 0 & 0 & 1 & -\tfrac13 & 0 & -\tfrac13 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & -\tfrac13 & 0 & -\tfrac13 & 1 & 0 & 0 & -\tfrac13 \\
0 & 0 & 0 & -\tfrac12 & 0 & 0 & 1 & -\tfrac12 & 0 \\
0 & 0 & 0 & 0 & -\tfrac13 & 0 & -\tfrac13 & 1 & -\tfrac13 \\
0 & 0 & 0 & 0 & 0 & -\tfrac12 & 0 & -\tfrac12 & 1 \\
\end{bmatrix}
\begin{bmatrix}\mathbf{q}_1\vphantom{\tfrac12} \\ \mathbf{q}_2\vphantom{\tfrac12} \\ \mathbf{q}_3\vphantom{\tfrac12} \\
\mathbf{q}_4\vphantom{\tfrac12} \\ \mathbf{q}_5\vphantom{\tfrac12} \\ \mathbf{q}_6\vphantom{\tfrac12} \\
\mathbf{q}_7\vphantom{\tfrac12} \\ \mathbf{q}_8\vphantom{\tfrac12} \\ \mathbf{q}_9\vphantom{\tfrac12} \end{bmatrix}
= \begin{bmatrix}
1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\
0\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12}
\end{bmatrix}
\]Inverting this matrix and solving for $\mathbf{q}$, here is the result in array form:
\[
\begin{array}{|c|c|c|} \hline
\mathbf{q}_1 = 6 & \mathbf{q}_2 = 5 & \mathbf{q}_3 = 6 \\ \hline
\mathbf{q}_4 = 5 & \mathbf{q}_5 = 0 & \mathbf{q}_6 = 5 \\ \hline
\mathbf{q}_7 = 6 & \mathbf{q}_8 = 5 & \mathbf{q}_9 = 6 \\ \hline
\end{array}
\]So if we start at node 1, it will take us on average 6 moves to reach node 5. It is by no means obvious that the expected hitting times should be integers, since they represent an average path length over an infinite number of possible paths. They just turn out to be integers in this case. The symmetry observed here also makes sense, and we could leverage it to simplify the problem to a Markov chain with only 3 states, but I’ll leave that as an exercise to the reader! (the reason I didn’t do this from the start is that I want to generalize easily to the multi-duck case, and the symmetries there are more complicated).

What about the variant where we start at the origin but we want to know how much time it will take us to return to the origin? In this case, we simply start counting from our second move, and we add 1 to the answer to account for the move we skipped. If we start at node 5, i.e. our initial distribution is $\mathbf{e}_5^\mathsf{T}$, then our distribution on the second turn is $\mathbf{e}_5^\mathsf{T} P$, and so the expected number of moves to start at 5 and return to 5 is: $\mathbf{e}_5^\mathsf{T} P \mathbf{q} + 1 = 6$.

Two ducks

At first glance, the two-duck version might seem much more challenging than the on-duck version. As we shall see, the problem is essentially the same once we look at it the right way. The key is to imagine a Markov chain where instead of the states being $\{1,2,\dots,9\}$, the states are $\{(1,1), (1,2), \dots, (9,9)\}$. In other words, there are 81 states, consisting of all possible pairs $(s_1,s_2)$, where $s_i \in \{1,2,\dots,9\}$ is the position of duck $i$.

What would be the transition matrix for such a Markov chain? The transition probability $(a_1, a_2) \to (b_1, b_2)$ is simply $P_{a_1,b_1} P_{a_2,b_2}$, i.e. the product of the respective transition probabilities for each duck. This means that if we order our states according to the first duck first (lexicographical ordering), then the transition matrix for this augmented Markov chain will be $P\otimes P$, where $\otimes$ is the Kronecker product.

So we can solve the problem in a very similar fashion. This time, our terminal set consists of all the nodes $(s_1, s_2)$ where $s_1=s_2$. There are 9 such nodes, and we can form the associated $\mathbf{t}$ by evaluating $\mathbf{t} = \mathbf{1}-\textrm{vec}(I)$, where $\mathrm{vec}$ is the vectorized version of the matrix, which you obtain by enumerating the columns. For example:
\[
\mathrm{vec} \begin{bmatrix}a_{11} & a_{12} & a_{13} \\
a_{21} & a_{22} & a_{23} \\
a_{31} & a_{32} & a_{33}
\end{bmatrix} =
\begin{bmatrix}
a_{11} \\ a_{21} \\ a_{31} \\ a_{12} \\ a_{22} \\ a_{32} \\ a_{13} \\ a_{23} \\ a_{33} \end{bmatrix}
\]the new equation for $\mathcal{q}$ is $\bigl( I-\textrm{diag}(\mathbf{t}) (P\otimes P) \bigr) \mathbf{q} = \mathbf{t}$ and the solution to the problem where you start at $(5,5)$ and want to know the expect time it takes for us to reach the terminal set (where the two ducks are on the same rock) can be computed following the same logic as in the one-duck case. That is, $\mathbf{e}_{(5,5)}^\mathsf{T} (P\otimes P) \mathbf{q} + 1$.

A cautionary note

While it may seem like this is straightforward algebra, there is a major issue I neglected to discuss: invertibility. When we solve the equation for $\mathbf{q}$ in the one-duck case, the equation looks like $A\mathbf{q} = \mathbf{t}$, and the solution is $\mathbf{q} = A^{-1}\mathbf{t}$. But what if $A$ isn’t invertible? This corresponds to cases where the Markov chain is not connected. So if you start on one island and the terminal set is on another island, you’ll never get there; the stopping time is infinite. That doesn’t happen with one duck, but it does happen with two ducks! This is because every time a duck moves to a new stone, the parity of the stone changes from odd to even or vice versa. So if one duck started on rock 1 and the other started on rock 2, they would never meet.

From a linear algebra standpoint, it just means that even though the matrix is not invertible, if we restrict our attention to the “island” (read: subspace) where the two ducks have the same parity, then the associated transition matrix for that island will be invertible and everything will be fine. From a practical standpoint, the equation above can still be solved; we just need to replace the inverse by a pseudoinverse, and ignore all the components of $\mathbf{q}$ for which the two components have different parity. More on this later…

Computation

I computed the result using the following Matlab code:

% transition matrix for one duck
P = [ 0 1/2 0 1/2 0  0  0  0  0
     1/3 0 1/3 0 1/3 0  0  0  0
      0 1/2 0  0  0 1/2 0  0  0
     1/3 0  0  0 1/3 0 1/3 0  0
      0 1/4 0 1/4 0 1/4 0 1/4 0
      0  0 1/3 0 1/3 0  0  0 1/3
      0  0  0 1/2 0  0  0 1/2 0
      0  0  0  0 1/3 0 1/3 0 1/3
      0  0  0  0  0 1/2 0 1/2 0 ];

T = kron(P,P);     % transition matrix
t = 1-vec(eye(9)); % terminal states

% starting distribution
s = kron([0 0 0 0 1 0 0 0 0],[0 0 0 0 1 0 0 0 0])';

% compute stopping time as a rational number
stop_time = s'*T*( (eye(81)-diag(t)*T)\t ) + 1
rats(stop_time)

The resulting expected stopping time was $\frac{363}{74}$, or approximately $4.905$ steps. We can also find an approximate answer via simulation. Here is what you get by simulating one million trials:

The staggered decrease of the probabilities is a real effect; it’s not a result of approximation error!

Many ducks

If we have $n$ ducks, we can clearly just generalize the approach used above. This time, the transition matrix will be $\underbrace{P\otimes \cdots \otimes P}_{n\text{ times}}$. The issue here is that our transition matrix will be quite large: $9^n \times 9^n$, to be precise. One way to shrink the number of states is to return to the idea of parity. The reason our transition matrix is large is because it describes all transitions between all possible duck configurations. Since we know all ducks start on the same node and that parity will be preserved, we only need to worry about the subset of states for which all ducks have the same parity. Namely, there will be $5^n$ odd states and $4^n$ even states. So we can restrict ourselves to a smaller $(5^n+4^n)\times(5^n+4^n)$ transition matrix. This is still large, but it’s an improvement. We will also leverage the fact that the transition matrix will be sparse, which will help with computation.

After all this work, we end up with the following result (for up to 6 ducks):

The values are as follows:

Number of ducks	Expected rendez-vous time
2	4.9054
3	18.4360
4	66.7420
5	237.3955
6	825.3364

It’s clear the expected rendez-vous time grows exponentially, but I haven’t found a practical way to approximate it or bound it. Unfortunately, it’s pretty difficult to extend my approach beyond $n=6$ ducks. To give you an idea of scale, our naive transition matrix would have had $9^6=531,441$ rows and columns. After our reduction procedure, we’re down to $5^6+4^6=19,721$. The transition matrix is sparse (about 98.5% of its entries are zeros), but that still leaves about 6 million nonzero entries, which makes computing $\mathbf{q}$ a challenge.

Further reductions are possible if we leverage symmetries. Let’s say we have 2 ducks. If we don’t do any reductions, there are $9^2=81$ states. If we reduce using even/odd parity, we drop to $5^2+4^2 = 41$.

From a stopping time standpoint, the only thing that matters is how many ducks are on each stone; it doesn’t matter which ducks are on which stones. This immediately allows us to reduce $5^n+4^n$ to ${n+4 \choose 4} + {n+3\choose 3}$, which is a dramatic improvement. In the two-duck case, this brings us to $25$.
Deeper symmetries occur as well. For example, a pair of ducks on 1 and 3 is equivalent to a pair of ducks on 7 and 9 by rotational symmetry. In the two-duck case, this brings us to $8$ states. Counting and bookkeeping these states might require some group theory…

Ultimately, these reductions would have a dramatic effect, reducing the scaling from exponential to polynomial. I didn’t implement any of these additional reductions because at this point (up to $n=6$) the exponential trend was clear. I didn’t think that extending the trend further would be particularly enlightening.

The lucky derby

In the spirit of the Kentucky Derby, this Riddler puzzle is about a peculiar type of horse race.

The bugle sounds, and 20 horses make their way to the starting gate for the first annual Lucky Derby. These horses, all trained at the mysterious Riddler Stables, are special. Each second, every Riddler-trained horse takes one step. Each step is exactly one meter long. But what these horses exhibit in precision, they lack in sense of direction. Most of the time, their steps are forward (toward the finish line) but the rest of the time they are backward (away from the finish line). As an avid fan of the Lucky Derby, you’ve done exhaustive research on these 20 competitors. You know that Horse One goes forward 52 percent of the time, Horse Two 54 percent of the time, Horse Three 56 percent, and so on, up to the favorite filly, Horse Twenty, who steps forward 90 percent of the time. The horses’ steps are taken independently of one another, and the finish line is 200 meters from the starting gate.

Handicap this race and place your bets! In other words, what are the odds (a percentage is fine) that each horse wins?

Here is my full derivation (long!):
[Show Solution]

Let’s take a look at a particular horse. Each step, suppose it moves forward with probability $q > \frac{1}{2}$ or backward with probability $1-q$. The horse starts at position $0$ and we would like to know the probability that it will reach position $2n$ in exactly $2k$ steps. The reason for the factor of $2$ is that to reach an even-numbered position (such as $200$, in the problem statement), the horse must take an even number of steps. Let’s call this probability $P(n,k,q)$.

This scenario is known as a 1D random walk, and past posts that deal with similar ideas include The Deadly Board Game and Gambler’s Ruin. This problem is a bit different, however, since we care about the distributions of the number of moves required (not just the expected number of moves).

Back to the task at hand: determining $P(n,k,q)$. If we reach the position $2n$ in $2k$ moves, then clearly we have $k \ge n$. Next, we can deduce that we must have moved forward $(k+n)$ times and backward $(k-n)$ times. The probability of this occurring is $q^{n+k}(1-q)^{n-k}$. The tricky part here is to deal with the order in which we take the steps. The last step we take must always be forward, so we will restrict our attention to arrange the remaining $2k-1$ steps. There is a total of $2k-1 \choose k-n$ ways to do this, but this is not the number we’re looking for! Some of these paths will reach $2n$ too early, so we will count these paths and subtract them.

The quantity we’re counting is closely related to the Catalan Numbers, and we can borrow a similar proof technique (namely André’s reflection method) to arrive at the solution. Here is the argument: imagine the path laid out as a “mountain range” where each of the $2k-1$ steps is a move to the right and either up or down depending on whether the step was forward or backward. The diagram below illustrates a valid path for $n=2$ and $k=5$. Note that we are allowed to go negative!

We would like to exclude “bad paths” that reach $n$ at some earlier point. Consider the first point at which this occurs, and reflect the remainder of the path about the line $y=2n$. This leads to a new path that has $(k+n)$ forward steps, $(k-n-1)$ backward steps, and ends at $2n+1$ instead of $2n-1$. See the figure below for an illustration of the reflection method.

There is a one-to-one correspondence between bad paths and unconstrained paths of length $2k-1$ that reach $2n+1$. So the total number of valid paths (total paths minus bad paths) is given by:
\[
{2k-1 \choose k-n}-{2k-1 \choose k-n-1}=\frac{n}{k}{2k \choose k-n}
\]The simplification follows from (these identities), but can also be easily derived using algebra. Putting everything together, the probability of a horse first reaching $2n$ in exactly $2k$ steps is given by the formula:

$\displaystyle
P(n,k,q) = \frac{n}{k}{2k \choose k-n} q^{k+n} (1-q)^{k-n},\quad\text{for }k=n,n+1,\dots
$

We can plot the distributions for our case of interest, which is $n=100$ and $q$ varies from $0.52$ (for Horse $1$) up to $0.90$ (for Horse $20$). Here is a plot of the distributions for some of the horses:

As we can see, the distributions are normal to very good approximation. Although they are challenging expressions to manipulate, the task is not beyond the capabilities of Mathematica! In fact, we can verify that $P$ is a legitimate probability mass function and we can also compute its mean and variance for the case $\tfrac{1}{2} \lt q \le 1$:
\begin{align}
\sum_{k=n}^\infty P(n,k,q) &= 1 & & \text{(sums to 1)}\\
\sum_{k=n}^\infty k\, P(n,k,q) &= \tfrac{n}{2q-1} & & \text{(computing the mean)}\\
\sum_{k=n}^\infty \bigl(k-\tfrac{n}{2q-1}\bigr)^2 P(n,k,q) &= \tfrac{2nq(1-q)}{(2q-1)^3} & & \text{(computing the variance)}
\end{align}Here is the mean and standard deviation for a few of the horses:
\begin{align}
\text{Horse 4} &\approx 1250 \pm 218\\
\text{Horse 8} &\approx 625 \pm 74\\
\text{Horse 12} &\approx 417 \pm 37\\
\text{Horse 16} &\approx 313 \pm 21\\
\text{Horse 20} &\approx 250 \pm 12
\end{align}Note: we must double the means and standard deviations from the formulas above since those formulas skip over the odd lengths (which are impossible). When making a continuous approximation, we want those spaces filled in!

One thing we notice right away is that the distributions are very well separated (it’s a log scale!). The only horses that stand any chance of beating Horse 20 are the horses that are close behind. Let’s plot the distributions one more time, on a linear scale, but only for the top horses:

If horse $i$ finishes in $h_i$ steps, computing the probability that horse $i$ wins amounts to finding the probability that $h_i > h_j$ for all $i\ne j$. We can approximate this as:
\begin{align}
\mathbb{P}(h_i > h_j\,\,\forall i \ne j)
&\approx \int_{-\infty}^\infty \mathbb{P}(h_i = x) \prod_{j \ne i} \mathbb{P}(h_j > x)\, \mathrm{d} x\\
&= \int_{-\infty}^\infty \frac{1}{\sigma_i}\varphi_i(z_i)\prod_{j \ne i} \left( 1-\Phi(z_j)\right) \, \mathrm{d} x\\
\end{align}where $\varphi$ and $\Phi$ are the probability density (pdf) and cumulative distribution (cdf) of the standard Normal distribution and $z_k = \tfrac{x-\mu_k}{\sigma_k}$ is the standard normal deviate computed for the $k^\text{th}$ horse ($\mu_k$ and $\sigma_k$ are the mean and standard deviations computed earlier).

This is a very messy integral for which there does not exist a closed-form solution, so we must resort to evaluating it numerically. For the solution, read on!

For the tl;dr, the answer is:
[Show Solution]

Horse number	Probability of winning the race
20	71.27%
19	21.54%
18	5.605%
17	1.268%
16	0.2595%
15	0.05085%
14	0.01018%
13	0.002212%

The deadly board game

This Riddler classic puzzle involves a combination of decision-making and probability.

While traveling in the Kingdom of Arbitraria, you are accused of a heinous crime. Arbitraria decides who’s guilty or innocent not through a court system, but a board game. It’s played on a simple board: a track with sequential spaces numbered from 0 to 1,000. The zero space is marked “start,” and your token is placed on it. You are handed a fair six-sided die and three coins. You are allowed to place the coins on three different (nonzero) spaces. Once placed, the coins may not be moved.

After placing the three coins, you roll the die and move your token forward the appropriate number of spaces. If, after moving the token, it lands on a space with a coin on it, you are freed. If not, you roll again and continue moving forward. If your token passes all three coins without landing on one, you are executed. On which three spaces should you place the coins to maximize your chances of survival?

Extra credit: Suppose there’s an additional rule that you cannot place the coins on adjacent spaces. What is the ideal placement now? What about the worst squares — where should you place your coins if you’re making a play for martyrdom?

Here is my solution:
[Show Solution]

Let’s suppose the track has length $N=1000$ (this will turn out not to matter!). We will begin by solving the simpler problem where we only get to place a single coin, and then we’ll solve the versions for two coins, three coins, etc.

One coin

Let’s call $p_k$ the probability that we will survive if the coin is placed on space $k$. In other words, $p_k$ is the probability that we land on the space labeled $k$. Clearly $p_1 = \tfrac{1}{6}$, but it gets more complicated for larger $k$. The key is to realize that we can compute $p_k$ recursively. For example, to compute some $p_k$, note that if our first roll is a one, we must now roll a total $k-1$ with our remaining rolls in order to survive. If our first roll is a two, we must roll a total of $k-2$ with our remaining rolls, and so on. So we conclude that:
\[
p_k = \frac{1}{6}\left( p_{k-1}+p_{k-2}+\dots+p_{k-6} \right)
\]This holds for all $k$, actually, so long as we choose our initial conditions appropriately. We know that $p_1 = \tfrac{1}{6}$, so we can achieve the desired recurrence by choosing:
\[
p_0 = 1,\quad\text{and}\quad p_{k} = 0\text{ for }k<0 \]So now we have a simple recursive way of computing all the $p_k$. In fact, we are dealing with a linear recurrence relation with constant coefficients, so we can find a formula for $p_k$ using the characteristic polynomial. The result turns out to be:
\[
p_k = \frac{1}{7}\left( 2 + \eta_1^k + \eta_2^k + \bar{\eta}_2^k + \eta_3^k + \bar{\eta}_3^k \right)
\]where $\eta_1\approx -0.670332$, $\eta_2\approx -0.375695 + 0.570175i$, and $\eta_3\approx 0.294195 + 0.668367i$. Here, $\bar{\eta}$ denotes the complex conjugate of $\eta$. Unfortunately we can’t find an exact expression because the characteristic polynomial is sixth-order, and even if we factor out the trivial root at 1, we are left with a quintic, which will not have a closed-form solution in general. Nevertheless, this formula tells us that the probability of landing on a coin placed very far away is $p_\infty = \frac{2}{7}$, Neat!

Ok, back to business. The formula for $p_k$ isn’t all that useful, but we still have an efficient way of computing $p_k$ exactly via the recursion. Here is a plot of what $p_k$ looks like:

Notice that in the limit, we approach $\tfrac{2}{7} \approx 0.285714$, as expected. We can read off the coin locations that either maximize survival, or minimize it (martyrdom). The results are shown in the table below.

Deadly Game with one coin	Optimal coin location	Exact survival probability	Approximate probability
Martyrdom	1	$\frac{1}{6}$	0.166667
Survival	6	$\frac{16807}{46656}$	0.360232

Two coins

Since the one-coin problem was difficult in the sense that we had to resort to computing all the probabilities and then picking the largest and smallest one, we can’t expect the two-coin problem to be any easier. That being said, we can leverage the work we’ve already done. Let’s call $q_{k,m}$ the probability that we will survive if the coins are placed on $k$ and $m$. We can compute $q$ by using the inclusion-exclusion principle. Namely,
\begin{align}
&P(\text{land on } k \text{ or } m) \\
&\qquad\qquad= P(\text{land on }k)+P(\text{land on }m)-P(\text{land on }k \text{ and } m)
\end{align}Landing on $k$ and $m$ requires landing on $k$ first (assuming $k\le m$) and then landing on $m$ given that we are on $k$. This is simply $p_{m-k}$. So in summary, if $k\le m$, we have:
\[
q_{k,m} = p_k+p_m-p_k p_{m-k}
\]Computing the $q$ values for each pair $(k,m)$, we can visualize the probabilities in 3D:

The landscape is even more complicated than before. The diagonal (corresponding to $k=m$) is sunken because this is when both coins are placed on the same square (much lower probability of survival!). The results are in the table below:

Deadly Game with two coins	Optimal coin locations	Exact survival probability	Approximate probability
Martyrdom	1, 2	$\frac{1}{3}$	0.333333
Survival	5, 6	$\frac{2401}{3888}$	0.617541

If we forbid adjacent coins, we obtain:

Martyrdom	1, 7	$\frac{16807}{46656}$	0.360232
Survival	4, 6	$\frac{4459}{7776}$	0.573431

Three coins

This case is similar to the two-coin case, except that the solution space is now three-dimensional and therefore much more difficult to visualize. If we let $r_{k,m,n}$ be the probability that we will survive if the coins are placed on $k$, $m$, and $n$, and we assume that $k\le m \le n$, we can compute the probability using inclusion-exclusion once more:
\[
r_{k,m,n} = p_k+p_m+p_n-p_kp_{m-k}-p_kp_{n-k}-p_mp_{n-m}+p_kp_{m-k}p_{n-m}
\]The results for this case are given in the table below.

Deadly Game with three coins	Optimal coin locations	Exact survival probability	Approximate probability
Martyrdom	1, 2, 7	$\frac{3697}{7776}$	0.475437
Survival	4, 5, 6	$\frac{343}{432}$	0.793981

If we forbid adjacent coins, we obtain:

Martyrdom	1, 3, 7	$\frac{3913}{7776}$	0.503215
Survival	6, 8, 10	$\frac{1198777}{1679616}$	0.713721

Note that it’s never in our best interest to place coins far apart from one another. All the “interesting” behavior (very large or very small probabilities) happens near the start. Take the three-coin case for example; if we place the coins far from the start and far from each other, each will get landed on with probability $\tfrac{2}{7}$ and they will be roughly independent. This means that:
\[
\textrm{P(survival)} \approx 1-(1-\tfrac{2}{7})^3 \approx 0.635569
\]and this places us squarely between the martyrdom and survival probabilities found in the table above. In other words, there is never any benefit to spreading out the coins, regardless of whether we are trying to win or lose… The $N=1000$ is a red herring!

Random walk with endpoints

The Riddler post for today is about a bar game where you flip a coin and move forward or backward depending on the result. Here is the problem:

Consider a hot new bar game. It’s played with a coin, between you and a friend, on a number line stretching from negative infinity to positive infinity. (It’s a very, very long bar.) You are assigned a winning number, the negative integer -X, and your friend is assigned his own winning number, a positive integer, +Y. A marker is placed at zero on the number line. Then the coin is repeatedly flipped. Every time the coin lands heads, the marker is moved one integer in a positive direction. Every time the coin lands tails, the marker moves one integer in a negative direction. You win if the coin reaches -X first, while your friend wins if the coin reaches +Y first. (Winner keeps the coin.)

How long can you expect to sit, flipping a coin, at the bar? Put another way, what is the expected number of coin flips in a complete game?

Here is my solution:
[Show Solution]

This problem is famously known as Gambler’s ruin and it’s an example of a one-dimensional random walk. In the typical Gambler’s ruin setup, you flip a coin repeatedly and either gain \$1 or lose \$1 depending on the outcome. You have a fixed starting budget and a desired goal amount. You keep playing until either you have reached your goal amount, or you run out of money. The bar game we’re looking at here is precisely Gambler’s ruin, but shifted so that you start with \$0 and end when you’re either up \$Y or you’re \$X in debt.

Here is one possible solution. Let’s define $E_n$ to be the expected number of coin flips in a complete game assuming the marker starts at location $n$ on the number line. We would like to find $E_0$. Suppose that at each flip, we move forward with probability $p$ and backward with probability $q$ (where $p+q=1$). If we are currently at $n$, then with probability $p$, we came from $n-1$ and with probability $q$ we came from $n+1$. Either way, we incurred one additional flip. We can therefore write the following recurrence relations for the expected number of coin flips:

\begin{align}
E_{-X} &= 0 \\
E_n &= p E_{n-1} + q E_{n+1} + 1 \qquad\text{for }-X < n < Y\\ E_Y &= 0 \end{align} In matrix form, these equations look like: \[ \begin{bmatrix} 1 & -q & 0 & \dots & 0 \\ -p & 1 & -q & \ddots & \vdots \\ 0 & -p & 1 & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & -q \\ 0 & \dots & 0 & -p & 1 \end{bmatrix} \begin{bmatrix} E_{-X+1} \\ E_{-X+2} \\ \vdots \\ E_{Y-2} \\ E_{Y-1} \end{bmatrix} = \begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \\ 1 \end{bmatrix} \] This system of equations has a nice structure; the matrix on the left is both tridiagonal and Toeplitz. One can solve this system of equations using a specialized form of Gaussian elimination, for example.

Although the equations have a closed-form solution, it is pretty messy. Instead, I’ll show how you can derive the answer by hand in the case $p=q=\tfrac{1}{2}$. The last equation can be solved for $E_{Y-1}$ in terms of $E_{Y-2}$:

\[
E_{Y-1} = \tfrac{1}{2} E_{Y-2} + 1
\]

We can then substitute this into the second equation and solve for $E_{Y-2}$ in terms of $E_{Y-3}$:

\[
E_{Y-2} = \tfrac{2}{3} E_{Y-3} + 2
\]

Substitute into the third equation and solve for $E_{Y-3}$ in terms of $E_{Y-4}$:

\[
E_{Y-3} = \tfrac{3}{4} E_{Y-4} + 3
\]

The pattern is now clear (you can prove it using induction if you like). Continuing in this fashion until we get to the starting point,

\[
E_{1} = (1- \tfrac{1}{Y}) E_{0} + Y – 1
\]

We can do something similar but starting from the other end. e.g. solving for $E_{-X+1}$ in terms of $E_{-X+2}$, and then $E_{-X+2}$ in terms of $E_{-X+3}$, and so on. This time, we obtain the result:

\[
E_{-1} = (1- \tfrac{1}{X}) E_0 + X – 1
\]

We can now combine the fruits of our labor with the only equation we haven’t used yet, which relates $E_0$, $E_1$, and $E_{-1}$. The result is that:

\begin{align}
E_0 &= \tfrac{1}{2} E_{-1} + \tfrac{1}{2} E_1 + 1 \\
&= \tfrac{1}{2}\left( (1- \tfrac{1}{X}) E_0 + X – 1 \right) + \tfrac{1}{2} \left( (1 – \tfrac{1}{Y}) E_0 + Y – 1 \right) + 1\\
&= \tfrac{1}{2}\left( 2- \tfrac{1}{X} – \tfrac{1}{Y} \right) E_0 + \tfrac{1}{2}(Y + X)
\end{align}

Rearranging and solving for $E_0$, we obtain:

$\displaystyle
E_0 = \frac{ Y + X }{ \tfrac{1}{X} + \tfrac{1}{Y} } = XY
$

A similar approach can be used to compute the probability that each player will win. To do this, let $P^Y_n$ be the probability that we’ll reach $Y$ first given that we are currently at $n$, and similarly define $P^X_n$. Then the recursions look like:

\begin{align}
P^Y_{-X} &= 0 \\
P^Y_n &= p P^Y_{n-1} + q P^Y_{n+1} \qquad\text{for }-X < n < Y\\ P^Y_Y &= 1 \end{align} Written as matrices, the only difference between this problem and that of finding expectation is the right-hand side of the equation! Using a similar recursive approach, we can solve the equations and we obtain:

$\displaystyle
P_0^X = \frac{Y}{X+Y} \qquad\text{and}\qquad P_0^Y = \frac{X}{X+Y}
$

So, practically speaking, if we must travel twice as far on the number line as our opponent, we are half as likely to win the game. Just for fun, I made some plots that show how the expected value and the probability of winning change if you’re not using a fair coin. The plots show the case where $X+Y = 50$. For example, the 40% curve means that the coin flips in your favor 40% of the time.

If you’re interested in a derivation of the solution for the case $p\ne \tfrac{1}{2}$, there is a well-written set of notes here.

and here is a much slicker solution, courtesy of Daniel Ross:
[Show Solution]