Get a feeling for how often it occurs that the oldest person alive died. ``` moose@pc08 ~ $ python birth.py -n 10000 > result.txt moose@pc08 ~ $ cat result.txt Day 1: Yes (Died in age of 70.01 years) Day 2: Yes (Died in age of 70.01 years) Day 3: Yes (Died in age of 70.01 years) Day 4: Yes (Died in age of 70.01 years) Day 25549: Yes (Died in age…

birth.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import random

""" Another day passes. Some people die, some survive, but everyone gets older.
	@param people dead people get removed from this list
	@return list of all people who died
"""
def updatePeople(people):
	deadAges = []
	deadList = []
	for i in xrange(len(people)):
		people[i][0] += 1
		people[i][1] -= 1
		if people[i][1] <= 0:
			deadList.append(i)

	deadList.sort()
	for i in reversed(deadList):
		deadAges.append(people.pop(i))

	return deadAges

""" Check if the oldest person died
	@param people all people that are still alive
	@param dead all people that just died today
	@return true iff no person is older than the oldes dead person
"""
def didOldestPersonDie(people, dead):
	peopleAgeMax = max(people, key=lambda person: person[0])
	deadAgeMax = max(dead, key=lambda zombie: zombie[0])
	return deadAgeMax[0] >= peopleAgeMax[0]

if __name__ == "__main__":
	from argparse import ArgumentParser

	parser = ArgumentParser()

	# Add more options if you like
	parser.add_argument("-n", dest="n", type=int, default=1000, 
						help="How many people are alive at the start?")
	parser.add_argument("-mu", dest="mu", type=int, default=70*365, 
						help="How long does the average human live (in days)?")
	parser.add_argument("-sigma", dest="sigma", type=float, default=0.7, 
						help="standard deviation for live expectancy")
	parser.add_argument("-t", "--maxDays", type=int, dest="t", default=200*365, 
						help="maximum simulation time")
	args = parser.parse_args()


	n = args.n
	mu = args.mu
	sigma = args.sigma
	maxDays = args.t

	people = []
	for i in xrange(n):
		# TODO: You should probably choose a better distributation than
		# normal distribution
		ttl = random.normalvariate(mu, sigma)
		currentAge = random.normalvariate(mu, sigma)

		while ttl < currentAge:
			ttl = random.normalvariate(mu, sigma)
			currentAge = random.normalvariate(mu, sigma)

		ttl -= currentAge

		people.append([currentAge, ttl]) # current age and time to live

	day = 0
	while len(people) > 0 and maxDays > day:
		# another day
		day += 1

		# some die
		deadAges = updatePeople(people)

		# and some get born
		# TODO: This should probably not be directly related to dead people
		for i in xrange(len(deadAges)):
			ttl = random.normalvariate(mu, sigma)
			people.append([0, ttl])

		if len(deadAges) > 0 and len(people) > 0 and didOldestPersonDie(people, deadAges):
			print("Day %i: Yes (Died in age of %.2f years)" % 
					(day, float(max(deadAges, key=lambda dead: dead[0])[0])/365))