Which Geyser Gushes First?

You arrive at Three Geysers National Park, where geysers $A$, $B$ and $C$ erupt at intervals of precisely two, four and six hours respectively. You have no idea how the eruptions are staggered: assume each geyser started at some independently random point in history. What are the probabilities that $A$, $B$ and $C$, respectively, will be the first to erupt after your arrival?

The Riddler, FiveThirtyEight(original post)

Solution

Because each geyser’s starting time is random and it erupts on a fixed cycle, the time you must wait for its next eruption is uniform on one full period. So the waits are independent with $$a \sim \mathcal U(0,2), \qquad b \sim \mathcal U(0,4), \qquad c \sim \mathcal U(0,6),$$ and geyser $A$ wins exactly when $a$ is the smallest of the three. Conditioning on $a$, the chance $b$ exceeds it is $(4-a)/4$ and the chance $c$ exceeds it is $(6-a)/6$, so $$\begin{aligned} \Pr(A\text{ first}) &= \int_0^2 \frac12 \cdot \frac{4-a}{4}\cdot\frac{6-a}{6}\,da\\ &= \frac{1}{48}\int_0^2 (4-a)(6-a)\,da = \frac{1}{48}\cdot\frac{92}{3} = \frac{23}{36}. \end{aligned}$$ For $B$ to win it must beat $a$, but $a \le 2$ always, so any $b > 2$ is already too late; only $b \in (0,2)$ can win, giving $$\Pr(B\text{ first}) = \int_0^2 \frac14 \cdot \frac{2-b}{2}\cdot\frac{6-b}{6}\,db = \frac{1}{48}\cdot\frac{32}{3} = \frac{8}{36}.$$ The three probabilities sum to $1$, so $C$ takes what is left: $$\Pr(A), \Pr(B), \Pr(C) = \boxed{\tfrac{23}{36},\ \tfrac{8}{36},\ \tfrac{5}{36}}.$$ The shortest cycle wins most often, as it should: $A$ erupts first nearly two-thirds of the time.

The computation

Replay the visit. Draw each geyser’s wait uniformly over its own period, see which is smallest, and tally the winners over many arrivals.

import numpy as np
rng = np.random.default_rng(0)
N = 10_000_000

a = rng.uniform(0, 2, N)      # wait until A's next eruption
b = rng.uniform(0, 4, N)
c = rng.uniform(0, 6, N)

print("A:", round((( a < b) & (a < c)).mean(), 4))
print("B:", round((( b < a) & (b < c)).mean(), 4))
print("C:", round((( c < a) & (c < b)).mean(), 4))
# A: 0.6391    (23/36 = 0.6389)
# B: 0.2222    ( 8/36 = 0.2222)
# C: 0.1387    ( 5/36 = 0.1389)