How Long Will The Bacterial Colony Last?

Riddler Express

In “bowling” a die you pinch two opposite faces, so the die lands on one of the other four faces, each with probability $\tfrac14$. You want to roll a $7$ or $11$ with two dice (standard rolling wins $2/9 \approx 22.2\%$). Bowling the dice one at a time, knowing the first result before the second, what is your best chance of a $7$ or $11$? Extra credit: now score $+1$ for a $7$ or $11$ and $-1$ for a $2$, $3$ or $12$; bowling to maximise the expected score, what is that average?

The Riddler, FiveThirtyEight, June 12, 2020(original post)

Solution

Opposite faces of a die sum to $7$, so a bowl pinches one of the pairs $\{1,6\}$, $\{2,5\}$ or $\{3,4\}$ and leaves the other four faces equally likely. Bowl the first die to keep your options open: pinching $\{3,4\}$ leaves $\{1,2,5,6\}$, so half the time you roll a $5$ or $6$. That matters because a first roll of $5$ can be completed to $7$ (with a $2$) or to $11$ (with a $6$), and a $6$ to $7$ or $11$ likewise, doubling the ways to win, whereas a first roll of $1,2,3,4$ can only reach $7$.

Given the first result, bowl the second die to keep the faces you need. Working through the four equally likely first outcomes and choosing the best second pinch each time, the chance of finishing on $7$ or $11$ is $$\boxed{\tfrac38 = 37.5\%},$$ a big jump from $22.2\%$. For the scoring version, the penalty faces $2,3,12$ change which first pinch is best (now $\{1,6\}$), and the optimal expected score works out to $$\boxed{\tfrac{5}{16} = 0.3125},$$ about three times the $\tfrac19$ of ordinary rolling.

The computation

Encode the bowl as a choice of which opposite pair to pinch (leaving four equiprobable faces). Optimise the second pinch for each first outcome, then the first pinch, for both the win-probability objective and the scoring objective.

from fractions import Fraction as F
pairs = [(1, 6), (2, 5), (3, 4)]
strats = [[f for f in range(1, 7) if f not in p] for p in pairs]
def best(score):
    return max(sum(max(sum(score(d1 + d2) for d2 in s2) / 4 for s2 in strats)
                   for d1 in s1) / 4 for s1 in strats)
win = best(lambda t: 1 if t in (7, 11) else 0)
ec  = best(lambda t: 1 if t in (7, 11) else (-1 if t in (2, 3, 12) else 0))
print(F(win).limit_denominator(), F(ec).limit_denominator())   # 3/8 5/16

Riddler Classic

Each bacterium of Riddlerium classicum either splits into two copies (probability $80\%$) or dies (probability $20\%$). Starting from a single bacterium, what is the probability the colony lasts forever? Extra credit: with split probability $p$ in general?

The Riddler, FiveThirtyEight, June 12, 2020(original post)

Solution

This is a branching process. Let $q$ be the probability that the colony dies out, starting from one bacterium. That first bacterium either dies at once (probability $1-p$) or splits into two independent colonies, each of which must itself die out (probability $q^2$). So $$q = (1-p) + p\,q^2.$$ This quadratic has roots $q = 1$ and $q = \dfrac{1-p}{p}$. Extinction takes the smaller root. For $p \le \tfrac12$ that is $q = 1$: extinction is certain. For $p > \tfrac12$ it is $q = \tfrac{1-p}{p}$, so the colony survives with probability $1 - q = \tfrac{2p-1}{p}$.

At $p = 0.8$, $q = \tfrac{0.2}{0.8} = \tfrac14$, so the everlasting probability is $$1 - \tfrac14 = \boxed{\tfrac34}.$$ In general, $$\Pr(\text{everlasting}) = \boxed{\max\!\left(0,\ \tfrac{2p-1}{p}\right)}.$$

The computation

Encode the process directly: simulate colonies generation by generation and measure how often one survives past a deep cutoff. Compare with the closed form $\tfrac{2p-1}{p}$.

import random
def survives(p, gens=200, trials=200000, seed=0):
    rng = random.Random(seed); lived = 0
    for _ in range(trials):
        pop = 1
        for _ in range(gens):
            pop = sum(2 if rng.random() < p else 0 for _ in range(pop))
            if pop == 0: break
            if pop > 5000: break           # treat as everlasting (runaway)
        lived += pop > 0
    return lived / trials
for p in (0.8, 0.6, 0.5):
    print(p, round(survives(p), 3), round(max(0, (2 * p - 1) / p), 3))
# 0.8 ~0.75 0.75 | 0.6 ~0.333 0.333 | 0.5 ~0.01 0.0 (critical: ~0 over a finite horizon)