Expected loops in shoelaces

In this note, we solve a problem from Introduction to Probability (Blitzstein, Hwang)

Question

There are 100 shoelaces in a box. At each stage, you pick two random ends and tie them together. Either this results in a longer shoelace (if the two ends came from different pieces), or it results in a loop (if the two ends came from the same piece). What are the expected number of steps until everything is in loops, and the expected number of loops after everything is in loops?

Solution

Start with small case, let $f(n)$ be the number of loops that can be formed. $\mathbb{E}[f(1)] = 1$, since with probability 1 the two ends of the only shoelace are chosen. Then,

$$ \mathbb{E}\left[ f(2) \right] = 2/{4\choose 2}\cdot (1+f(1)) + \left( 1-2/{4\choose 2} \right)\cdot f(1) = f(1) + 1/3 = 1+1/3 $$

Pattern:

$$ \mathbb{E}[f(n)] = 1+1/3+1/5+\cdots +1/(2n-1) $$

Induction:

$$ \begin{align*} \mathbb{E}[f(n)] &= \frac{n}{{n\choose 2}}\cdot (1+f(n-1)) + \left( 1-\frac{n}{{n \choose 2}} \right)f(n-1) \\ &= f(n-1) + \frac{n}{{2n\choose 2}} \\ &= f(n-1) + \frac{1}{2n-1} \end{align*} $$

We provide a few codes to simulate this result.

# analytic
def expected(n):
    res = 0
    for i in range(1, n+1):
        res += (1 / (2 * i - 1))
    return res 

#print(expected(100))  # 3.28434...

# simulation
import numpy as np
def simulate_outer(n):
    memo = {} 
    def simulate(n):
        """ generates a stochastic ooutput in n noodles"""
        if n == 1: return 1
        if n in memo: return memo[n]

        # generate a probability 
        a = np.random.uniform() 
        if a <= 1/(2*n-1):
            return 1+simulate(n-1)
        else:
            return simulate(n-1)
    return simulate(n)

# test
#n = 100
#nmc = int(1e4)
#print(sum([simulate_outer(n) for _ in range(nmc)]) / nmc) # 3.2872

# can also simulate noodles themselves
def simulate_noodles(n):
    """simulate the number of loops formed by randomly connecting n noodles"""
    # ends -> noodle
    ends = {i: None for i in range(2*n)}
    # assign noodles
    for i in range(n):
        # (the other end, which noodle)
        ends[2*i] = (2*i+1, i)
        ends[2*i+1] = (2*i, i)
    loops = 0
    k = n
    while k > 0:
        left, right = np.random.choice(
            list(ends.keys()), 
            2, replace=False
        )
        # if belong to the same noodle
        if ends[left][1] == ends[right][1]:
            loops += 1
            ends.pop(left)
            ends.pop(right)
        # if does not belong to same noodle, combine
        else:
            other_end1, noodle1 = ends.pop(left)
            other_end2, noodle2 = ends.pop(right) 
            ends[other_end1] = (other_end2, noodle1)
        k -= 1
    return loops