Lindy effect in surnames problem

I have a challenge for you: Assume you have a population of 100 000, with distinct surnames, and the same number of males and females. They randomly pair with each other, and make next population. Each pair has random number of children, from 0 to 5. Child takes surname from one of the parents. The gender of each child is random. Then population is limited to 100 000. How many distinct surnames there will be after 100 000 iterations? What is the relationship between the time how long a surname has survived, and what is its expected survival time?

If you are good at mathematics, you can find answers with complex formulas. However, I suggest simulating this situation, and see the results.

My solution

I used Kotlin, to find out the result. With every generation, every surname has a change of being eliminated from the pool, so it seems intuitive, that over time there will be only one surname left. The question is, how long it is going to take. I will model a person as a class with a surname (I used Int, because why not) and a biological gender (for this exercise, I assumed Boolean). To make a population of 100_000 people, I used List builder, and made a random people, each with a different surname. Then for each iteration, I shuffle my population, limit it, and randomly assign males to females. Then each pair have a random number of children, from 0 to 5, what makes the new population. Then I print the iteration number and the number of distinct surnames.

import kotlin.random.Random

class Person(val surname: Int, val male: Boolean)

fun newPerson(surname: Int) = Person(
    surname = surname,
    male = Random.nextBoolean()
)

fun main() {
    var population = List(100_000) { newPerson(it) }
    for (i in 1..100_000) {
        val (males, females) = population.shuffled()
            .take(100_000)
            .partition { it.male }
        val pairs = males.zip(females)
        population = pairs.flatMap { (m, _) ->
            val childrenNum = Random.nextInt(6)
            List(childrenNum) { newPerson(m.surname) }
        }
        println("Iteration: $i")
        println("Surnames num: ${population.distinctBy { it.surname }.size}")
    }
}

It mu case, it took around 6000 iterations, until only two surnames were left. Fights for full domintation took around 2500 more iterations. It seems nearly impossible, to have more than one surname after 100 000 iterations. It seems quite clear, that the longer a surname survived, the less surnames it has to fight with, the longer expected survival time it has. This perfecty demonstrates the Lindy effect. Let's see it in our simutation.

We will first see, what is the distribution of the number of surnames, for that, we will note how many

import kotlin.random.Random

class Person(val surname: Int, val male: Boolean)

fun newPerson(surname: Int) = Person(
    surname = surname,
    male = Random.nextBoolean()
)

fun main() {
    var population = List(100_000) { newPerson(it) }
    var surnamesLost = mapOf<Int, Int>()
    var prevSurnamesCount = population.size

    for (i in 1..1_000_000) {
        val (males, females) = population.shuffled()
            .take(100_000)
            .partition { it.male }
        val pairs = males.zip(females)
        population = pairs.flatMap { (m, _) ->
            val childrenNum = Random.nextInt(6)
            List(childrenNum) { newPerson(m.surname) }
        }
        val surnamesCount = population.distinctBy { it.surname }.size
        val surnamesLostCount = prevSurnamesCount - surnamesCount
        if (surnamesLostCount > 0) {
            surnamesLost += i to surnamesLostCount
        }
        prevSurnamesCount = surnamesCount
        if (i % 1000 == 0) {
            println(i)
            println(surnamesCount)
            println(surnamesLost)
        }
        if (surnamesCount <= 1) {
            break
        }
    }
    println(surnamesLost) // {1=58396, 2=14483, 3=6725, 4=3858, 5=2636, 6=1904, 7=1430,
    // 8=1127, 9=872, 10=743, 11=660, 12=543, 13=451, 14=403, 15=351, 16=304, 17=286,
    // 18=241, 19=197, 20=225, 21=182, 22=163, 23=155, 24=152, 25=126, 26=107, 27=110,
    // 28=106, 29=91, 30=101, 31=79, 32=86, 33=90, 34=62, 35=65, 36=79, 37=54, 38=66,
    // 39=52, 40=49, 41=45, 42=37, 43=57, 44=55, 45=37, 46=31, 47=50, 48=39, 49=36,
    // 50=41, 51=35, 52=27, 53=43, 54=31, 55=30, 56=19, 57=26, 58=29, 59=17, 60=22,
    // 61=20, 62=19, 63=23, 64=22, 65=20, 66=22, 67=18, 68=22, 69=18, 70=22, 71=14,
    // 72=19, 73=22, 74=19, 75=8, 76=15, 77=16, 78=21, 79=16, 80=10, 81=22, 82=15,
    // 83=16, 84=10, 85=13, 86=16, 87=9, 88=9, 89=8, 90=11, 91=10, 92=17, 93=9, 94=8,
    // 95=8, 96=6, 97=4, 98=8, 99=6, 100=2, 101=10, 102=6, 103=10, 104=6, 105=11, 106=6,
    // 107=6, 108=4, 109=2, 110=4, 111=6, 112=11, 113=9, 114=17, 115=7, 116=8, 117=5,
    // ...
    // 3421=1, 3431=1, 3802=1, 4168=1, 4807=1, 4870=1, 5375=1, 5598=1, 6755=1, 6887=1,
    // 6919=1, 8170=1, 8180=1, 8460=1, 8788=1}
}

I started that in REPL, to be able to operate on the result surnamesLost variable, without recalculating it every time.

Notice, that the algorithm does not include the first iteration, when over half of surnames are lost. I decided to exclude it, because it is not a natural situation, and it will affect other estimations too much.

What is the average survival rate after the first round? This can be calculated by waighted average of iterations a surname survived. We will exclude the last surname, that will live forever.

println(
    surnamesLost.toList().sumOf { (i, num) -> i * num }.toFloat()
    / surnamesLost.toList().sumOf { (i, num) -> num }
) // 10.163172

How does it change, if we ignore first 10 rounds?

val survivalRateAfter10rounds = surnamesLost
    .filter { (i) -> i > 10 }
    .mapKeys { (i) -> i - 10 }
    .toList()
println(
    survivalRateAfter10rounds.sumOf { (i, num) -> i * num }.toFloat()
    / survivalRateAfter10rounds.sumOf { (i, num) -> num }
) // 96.74638

Only 10 round has passed, and now the expected survival time is nearly 100 iterations. What is it after 100 iterations?

val survivalRateAfter100rounds = surnamesLost
    .filter { (i) -> i > 100 }
    .mapKeys { (i) -> i - 100 }
    .toList()
println(
    survivalRateAfter100rounds.sumOf { (i, num) -> i * num }.toFloat()
    / survivalRateAfter100rounds.sumOf { (i, num) -> num }
) // 627.6806

As you can see, the longer a surname survived, the longer expected survival rate it has. Lindy effect. Let's see its distribution for 1000 iterations.

var expectedSurvivalRates = listOf<Double>()
for (iterations in 0..1000) {
    val survivalRate = surnamesLost
        .filter { (i) -> i > iterations }
        .mapKeys { (i) -> i - iterations }
        .toList()
    expectedSurvivalRates += survivalRate.sumOf { (i, num) -> i * num }.toDouble() / survivalRate.sumOf { (i, num) -> num }
}
println(expectedSurvivalRates) // [10.163171631716317, 21.93593794886527, 32.097775820303895, 41.19257134522826, 50.16815340736695, 58.42950890603591, 66.8324102047667, 74.76022109978081, 82.31929376728355, 90.0440466376163, 96.74638515674984, 103.72303853617964, 110.62009255112703, 117.64578660685677, 124.41951469583049, 131.15144361721445, 137.65008642212408, 143.71581498687135, 151.07675860594398, 158.12322595179094, 165.47330960854092, 173.3570892723181, 179.54050533993228, 185.9010035259018, 193.67329545454547, 200.1209796400118, 206.4935740514076, ..., 3714.4891304347825, 3713.4891304347825]

Estimated expected survival time, based on the numbers of survived rounds so far.

Clearly, together with surviving this iteration, the chance to survive another iteration is growing. How much? We can calculate and visualize it.

val expectedSurvivalRatesDiff = expectedSurvivalRates
    .zipWithNext { prev, next -> next - prev }

The growth of estimated expected survival time, based on the number of survived rounds so far. The graph is not beautiful, because it is based on a single simulation.

...