I was curious a few weeks ago about how surnames come and go (mostly on the go). I wanted to know how quickly surnames die out, because (in my simplified simulation world) once the last person with a surname dies, it's gone forever. I started writing a small simulation to see this happen in "real time". I started with a simple "person" representation that behaved exactly like all the other "persons" in the simulation. However, as I got into the project, I started coming up with certain traits, since not everyone lives to the same age, not everyone has the same number of children, not everyone pairs up to even have children, etc.

A few assumptions:

  • I didn't base this simulation on any actual population statistics, though I think it does a decent job at producing some interesting results.
  • Surnames do not regenerate. The starting population each gets a unique "name" (an integer one greater than the last), and that forms the maximum amount of names in the simulation.
  • This is paternally biased. To keep things simple, females of the population take the male surname when paired. Surnames are never combined (see point 2), or modified in any way. This simulation would work the same (mostly?) if maternally biased, though I suppose at that point males are only needed for population growth (topic not relevant for this simulation, just some side thoughts).
  • Surnames are inherited from the parent. I didn't originally consider this an assumption, but according to Toby Speight, this is not always true in every culture, with Iceland being the example.
  • People die when they reach their maximum lifespan, but can also die early due to "accidents", being more prone the old they are.
  • When a person with a partner dies, the surviving partner does not gain a new partner.
  • Exponential distributions are used here to help move away from uniform distributions, since I found that uniform populations die really easily. An example would be that a Person is more likely to have a life span of 80 than 20, but that drops off significantly past 80.

Therefore, I came up with the following Person object to act as the individual units of my simulation:

# Person.py

import random
import numpy as np

male = 1
female = 2

# these reduce performance contributions from lookups by half, and because of how often they are called, any bit helps
# # $ python -m timeit -s "import numpy as np" "np.random.randint(0, 1000)"
# # 10000000 loops, best of 3: 0.168 usec per loop
# # $ python -m timeit -s "import numpy as np; npRandInt = np.random.randint" "npRandInt(0, 1000)"
# # 10000000 loops, best of 3: 0.0814 usec per loop
npRandInt = np.random.randint
npRandUni = np.random.uniform
npRandEx = np.random.exponential

class Person(object):
    currentName = 0
    deaths = 0
    maxAgeDistribution = []
    maxChildDistribution = [0, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5]

    def __init__(self, parent=None):
        if parent is None:
            Person.currentName += 1
            self.surname = Person.currentName
            self.surname = parent.surname

        self.gender = random.choice([male, female])

        self.currentChildCount = 0
        self.currentAge = 0

        self.maxChildren = random.choice(Person.maxChildDistribution)  # cap how many children one person can have
        self.maxAge = random.choice(Person.maxAgeDistribution)

        self.childDifficulty = npRandInt(0, 100 + 1) / 100.0  # how difficult is it to have a child
        self.childYearsStart = npRandInt(10, 15 + 1)
        self.childYearsEnd = npRandInt(30, 50 + 1)

        self.partner = None  # some other person object
        self.withChild = False
        self.hasHadPartner = False  # give a way to check if they have had a partner before, so if one dies, they don't get another

        # determine age person is eligible to gain a partner. Having this too high (on average) leads to excessive population death
        if self.childYearsStart < self.maxAge:
            self.partnerEligibleAge = self.childYearsStart + int(self.childYearsStart * npRandEx(0.2))

            if self.partnerEligibleAge > self.maxAge:
                self.partnerEligibleAge = self.maxAge
            self.partnerEligibleAge = self.maxAge

    def __repr__(self):
        return 'Person(Gender: %6s, Surname: %5d, Age: %2d, MaxAge: %2d, Children: %1d, MaxChildren: %1d, Partner: %5s, PartnerEligibleAge: %2d)' % \
               ("Male" if self.isMale() else "Female", self.surname, self.currentAge, self.maxAge, self.currentChildCount, self.maxChildren,
                "True" if self.partner is not None else "False", self.partnerEligibleAge)

    def gainPartner(self, partner):  # gain a partner if both persons are eligible, female taking the surname of the other partner
        if self.isPartnerEligible():
            if self.gender == female:
                self.surname = partner.surname

            self.partner = partner
            self.hasHadPartner = True

    def age(self):
        self.currentAge += 1

        if self.currentAge >= self.maxAge:  # if they outlive their lifespan, they die
            return self.die()

        percentAge = 1.0 - ((self.maxAge - self.currentAge) / float(self.maxAge))  # get their current age as a percent of their max age

        # get a random [0, 1] with an exponential distribution.
        # Magic numbers were determined experimentally, seeing what allowed population growth, and then by checking the distribution using pyplot
        randDeathPercent = 1.0 - npRandEx(0.05) * 2.0

        # random accidents happen, so even if they haven't reached their max age,
        # there is some risk of death, especially as they get closer to their max age
        return self.die() if percentAge > randDeathPercent else False  # perform self.die, which returns True, if they have died

    def gainChild(self):  # gain a child if eligible and passes a "difficulty" check
        if self.isChildEligible():
            chance = 1.0 - (self.childDifficulty * self.partner.childDifficulty)
            if chance > npRandUni(0, 1):
                self.withChild = True

    def haveChild(self, modifier=2.0):  # produce child if they currently have one and pass a modified "difficulty" check
        if self.withChild is True:
            self.withChild = False  # pass or fail, they will no longer have a child

            chance = 1.0 - self.childDifficulty
            if chance * modifier > npRandUni(0, 1):  # if this modifier is set too low, not enough children are created and the population dies
                self.currentChildCount += 1
                return Person(parent=self)  # return a new person, giving it the last name of the parent

        return None  # if we get this far, a child fails to be created

    def die(self):
        # if they have a partner, remove themselves from that person
        if self.partner is not None:
            self.partner.partner = None
        Person.deaths += 1  # add to the death stats
        return True

    def hasPartner(self):
        return self.partner is not None

    def isChildEligible(self):  # check if person can have a child
        # Criteria:
            # only females can have children
            # they must have a partner
            # they haven't had more than their max children
            # they are within child bearing years, and so is their partner
            # they don't currently have a child
        return self.gender == female and \
               self.partner is not None and \
               self.currentChildCount < self.maxChildren and \
               self.childYearsStart <= self.currentAge <= self.childYearsEnd and \
               self.partner.childYearsStart <= self.partner.currentAge <= self.partner.childYearsEnd and \
               self.withChild is False

    def isPartnerEligible(self):  # check if person can gain a partner
        # Criteria:
            # they don't currently have a partner
            # they are old enough to have a partner
            # they haven't had a partner in the past
        return self.partner is None and \
               self.currentAge >= self.partnerEligibleAge and \
               self.hasHadPartner is False

    # helper functions to check for gender.
    # Performance profiling indicates these calls are expensive over time, so it is better to check this directly
    def isMale(self):
        return self.gender == male

    def isFemale(self):
        return self.gender == female

Now I needed to actually run the simulation. After a lot of performance checking, both profiling and seeing the results of the simulation, I was able to come up with the following main application, with main handling the bulk of the administrative work, and sim doing the actual simulation:

# main.py

import random
import numpy as np
import cProfile
import operator
import datetime
import math
import matplotlib.pyplot as plt
import cPickle as pickle
import person
from person import Person
from stats import Stats
import sys
import uuid
import time

npRandInt = np.random.randint
npRandUni = np.random.uniform
npRandEx = np.random.exponential

fileId = str(uuid.uuid1())

def getTimeDelta(start, stop):
    delta = stop - start
    return (delta.microseconds + (delta.seconds + delta.days * 24 * 3600) * 10 ** 6) / (10 ** 6 * 1.0) * 1000.0

def setupAgeDistribution(distribution):
    for i in range(10000):

    for i in range(2000):

# get a random age in range [0, 80]
# Magic numbers were determined experimentally,
# seeing what allowed population growth,
# and then by checking the distribution using pyplot
def getRandomAgeTo80():
    ran = 80.0 - npRandEx(50.0) * .4
    while ran < 0:
        ran = 80.0 - npRandEx(50.0) * .4
    return ran

# get a random age in range [80, 100]
# Magic numbers were determined experimentally,
# seeing what allowed population growth,
# and then by checking the distribution using pyplot
def getRandomAgeFrom80():
    ran = 80 + npRandEx(50.0) * .1
    while ran > 100.0:
        ran = 80 + npRandEx(50.0) * .1
    return ran

def countItems(aList):  # no access to collections.Counter, so use this instead
    array = np.bincount(np.array(aList))
    nonZero = np.nonzero(array)[0]
    return dict(zip(nonZero, array[nonZero]))

def gainChildren(people):  # have each person try to gain a child, letting the eligibility be determined by the object
    for p in people:

def gainPartners(people):
    # get a list of eligible partners, split by gender
    mPartners = [x for x in people if x.gender == person.male and x.isPartnerEligible()]
    fPartners = [x for x in people if x.gender == person.female and x.isPartnerEligible()]

    numPartners = min([len(mPartners), len(fPartners)])  # get the smaller count of partners

    if numPartners > 0:
        numChoices = npRandInt(0, numPartners + 1)  # select a number of them to gain a partner

        if numChoices > 0:
            # shuffle both lists, splice to make lists the same length, zip them together, and pair up the zipped pairs.
            # Much faster than using random.choice.

            pairs = zip(mPartners[:numChoices], fPartners[:numChoices])

            for m, f in pairs:

def haveChildren(people):
    newCount = 0
    for p in people.copy():
        if p.withChild is True:
            newChild = p.haveChild(modifier=1.0)
            if newChild is not None:
                newCount += 1
    return newCount

def ageSim(people):
    for p in people.copy():
        if p.age():  # if they ded

def sim(people):
    ageSim(people)  # age the whole sim, removing those that die
    gainPartners(people)  # pair up eligible people objects
    newCount = haveChildren(people)  # have children before gaining them to "develop" them for a year
    gainChildren(people)  # gain new children to attempt being created next year

    return newCount  # return how many total new people have been added

def main(people=None, stats=None):  # arguments allow resuming from serialized (saved) data

    if people is None and stats is None:
        # set up the stats object, setting the initial population count, and how often to print the stats
        stats = Stats(oldDate=datetime.datetime.now(), startingPopulation=10000, statYear=100.0)
        people = set([Person() for i in range(stats.startingPopulation)])  # create a new set of random people

    stats.oldDate = datetime.datetime.now()

        for year in range(stats.year, stats.year + 500 + 1):
            stats.year = year

            if len(people) == 0:
                print "Population has died! Final stats below:\n"
                printStats(people, stats)

            if year % stats.statYear == 0:
                printStats(people, stats)
            stats.newPersonCount += sim(people)

    except KeyboardInterrupt:  # capture any ctrl-c, and print stats on exit
        printStats(people, stats)

        if len(people) > 0:
            sortedList = sorted(countItems([p.surname for p in people]).items(), key=operator.itemgetter(1))  # sort list by surname popularity
            if len(sortedList) < 1000:  # only print it if the sorted list isn't too long. Console has an issue if the line is too long.
                sortedList.reverse()  # display the list in descending order
                print sortedList


    # plot(people, stats)

def printStats(people, stats):
    statYear = stats.statYear
    if stats.year % statYear != 0:
        statYear = stats.year % stats.statYear

    delta = getTimeDelta(stats.oldDate, datetime.datetime.now())
    timeTaken = math.floor(delta / 1000.0)

    if timeTaken < 1:
        time.sleep(1)  # don't spam the console output

    population = len(people)

    meanAge = math.floor(np.mean([person.currentAge for person in people])) if len(people) > 0 else 0

    peopleDelta = stats.newPersonCount - stats.oldPeople
    stats.oldPeople = stats.newPersonCount

    surnameCount = len(countItems([p.surname for p in people])) if len(people) > 0 else 0

    deathsDelta = Person.deaths - stats.oldDeaths
    stats.oldDeaths = Person.deaths

    netPop = (peopleDelta / statYear) - (deathsDelta / statYear)

    pairCount = len([True for person in people if person.partner is not None]) if len(people) > 0 else 0

    print datetime.datetime.now()
    print "Time taken: %d seconds (%d milliseconds / year)" % (timeTaken, math.floor(delta / statYear))
    print "Current population size at year %d: %d" % (stats.year, population)
    print "Current mean age: %d" % meanAge
    print "Current new people: %d" % stats.newPersonCount
    print "Current deaths: %d" % Person.deaths
    print "Current delta new people: %d (%d people / year)" % (peopleDelta, (peopleDelta / statYear))
    print "Current delta deaths: %d (%d deaths / year)" % (deathsDelta, (deathsDelta / statYear))
    print "Current surname count: %d" % surnameCount
    print "Current partner percentage: %d" % (pairCount / float(population) * 100) if len(people) > 0 else 0
    print "Net population per year: %d" % netPop
    print ""
    stats.oldDate = datetime.datetime.now()

    if SERIALIZE:  # controls whether serialization is enabled, useful to turn off for when debugging or testing
        serialize(people, stats, fileId)

def serialize(people, stats, id=""):
    with open("people%s.dat" % id, "wb") as f:
        pickle.dump(people, f, pickle.HIGHEST_PROTOCOL)

    with open("stats%s.dat" % id, "wb") as f:
        pickle.dump(stats, f, pickle.HIGHEST_PROTOCOL)

def deserialize(people=True, stats=True, id=""):
    peopleData = None
    statsData = None

    if people:
        with open("people%s.dat" % id, "rb") as f:
            peopleData = pickle.load(f)

    if stats:
        with open("stats%s.dat" % id, "rb") as f:
            statsData = pickle.load(f)

    return peopleData, statsData

def plot(people, stats):  # experimental plotting function, currently need a rework using subplots
    numplots = len(stats.getStats())
    colormap = plt.cm.gist_ncar
    plt.gca().set_color_cycle([colormap(i) for i in np.linspace(0, 0.9, numplots)])

    labels = []
    for stat, label in stats.getStats():


    plt.legend(labels, ncol=2, loc='upper center', bbox_to_anchor=[0.5, 1.1], columnspacing=1.0,
               labelspacing=0.0, handletextpad=0.0, handlelength=1.5, fancybox=True, shadow=True)

    # plt.yscale('log')

    # Helpful diagnostic histogram for various stats

    # bins = np.arange(0, 100+1, 1)
    # plt.hist([person.maxAge for person in people], bins=bins)
    # plt.show()

if __name__ == '__main__':
    # deserializing usually happens here, but I have excluded that portion here

I noticed that I needed a way to keep track of stats over time, so I could see this as plots and not just printed stats. I wrote a helper class Stats to do this:

# stats.py

class Stats(object):
    def __init__(self, oldDate, startingPopulation, statYear):
        self.timeTaken = []
        self.population = []
        self.meanAge = []
        self.newPeople = []
        self.deaths = []
        self.deltaNewPeople = []
        self.surnames = []
        self.deltaDeaths = []
        self.netPopDelta = []
        self.startingPopulation = startingPopulation
        self.year = 0
        self.statYear = statYear
        self.oldDate = oldDate
        self.newPersonCount = 0
        self.oldDeaths = 0
        self.oldPeople = 0

    def getStats(self):  # convenience function to get all plot-able stats as an iterable with their labels attached
        return [
            (self.timeTaken, "Time Taken"),
            (self.population, "Pop Count"),
            (self.meanAge, "Mean Age"),
            (self.newPeople, "New People"),
            (self.deaths, "Deaths"),
            (self.deltaNewPeople, "Delta New People"),
            (self.deltaDeaths, "Delta Deaths"),
            (self.surnames, "Surnames"),
            (self.netPopDelta, "Net Pop Change")

After a lot of trial and error and diagnostic plots, I feel I finally got this working decently. By no means do I claim that the results are statistically accurate for any real world population, but for what I set out to explore, the results seem reasonably OK, especially since The Atlantic notes that for China, with 1.3 billion people, 87% of the population shares only 100 unique surnames.

For a country of 1.3 billion people, there is a remarkably small number of common last names in China. An estimated 87 percent of the population shares one of 100 surnames, and more than one in five Chinese citizens is surnamed Li, Wang, or Zhang: more than 275 million people in all.

A sample run was done, simulating 500 years, printing stats every 100 years:

2018-03-08 15:59:02.250965
Time taken: 0 seconds (0 milliseconds / year)
Current population size at year 0: 10000
Current mean age: 0
Current new people: 0
Current deaths: 0
Current delta new people: 0 (0 people / year)
Current delta deaths: 0 (0 deaths / year)
Current surname count: 10000
Current partner percentage: 0
Net population per year: 0

2018-03-08 15:59:11.135088
Time taken: 8 seconds (88 milliseconds / year)
Current population size at year 100: 26898
Current mean age: 25
Current new people: 48750
Current deaths: 31852
Current delta new people: 48750 (487 people / year)
Current delta deaths: 31852 (318 deaths / year)
Current surname count: 2044
Current partner percentage: 51
Net population per year: 168

2018-03-08 15:59:24.385102
Time taken: 13 seconds (132 milliseconds / year)
Current population size at year 200: 32645
Current mean age: 25
Current new people: 113218
Current deaths: 90573
Current delta new people: 64468 (644 people / year)
Current delta deaths: 58721 (587 deaths / year)
Current surname count: 1085
Current partner percentage: 52
Net population per year: 57

2018-03-08 15:59:41.181756
Time taken: 16 seconds (167 milliseconds / year)
Current population size at year 300: 39817
Current mean age: 25
Current new people: 191768
Current deaths: 161951
Current delta new people: 78550 (785 people / year)
Current delta deaths: 71378 (713 deaths / year)
Current surname count: 806
Current partner percentage: 53
Net population per year: 71

2018-03-08 16:00:02.352846
Time taken: 21 seconds (211 milliseconds / year)
Current population size at year 400: 47657
Current mean age: 25
Current new people: 286314
Current deaths: 248657
Current delta new people: 94546 (945 people / year)
Current delta deaths: 86706 (867 deaths / year)
Current surname count: 659
Current partner percentage: 50
Net population per year: 78

2018-03-08 16:00:28.200274
Time taken: 25 seconds (257 milliseconds / year)
Current population size at year 500: 58446
Current mean age: 25
Current new people: 400610
Current deaths: 352164
Current delta new people: 114296 (1142 people / year)
Current delta deaths: 103507 (1035 deaths / year)
Current surname count: 578
Current partner percentage: 53
Net population per year: 107

Now to why I'm here on Code-Review in the first place with this monster. I feel I have whittled down some of the major performance issues, but fresh eyes tend to see things that I miss. I'm also looking to make things more Pythonic and readable. One thing I originally forgot to mention explicitly, though it does appear in the code, is that I want to have a more robust plotting scheme, probably utilizing subplots. I haven't quite gotten around to making that happen, so if you see a better way, please share!

I don't claim to be a python guru by any means, but StackOverflow goes a long way to help learn a new language. Feel free to discuss anything you find odd, and if anything is unclear as to what it's doing, feel free to ask.

In terms of readability, I am aware that one of the things PEP8 specifies is snake_case but I'm originally a developer and old habits die hard (and I hate typing underscores). I also added a bunch of comments explaining my rationale behind why things appear as they do, which I feel gums things up a little, but I understand that everyone is not in my head.

Some side notes about timings below:

$ python -m timeit -s "import numpy as np" "np.random.randint(0, 1000)"
10000000 loops, best of 3: 0.164 usec per loop
$ python -m timeit -s "import random" "random.randint(0, 1000)"
1000000 loops, best of 3: 1.52 usec per loop
$ python -m timeit -s "import random" "random.uniform(0, 1)"
1000000 loops, best of 3: 0.489 usec per loop
$ python -m timeit -s "import numpy as np" "np.random.uniform(0, 1)"
1000000 loops, best of 3: 0.219 usec per loop

timeit.timeit("rand(0, 1000)", setup="from numpy.random import uniform as rand")
timeit.timeit("rand(0, 1000)", setup="from random import uniform as rand")

Edit from the future: I've converted this over to a C# project and posted it on github

  • \$\begingroup\$ Note, this should be a complete project that can be run on your local system. It does, however, hit a single thread pretty hard, so keep that in mind. \$\endgroup\$
    – Drise
    Mar 8, 2018 at 23:30
  • \$\begingroup\$ You missed the zeroth assumption - surnames are inherited (not true in many cultures, and we're explicitly not interested in those). \$\endgroup\$ Mar 9, 2018 at 11:29
  • \$\begingroup\$ @TobySpeight By inherited, you mean from parent to child? self.surname = parent.surname? \$\endgroup\$
    – hjpotter92
    Mar 9, 2018 at 11:37
  • \$\begingroup\$ Yes - that doesn't happen in Iceland, for example. I'm not claiming that's a problem - it's just outside the applicability of the research. And this problem is an interesting thing to simulate, so thank you for sharing it with us! \$\endgroup\$ Mar 9, 2018 at 11:51
  • 1
    \$\begingroup\$ I don't quite understand why C++ draws you away from snake_case; I'd have thought the opposite, given that that's the Standard Library's identifier convention! \$\endgroup\$ Mar 9, 2018 at 11:54

2 Answers 2


In terms of readability, I am aware that one of the things PEP8 specifies is snake_case

Since you're aware of PEP8, ignoring the naming, there are still a few more points to consider from the guide, such as:

  1. import order of modules
  2. Writing comments
  3. Having docstrings (I would heavily suggest this)

In a lot of places, I notice that you have snippets like

if self.gender == female:
[x for x in people if x.gender == person.male and x.isPartnerEligible()]

etc., despite your self and x having available the following definition:

def isMale(self):
    return self.gender == male

Instead of defining male and female as global integers, make use of an Enum which has been backported for older python versions. Or, you can just have True/False value (since it is a boolean in your simulated world). This helps with the helpers as well:

def isMale(self):
    return self.gender is MALE  # MALE might be True or False

def isFemale(self):
    return not self.isMale()

Variable names like haveChild, gainChild are a bit confusing (might be just me), use a more familiar terminology; for eg. isPregnant, giveBirth etc.

  • \$\begingroup\$ As a (micro?) optimization, I noticed function calls take orders of magnitude longer than just checking the value directly. It pained me to switch from p.isMale() to p.gender == male. Some stats (on a dummy class): checking directly t.testVal == 5 yields 0.0534 usec, whereas checking indirectly t.checkVal() yields 0.235 usec. Since this was a high traffic function, the nearly 5x speed up was needed. \$\endgroup\$
    – Drise
    Mar 9, 2018 at 15:28
  • 1
    \$\begingroup\$ @Drise: so instead of an p.is_male() method, you could simply define p.is_male as an attribute instead of having a gender attribute. You could even set p.is_female = not p.is_male for ease of use. \$\endgroup\$ Mar 13, 2018 at 8:11
  • \$\begingroup\$ That sounds like a good plan. \$\endgroup\$
    – Drise
    Mar 13, 2018 at 15:09

One thing that will speed things up a lot is to not use numpy without arrays. Specifically, when you use nprandint, you are using it to generate one random number. Instead, you should use random.randint, which is about 10% faster. Similarly, for uniform random is 3x faster for uniform, and 5x faster for exponential.

I got these timings as follows:

timeit.timeit('rand(0,1000)', setup='from numpy.random import uniform as rand')
timeit.timeit('rand(0,1000)', setup='from random import uniform as rand')
  • 1
    \$\begingroup\$ Oddly enough, when I did some of my performance testing (timeit on terminal) I found that numpy was significantly faster. I don't think it would just be present on my machine, so if that's the case, then that's very interesting. I'll do some investigation tomorrow and report my findings. Also, does python have a built in exponential distribution? \$\endgroup\$
    – Drise
    Mar 9, 2018 at 3:16
  • 1
    \$\begingroup\$ So I started with randint, timeit reports that np.random.randint(0, 1000) takes 0.164 usec, whereas random.randint(0, 1000) takes 1.52 usec. For uniform, np.random.uniform(0, 1) takes 0.219 usec, whereas random.uniform(0, 1) takes 0.489 usec. I don't see any increase using random, which is why I converted to using numpy's random facilities. \$\endgroup\$
    – Drise
    Mar 9, 2018 at 15:37
  • \$\begingroup\$ random.expovariate gives exponential distribution. Also, how were you testing? It is very easy to misuse timeit and get wrong results. \$\endgroup\$ Mar 9, 2018 at 17:01
  • \$\begingroup\$ I've moved my responses to the last section of the question for better formatting. \$\endgroup\$
    – Drise
    Mar 9, 2018 at 18:21
  • 1
    \$\begingroup\$ yeah, my timings have no relation to yours at all. That's bizzare \$\endgroup\$ Mar 9, 2018 at 19:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.