# Simulating surname convergence in a population

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
else:
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.partnerEligibleAge > self.maxAge:
self.partnerEligibleAge = self.maxAge
else:
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

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 \

# 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())
SERIALIZE = False

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):
distribution.append(int(getRandomAgeTo80()))

for i in range(2000):
distribution.append(int(getRandomAgeFrom80()))

# 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

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:
p.gainChild()

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.
random.shuffle(mPartners)
random.shuffle(fPartners)

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

for m, f in pairs:
m.gainPartner(f)
f.gainPartner(m)

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
people.remove(p)

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
setupAgeDistribution(Person.maxAgeDistribution)
random.shuffle(Person.maxAgeDistribution)
random.shuffle(Person.maxChildDistribution)

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()

try:
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)
break

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

sys.exit(0)

# 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
stats.timeTaken.append(timeTaken)

population = len(people)
stats.population.append(population)

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

stats.newPeople.append(stats.newPersonCount)
peopleDelta = stats.newPersonCount - stats.oldPeople
stats.oldPeople = stats.newPersonCount
stats.deltaNewPeople.append(peopleDelta)

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

stats.deaths.append(Person.deaths)
deathsDelta = Person.deaths - stats.oldDeaths
stats.oldDeaths = Person.deaths

netPop = (peopleDelta / statYear) - (deathsDelta / statYear)
stats.netPopDelta.append(netPop)

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:

if stats:
with open("stats%s.dat" % id, "rb") as 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.plot(stat)
labels.append(label)

plt.grid(True)

plt.legend(labels, ncol=2, loc='upper center', bbox_to_anchor=[0.5, 1.1], columnspacing=1.0,

# plt.yscale('log')
plt.show()

# 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
main()


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.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.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")
0.1250929832458496
timeit.timeit("rand(0, 1000)", setup="from random import uniform as rand")
0.47645998001098633

• 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. – Drise Mar 8 '18 at 23:30
• You missed the zeroth assumption - surnames are inherited (not true in many cultures, and we're explicitly not interested in those). – Toby Speight Mar 9 '18 at 11:29
• @TobySpeight By inherited, you mean from parent to child? self.surname = parent.surname? – hjpotter92 Mar 9 '18 at 11:37
• 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! – Toby Speight Mar 9 '18 at 11:51
• 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! – Toby Speight Mar 9 '18 at 11:54

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
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.

• 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. – Drise Mar 9 '18 at 15:28
• @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. – Eric Duminil Mar 13 '18 at 8:11
• That sounds like a good plan. – Drise Mar 13 '18 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')
0.7485988769913092
timeit.timeit('rand(0,1000)', setup='from random import uniform as rand')
0.2676771300029941

• 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? – Drise Mar 9 '18 at 3:16
• 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. – Drise Mar 9 '18 at 15:37
• random.expovariate gives exponential distribution. Also, how were you testing? It is very easy to misuse timeit and get wrong results. – Oscar Smith Mar 9 '18 at 17:01
• I've moved my responses to the last section of the question for better formatting. – Drise Mar 9 '18 at 18:21
• yeah, my timings have no relation to yours at all. That's bizzare – Oscar Smith Mar 9 '18 at 19:23