In [4]:
import sys
import os
import math
import importlib
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import random
%matplotlib inline
In [36]:
def simulate(initialRatio):
rand = random.Random(1337)
nprand = np.random.RandomState(1337)
populationSize = 5000 # small town
# model # of children as Poisson:
# 2.5 avg children, and for Poisson, mean=variance
# note that the 2.5 number is for present day - that's a conservative choice for it
# since in past it's much higher:
# https://ourworldindata.org/fertility-rate#the-number-of-children-per-woman-over-the-very-long-run
numChildrenMeanVariance = 2.5
# initialize the initial generation
initialSeeds = rand.sample(range(populationSize), int(populationSize * initialRatio))
population = [False]*populationSize
for i in initialSeeds:
population[i] = True
result = [(populationSize, len(initialSeeds))]
# start simulating
generations = 10
for generation in range(generations):
# shuffle the array
newpopulation = []
rand.shuffle(population)
popWithTrait = 0
# for each (randomly shuffled) pair:
endval = len(population) & ~1
for i in range(0, endval, 2):
# generate number of children
# if any parent has it, child has it
hasIt = population[i] or population[i+1]
numChildren = int(nprand.poisson(numChildrenMeanVariance))
newpopulation = newpopulation + [hasIt]*numChildren
if hasIt:
popWithTrait += numChildren
print(len(newpopulation), popWithTrait)
result.append((len(newpopulation), popWithTrait))
population = newpopulation
return result
In [37]:
ratios = [0.01, 0.05, 0.10, 0.15, 0.20]
simresults = []
for ratio in ratios:
print("starting", ratio)
simresults.append(simulate(ratio))
In [42]:
for i in range(len(ratios)):
simresult = simresults[i]
ratio = [a[1]/a[0] for a in simresult]
plt.plot(list(range(len(simresult))), ratio, label=str(int(ratios[i]*100)) + "%")
plt.legend()
Out[42]:
In [49]:
def simulate_const(initialRatio):
rand = random.Random(1337)
populationSize = 2726150 # matches 1809 pop
# initialize the initial generation
initialSeeds = rand.sample(range(populationSize), int(populationSize * initialRatio))
population = [False]*populationSize
for i in initialSeeds:
population[i] = True
result = [(populationSize, len(initialSeeds))]
# start simulating
generations = 10
for generation in range(generations):
# shuffle the array
newpopulation = [False]*len(population)
rand.shuffle(population)
popWithTrait = 0
# for each (randomly shuffled) pair:
endval = len(population) & ~1
for i in range(0, endval, 2):
# generate number of children
# if any parent has it, child has it
hasIt = population[i] or population[i+1]
if hasIt: # changed from above!
newpopulation[i] = newpopulation[i+1] = True
popWithTrait += 2
print(len(newpopulation), popWithTrait)
result.append((len(newpopulation), popWithTrait))
population = newpopulation
return result
simresults_const = []
for ratio in ratios:
print("starting", ratio)
simresults_const.append(simulate_const(ratio))
In [50]:
for i in range(len(ratios)):
simresult = simresults_const[i]
ratio = [a[1]/a[0] for a in simresult]
plt.plot(list(range(len(simresult))), ratio, label=str(int(ratios[i]*100)) + "%")
plt.legend()
Out[50]:
In [51]:
# Just for fun: what if only one child gets the trait? (matches Guardian)
def simulate_onlyone(initialRatio):
rand = random.Random(1337)
populationSize = 2726150 # matches 1809 pop
# initialize the initial generation
initialSeeds = rand.sample(range(populationSize), int(populationSize * initialRatio))
population = [False]*populationSize
for i in initialSeeds:
population[i] = True
result = [(populationSize, len(initialSeeds))]
# start simulating
generations = 10
for generation in range(generations):
# shuffle the array
newpopulation = [False]*len(population)
rand.shuffle(population)
popWithTrait = 0
# for each (randomly shuffled) pair:
endval = len(population) & ~1
for i in range(0, endval, 2):
# generate number of children
# if any parent has it, child has it
hasIt = population[i] or population[i+1]
if hasIt: # changed from above!
newpopulation[i] = True
popWithTrait += 1
print(len(newpopulation), popWithTrait)
result.append((len(newpopulation), popWithTrait))
population = newpopulation
return result
simresults_onlyone = []
for ratio in ratios:
print("starting", ratio)
simresults_onlyone.append(simulate_onlyone(ratio))
In [52]:
for i in range(len(ratios)):
simresult = simresults_onlyone[i]
ratio = [a[1]/a[0] for a in simresult]
plt.plot(list(range(len(simresult))), ratio, label=str(int(ratios[i]*100)) + "%")
plt.legend()
Out[52]:
In [ ]: