AM115 G-HW2 Q1
from matplotlib import pyplot as plt
from matplotlib import animation
import numpy as np
import pandas as pd
import seaborn as sns
from numba import jit
%matplotlib inline
R=1
K=50
dimP=10*K
mu=R/K
lambda_var=R*(1-(1/K))
print(lambda_var)
0.98
@jit
def sim_population(nsteps):
N=np.zeros(nsteps)
N[0]=3
t=0
dt=10**(-3)
while t < nsteps-1 and N[t] > 0:
r=np.random.rand(1)
if r < lambda_var * N[t] * dt:
N[t+1]=N[t]+1
elif r < (lambda_var*N[t]+mu*N[t]*(N[t]-1)) * dt:
N[t+1] = N[t]-1
else:
N[t+1] = N[t]
t=t+1
return N
realizations = np.array([
sim_population(10000) for _ in range(5000)
])
realizations_df = pd.DataFrame(realizations)
for i in range(0, 10000, 1000):
pd.Series(realizations[:, i + 99]).plot.hist(
bins=30, figsize=(9, 1), xlim=(0, 85),
title=f'Distribution after step {i+99}'
)
plt.show()
realzn_05 = realizations_df.quantile(0.05)
realzn_95 = realizations_df.quantile(0.95)
step_ts = realizations_df.columns
ax = realizations_df.mean().plot(
figsize=(9, 4), lw=4,
title='blue line: mean; green line: median; blue shades: 5%-95% conf interval',
xlabel='Time Step', ylabel='Population'
)
ax.plot(step_ts, realizations_df.median(), lw=2, alpha=0.8)
ax.plot(step_ts, realzn_05, lw=1, color='k')
ax.plot(step_ts, realzn_95, lw=1, color='k')
ax.fill_between(step_ts, realzn_05, realzn_95, alpha=0.3)
realizations_df.skew().plot(
figsize=(9, 4), lw=2, title='Distribution Skewness',
xlabel='Time Step'
)