Finding the luckiest coin

I have in my possession 1 million fair coins. Before you ask, these are not legal tender. Among these, I want to find the “luckiest” coin.

I first flip all 1 million coins simultaneously (I’m great at multitasking like that), discarding any coins that come up tails. I flip all the coins that come up heads a second time, and I again discard any of these coins that come up tails. I repeat this process, over and over again. If at any point I am left with one coin, I declare that to be the “luckiest” coin.

But getting to one coin is no sure thing. For example, I might find myself with two coins, flip both of them and have both come up tails. Then I would have zero coins, never having had exactly one coin.

What is the probability that I will — at some point — have exactly one “luckiest” coin?

Let $N_t$ be the number of remaining coins after $t$ flips, removing any coins that were tails. We start with $N_0 = 1,000,000.$ We could set this entire problem up to be a Markov chain with binomial transition probabilities, however, since we will want to keep track of the “luckiest” we would be better served making both $N=0$ and $N=1$ into absorbing states. That is, let $P = \{ P_{ij} \}_{i,j=1}^{N_0}$ be the transition probability matrix then $$P_{ij} = \mathbb{P} \{ N_{t+1} = j \mid N_t = i \} = \begin{cases} 1, &\text{if $i=j=0$ or $i=j=1$;}\\ \binom{i}{j} 2^{-i}, &\text{if $j \geq 2$, $0 \leq i \leq j$;}\\ 0, &\text{otherwise.}\end{cases}$$

If we define the transient transition matrix $Q$ as the portion of the transition matrix associated with states $N \in \{2, \dots, N_0\},$ then the fundamental matrix of our Markov chain is given by $N=(I-Q)^{-1}.$ Further, letting $r = (P_{21},P_{31},\dots, P_{N_01})^T = (1/2, 3/2^3, \dots, N_0/2^{N_0})^T$, then the vector $b = \{b_i\}_{i=2}^{N_0}$ given by $b = Nr$ is gives the probability of starting at node $N=i \in \{2, \dots, N_0\}$ and ending up with a single luckiest coin. This is all nice in theory, but even though this just requires solving a triangular system of linear equations, given the space and time requirements of actually calculating the answer, it is somewhat prohibitive.

Instead, we can either (a) estimate the answer using simulations or (b) rely on estimates with some sufficiently large $N_0.$ Both of these methods give the same result of roughly $72\%$ that I will find a luckiest coin.

In [1]:

import numpy as np
from scipy.stats import binom

N0 = 1000000
def simulate(N0=N0):
    Niter = N0
    while Niter > 1:
        T = binom.rvs(Niter, 0.5, size=1)
        Niter -= T
    return Niter

num_sims = 10000
output = []
for i in range(num_sims):
    out = simulate()
    output.append(float(out))
print('{} simulations with N0={} yield an estimated probability of {}'.format(num_sims, N0, sum(output) / num_sims))


## Running matrix algebra with (N0-2) x (N0-2) matrices is dicey for N0 = 1e6, so let's
## scale it down and appeal to the fact that the answer should not change much for any N0
## as N0 -> inf
M = 5000
Q = np.eye(M-1)
for i in range(M-1):
    for j in range(i+1):
        Q[i,j] = binom.pmf(j+2, i+2, 0.5)
R = np.array([binom.pmf(1, i, 0.5) for i in range(2,M+1)])
B = np.linalg.solve(np.eye(M-1) - Q, R)
print('Using the fundamental matrix of the Markov chain gives a theoretical value of {}, for an initial number of coins of N0={}'.format(B[-1],M))

10000 simulations with N0=1000000 yield an estimated probability of 0.7181
Using the fundamental matrix of the Markov chain gives a theoretical value of 0.7213410794890949, for an initial number of coins of N0=5000

In [ ]: