CoCalc -- 3044.sagews

All published worksheets from http://sagenb.org

Project: sagenb.org published worksheets

Views: ¹⁶⁸⁷⁴²
Image: ubuntu2004

birth_rate=0.5
catastrophe_birth_rate=birth_rate*.6
death_rate=.1
catastrophe_death_rate=death_rate*1.25

catastrophe_birth_rate-catastrophe_death_rate

0.175000000000000

def f(x):
    return x*(1+birth_rate-death_rate)

f(100)

140.000000000000

# Construct a list of populations
pop=[100]
for i in [1..20]:
    pop.append(f(pop[-1]))
list_plot(pop,plotjoined=True)

import scipy.stats
bernoulli_dist=scipy.stats.bernoulli(.04)

[bernoulli_dist.rvs() for _ in range(100)]

[0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0]

def f(x):
    if bernoulli_dist.rvs() == 0:
        return x*(1+birth_rate-death_rate)
    else:
        return x*(1+catastrophe_birth_rate-catastrophe_death_rate)

# Construct a list of populations
pop=[]
for j in range(10):
    pop1=[100]
    for i in [1..10]:
        pop1.append(f(pop1[-1]))
    pop.append(pop1)
pop
sum([list_plot(p,plotjoined=True) for p in pop])

html.table(zip(*pop))

100	100	100	100	100	100	100	100	100	100
127.827645690169	135.520933632463	137.598362118643	146.756495889536	137.817422222928	141.616661273534	154.887101933696	133.559314024952	123.162535541558	134.518193117627
193.656462523654	195.124565604085	188.301976204673	195.940469890338	201.670383054670	202.523231082640	225.017691317239	179.532989422411	165.242745115650	171.270971480620
278.186033319413	251.091511973860	248.949163043783	264.111924626650	274.795973898064	266.724791651686	292.896474693648	280.739257700192	233.497841126388	240.186656820459
379.571379473715	407.208464761504	326.291110873411	422.434059289027	423.430715009156	347.250210171534	389.019961620823	381.507603995059	342.666366524243	360.545716285012
525.922473848415	627.475005622336	492.725360205785	583.596445104336	595.463035594874	472.323205179581	567.513303470037	498.423155973587	488.645286466954	504.467873129729
692.307591306724	846.927654108813	693.820489590030	861.423399128817	782.837032636771	717.023911925502	763.962081067512	661.399130843303	698.244140377754	676.071753637208
990.181390260939	1235.24158758905	976.558554730306	1080.33669954585	1016.25056161721	885.228502450429	1057.22584496739	933.627086292145	904.978834041131	939.834546732937
1492.60780650675	1659.43005112830	1333.21067166941	1541.57854110041	1480.63423389602	1239.02002688985	1404.86535521309	1413.76874887570	1207.80327101258	1324.76188272373
2013.97729779576	2365.59686511154	2011.48019128964	2215.01639980188	2258.10745430564	1735.14080468801	1875.73242004762	2134.41054764020	1686.12008700979	1994.51662279211
3195.75954927205	3404.99830369928	3065.37838479668	2977.49752885509	3296.11905350034	2452.29834678994	2618.59391245239	2931.69085898117	2129.04375314629	2764.55817803471

Now let's try varying the birth and death rates with a normal distribution.

import scipy.stats
birth_rate_dist=scipy.stats.norm(0.5, 0.03)
death_rate_dist = scipy.stats.norm(0.1,0.08)

A couple of sample birth rates, and a couple of sample death rates.

birth_rate_dist.rvs(10)

array([ 0.54570469,  0.52116238,  0.45963489,  0.49774404,  0.50283972,
        0.51396968,  0.48677194,  0.472691  ,  0.43738911,  0.45034876])

death_rate_dist.rvs(10)

array([ 0.2086821 ,  0.10122123,  0.02327865,  0.23890433,  0.15732268,
        0.03491447,  0.28749208,  0.12381044,  0.15774519,  0.02695242])

plot(birth_rate_dist.pdf, 0, 1) + plot(death_rate_dist.pdf, -1, 1, color='red')

def f(x):
    return x*(1+birth_rate_dist.rvs()-death_rate_dist.rvs())

# Construct a list of populations
pop=[]
num_trials=100
starting_population=100
num_years=10
for j in range(num_trials):
    pop1=[starting_population]
    for i in [1..num_years]:
        pop1.append(f(pop1[-1]))
    pop.append(pop1)

Here is a table of the first 5 populations. The "zip" thing is just to make the table so that years are the rows of the tables instead of the columns of the table (i.e., it takes the transpose of the matrix)

html.table(zip(*pop[:5]))

100	100	100	100	100
139.134935321599	139.783780237118	132.758965646092	136.043913588454	142.759376457081
198.649441546414	193.217165964711	188.680322531604	192.625898484487	186.628869680903
276.120552247139	273.086679191364	264.836102674358	298.910747293408	264.434000742601
411.644515127575	359.801312829274	363.267899724298	394.401556209514	360.723371546198
601.785508869635	483.616638874819	482.832579018686	516.565131105165	490.510125502479
904.841704316214	615.554911068894	633.202926806810	743.969532342385	675.076610863102
1260.55111306693	891.913482794109	880.105016325925	1083.68128196380	862.768570665592
1688.47864266620	1248.32487759989	1330.72493014679	1579.08216593386	1111.02703178465
2227.47063239747	1839.28881733030	1600.10588925458	2302.34649546089	1618.97999717236
3327.18222654253	2668.59704057534	2314.84331170542	3673.51746701349	2275.24788361626

sum([list_plot(p,plotjoined=True) for p in pop])

import matplotlib.pyplot as plt

This just makes the list so that year_population[i] gives the populations of year i for each simulation.

year_population=zip(*pop)

year_population[0] # population of year 0 for each run

(100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100)

year_population[1] # population of year 1 for each run

(139.134935321599, 139.783780237118, 132.758965646092, 136.043913588454, 142.759376457081, 133.885283691328, 148.226556697086, 126.880026606073, 135.213510039784, 139.348026285354, 143.361446137340, 131.892122031762, 122.650334572014, 140.953881760382, 144.516517650291, 130.275781408026, 126.860723078949, 135.737874297401, 139.198577789482, 142.167609886886, 149.849628815562, 129.762264477777, 133.992231284889, 113.730643434116, 153.617860781787, 146.458859444053, 135.485454612436, 135.158913430591, 154.466626599108, 143.498789205192, 145.000025311566, 136.326086989476, 146.915740733372, 152.426435575132, 141.399953689912, 131.161659751253, 137.572279965885, 155.298285260022, 143.875651923280, 141.189129655305, 127.738745381908, 141.167192657279, 123.604102652840, 128.464635079208, 136.057893655640, 149.931677109767, 134.982018945917, 138.349280702407, 136.203566829550, 138.424484858375, 134.889300473408, 139.103223125291, 153.777466820403, 142.100323852806, 139.731316397147, 136.459690226559, 135.153647256176, 131.435611984491, 132.750256388951, 134.637565798292, 155.269106184037, 142.569627882389, 135.575273030306, 137.815595191047, 138.062216570851, 145.619419602813, 135.859551097522, 142.271004103355, 152.879770872037, 141.766109677045, 148.658182453907, 135.744168406247, 138.734422552652, 140.530748720943, 144.861921955497, 139.251448893494, 131.690366379413, 125.803201859373, 133.702775946830, 144.789609547350, 149.377376182883, 143.931836949578, 144.800932588664, 138.784187344468, 140.635080844022, 144.832054733102, 127.805760943844, 131.742685359520, 145.745036981688, 153.343088906196, 137.082669418678, 124.506384690711, 136.690071305834, 135.518580133435, 136.933917243670, 141.569111490367, 129.911053551910, 124.994951777877, 134.881668660948, 132.030271871814)

year_population[-1] # population of the last year of each run

(3327.18222654253, 2668.59704057534, 2314.84331170542, 3673.51746701349, 2275.24788361626, 2523.37241443004, 2166.47346849146, 2200.10406436381, 2874.74998565683, 2930.29477936063, 2388.77660003879, 2560.92624274997, 2794.48209699705, 2869.99220751551, 3563.27264630089, 2430.78441392426, 3195.26250843154, 2388.17213894411, 2579.23041342220, 1672.60611810455, 3243.88895084611, 3298.42808360321, 2550.05699381109, 2672.02798664290, 3992.84471130146, 2157.00320555761, 2434.03566317628, 3374.29801416882, 3783.62146602039, 3110.68389676390, 2825.68402185595, 2524.09698940108, 2912.72505175626, 3305.69100667337, 2158.19564956439, 2950.25217987774, 2076.91430577845, 3444.79725610605, 3050.36489954344, 2981.57023381918, 2544.05129143295, 2892.77272475061, 2578.42324554406, 3065.71951641101, 2071.71685662934, 3055.78866285866, 2857.44189126427, 1598.97459467169, 3368.61743153454, 3305.23073493070, 3164.39705365245, 2490.03583587049, 3854.70541486060, 3138.12749605614, 2954.16346906145, 2887.78911100064, 3129.72703849244, 2587.44015894855, 3162.94993932421, 2089.97558844329, 3333.05754872896, 2946.13716052072, 4096.07599036791, 3344.56253850379, 2962.38981423110, 3320.81548475485, 4233.14640648488, 2693.10391564959, 3028.54972651927, 3126.48215050352, 3419.11313799296, 3032.06146867666, 2661.97036567872, 2505.95586250678, 2663.84827810512, 1965.47323694509, 2752.52446607466, 2238.99120104620, 3035.71825311736, 3174.96289121128, 3403.64880043544, 2844.79208332133, 3943.86304755817, 3104.53099363627, 2892.61776581906, 2853.64227258368, 2665.95155556990, 2275.29365276608, 3864.82892900084, 2564.90683774527, 3166.31001328293, 2702.47944171152, 2563.97218609187, 2776.73276596825, 2726.65236668962, 2379.11132306736, 3003.06771669461, 3104.46258181592, 3046.07938455707, 3699.65700956455)

plt.figure()
n,bins,_=plt.hist(year_population[-1],bins=8)
print "frequencies: ",n
print "bin boundaries: ",list(bins)
plt.title("Ending Population")
plt.xlabel("Population")
plt.ylabel("Frequency")
plt.savefig('test.png')

frequencies:  [ 2  9 19 22 26 13  5  4]
bin boundaries:  [1598.9745946716864, 1928.2460711483361, 2257.517547624986, 2586.7890241016357, 2916.0605005782854, 3245.3319770549351, 3574.6034535315848, 3903.8749300082345, 4233.1464064848842]

plt.figure()
bp=plt.boxplot(year_population, positions=range(len(year_population)))
plt.title("Populations")
plt.xlabel("Year")
plt.ylabel("Population")
plt.savefig('test.png')