Subtracting random number does not match average I should get

Question

Suppose you organize a competition with 25 million participants. At each round, a random number from the people remaining is eliminated. How many rounds should we expect to get 5 participants or less ? Now 25M/2^22 is the closest number higher than 5 so 22 rounds on average. However, when verifying this in Python, I get that the distribution is more accumulated around 15-16 rather than 22-23. (I am not very familiar with Python so I know I should use dictionaries to count stuff etc.) Is there a mistake in my code ? Because theory cannot lie.

import random

lst=list()
for j in range(1,10000):
    i=0
    X=25000000
    while X > 5:
        i = i+1
        elim=(random.randint(0,X))
        X = X - elim
    lst.append(i)
    
for i in range(1,33):
    print( i, "appears", lst.count(i), "times")

print("22 appears", lst.count(22)*100/len(lst), "% of the time")
print("15 appears", lst.count(15)*100/len(lst), "% of the time")

Answer 1

Expect about 16.33 rounds.

Let e[X] be the expected number of rounds for starting with X participants.

X <= 5:
  e[X] = 0

X > 5:
  e[X] = 1 + (e[X]/(X+1) + ... + e[0]/(X+1))
  e[X]*(X+1) = (X+1) + (e[X] + ... + e[0])
  e[X]*X = (X+1) + (e[X-1] + ... + e[0]
  e[X] = 1 + 1/X + (e[X-1] + ... + e[0])/X

In Python:

total = 0
for X in range(6, 25_000_001):
    e = 1 + 1/X + total/X
    total += e
print(e)

Output:

16.32826873439911

Running your experiment three times, printing the average sum(lst) / len(lst):

16.352835283528353
16.27162716271627
16.315431543154315