#!/usr/bin/env python # coding: utf-8 # In[4]: import sys import os import math import importlib import matplotlib import matplotlib.pyplot as plt import numpy as np import random get_ipython().run_line_magic('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() # 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() # 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() # In[ ]: