How can I model flips until N successes?

Question

You and I decide to play a game where we take turns flipping a coin. The first player to flip 10 heads in total wins the game. Naturally, there is an argument about who should go first.

Simulations of this game show that the player to flips first wins 6% more than the player who flips second (the first player wins approx 53% of the time). I'm interested in modelling this analytically.

This isn't a binomial random variable, as there are no fixed number of trials (flip until someone gets 10 heads). How can I model this? Is it the negative binomial distribution?

---

So as to be able to recreate my results, here is my python code:

import numpy as np
from numba import jit


@jit
def sim(N):

    P1_wins = 0
    P2_wins = 0

    for i in range(N):

        P1_heads = 0
        P2_heads = 0
        while True:

            P1_heads += np.random.randint(0,2)

            if P1_heads == 10:
                P1_wins+=1
                break

            P2_heads+= np.random.randint(0,2)
            if P2_heads==10:
                P2_wins+=1
                break
    return P1_wins/N, P2_wins/N


a,b = sim(1000000)

Answer 1

The distribution of the number of tails before achieving $10$ heads is Negative Binomial with parameters $10$ and $1/2$. Let $f$ be the probability function and $G$ the survival function: for each $n\ge 0$, $f(n)$ is the player's chance of $n$ tails before $10$ heads and $G(n)$ is the player's chance of $n$ or more tails before $10$ heads.

Because the players roll independently, the chance the first player wins with rolling exactly $n$ tails is obtained by multiplying that chance by the chance the second player rolls $n$ or more tails, equal to $f(n)G(n)$.

Summing over all possible $n$ gives the first player's winning chances as

$$\sum_{n=0}^\infty f(n)G(n) \approx 53.290977425133892\ldots\%.$$

That is about $3\%$ more than half the time.

In general, replacing $10$ by any positive integer $m$, the answer can be given in terms of a Hypergeometric function: it is equal to

$$1/2 + 2^{-2m-1} {_2F_1}(m,m,1,1/4).$$

When using a biased coin with a chance $p$ of heads, this generalizes to

$$\frac{1}{2} + \frac{1}{2}(p^{2m}) {_2F_1}(m, m, 1, (1 - p)^2).$$

---

Here is an R simulation of a million such games. It reports an estimate of $0.5325$. A binomial hypothesis test to compare it to the theoretical result has a Z-score of $-0.843$, which is an insignificant difference.

n.sim <- 1e6
set.seed(17)
xy <- matrix(rnbinom(2*n.sim, 10, 1/2), nrow=2)
p <- mean(xy[1,] <= xy[2,])
cat("Estimate:", signif(p, 4), 
    "Z-score:", signif((p - 0.532909774) / sqrt(p*(1-p)) * sqrt(n.sim), 3))

Answer 2

We can model the game like this:

• Player A flips a coin repeatedly, getting results $A_1, A_2, \dots$ until they get a total of 10 heads. Let the time index of the 10th heads be the random variable $X$.
• Player B does the same. Let the time index of the 10th heads be the random variable $Y$, which is an iid copy of $X$.
• If $X \le Y$, Player A wins; otherwise Player B wins. That is, \begin{align} \Pr(A\text{ wins})&= \Pr(X \ge Y) = \Pr(X > Y) + \Pr(X = Y)\\ \Pr(B\text{ wins})&= \Pr(Y > X) = \Pr(X > Y). \end{align}

The gap in the win rates is thus $$ \Pr(X = Y) = \sum_k \Pr(X = k, Y = k) = \sum_k \Pr(X = k)^2 .$$

As you suspected, $X$ (and $Y$) are distributed essentially according to a negative binomial distribution. Notations for this vary, but in Wikipedia's parameterization, we have heads as a "failure" and tails as a "success"; we need $r = 10$ "failures" (heads) before the experiment is stopped, and success probability $p = \tfrac12$. Then the number of "successes," which is $X - 10$, has $$\Pr(X - 10 = k) = \binom{k + 9}{k} 2^{-10 - k},$$ and the collision probability is $$ \Pr(X = Y) = \sum_{k=0}^\infty \binom{k + 9}{k}^2 2^{-2k - 20} ,$$ which Mathematica helpfully tells us is $\frac{76\,499\,525}{1\,162\,261\,467} \approx 6.6\%$.

Thus Player B's win rate is $\Pr(Y > X) \approx 46.7\%$, and Player A's is $\frac{619\,380\,496}{1\,162\,261\,467} \approx 53.3\%$.

Answer 3

Let $E_{ij}$ be the event that the player on roll flips i heads before the other player flips j heads, and let $X$ be the first two flips having sample space $ \{ hh,ht,th,tt\}$ where h means heads and t tails, and let $p_{ij} \equiv Pr(E_{ij})$.

Then $p_{ij}=Pr(E_{i-1j-1}|X=hh)*Pr(X=hh)+Pr(E_{i-1j}|X=ht)*Pr(X=ht)+Pr(E_{ij-1}|X=th)*Pr(X=th)+Pr(E_{ij}|X=tt)*Pr(X=tt)$

Assuming a standard coin $Pr(X=*)=1/4$ means that $p_{ij}=1/4*[p_{i-1j-1}+p_{i-1j}+p_{ij-1}+p_{ij}]$

solving for $p_{ij}$, $= 1/3*[p_{i-1j-1}+p_{i-1j}+p_{ij-1}]$

But $p_{0j}=p_{00}=1$ and $p_{i0}=0$, implying that the recursion fully terminates. However, a direct naive recursive implementation will yield poor performance because the branches intersect.

An efficient implementation will have complexity $O(i*j)$ and memory complexity $O(min(i,j))$. Here's a simple fold implemented in Haskell:

Prelude> let p i j = last. head. drop j $ iterate ((1:).(f 1)) start where
  start = 1 : replicate i 0;
  f c v = case v of (a:[]) -> [];
                    (a:b:rest) -> sum : f sum (b:rest) where
                     sum = (a+b+c)/3 
Prelude> p 0 0
1.0
Prelude> p 1 0
0.0
Prelude> p 10 10
0.5329097742513388
Prelude> 

UPDATE: Someone in the comments above asked whether one was suppose to roll 10 heads in a row or not. So let $E_{kl}$ be the event that the player on roll flips i heads in a row before the other player flips i heads in a row, given that they already flipped k and l consecutive heads respectively.

Proceeding as before above, but this time conditioning on the first flip only, $p_{k,l} = 1-1/2*[p_{l,k+1}+p_{l,0}]$ where $p_{il}=p_{ii}=1, p_{ki}=0$

This is a linear system with $i^2$ unknowns and one unique solution.

To convert it into an iterative scheme, simply add an iterate number $n$ and a sensitivity factor $\epsilon$:

$p_{k,l,n+1} = 1/(1+\epsilon)*[\epsilon*p_{k,l,n} +1-1/2*(p_{l,k+1,n}+p_{l,0,n})]$

Choose $\epsilon$ and $p_{k,l,0}$ wisely and run the iteration for a few steps and monitor the correction term.