Load ilmtools, and any other necessary libraries
import ilmtools as ilm
import time
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Generate a population (n=100) with a uniform spatial distribution with bounds x=(0,10), y=(0,10)
n=100
xbounds=[0,10]
ybounds=[0,10]
pop = ilm.unif_indvs(n, xbounds[0], xbounds[1], ybounds[0], ybounds[1], 'example_pop')
Propagate an infectious disease over 10 days through the previously defined population with alpha=2, beta=15, from a single initial infection
initial_infections = 1
timesteps = 10
true_alpha = 2
true_beta = 15
infect_db = ilm.si_model(pop, initial_infections, timesteps, true_alpha, true_beta)
Visualize outbreak at any time point, e.g. day 10
ilm.plot_si(pop, infect_db, 10)
Visualize epidemic curve
ilm.epi_curve(infect_db)
Define prior distributions for SI model parameters alpha and beta, and perform Bayesian inference using MCMC (initial run of 1000 iterations), and report the time required for computation
def prior_alpha(x):
return ilm.pdf_gamma(2, 1, x)
init_alpha1 = np.random.uniform(0, 10)
def prior_beta(x):
return ilm.pdf_gamma(10, 1, x)
init_beta1 = np.random.uniform(10, 20)
niter1=1000
transition_cov = [[1, 0], [0,3]]
start_time1=time.time()
mcmc1 = ilm.si_infer(pop, infect_db, prior_alpha, init_alpha1, prior_beta, init_beta1, niter1, transition_cov)
required_time1=time.time()-start_time1
required_time1
1535.2476561069489
Estimate optimal proposal covariance matrix (Rosenthal), and run 1000 more MCMC iterations, initiated at the last value in the initial run
opt_cov=np.cov(mcmc1.alpha, mcmc1.beta)*((2.38**2.)/2.)
niter2=1000
init_alpha2 = mcmc1.alpha[999]
init_beta2 = mcmc1.beta[niter1-1]
start_time2=time.time()
mcmc2 = ilm.si_infer(pop, infect_db, prior_alpha, init_alpha2, prior_beta, init_beta2, niter2, opt_cov)
required_time2=time.time()-start_time2
required_time2
1519.4822149276733
Now, produce traceplots of alpha, beta, and posterior density
plt.plot(mcmc2.alpha)
[<matplotlib.lines.Line2D at 0x11043b9d0>]
plt.plot(mcmc2.beta)
[<matplotlib.lines.Line2D at 0x1104c0fd0>]
plt.plot(mcmc2.density)
[<matplotlib.lines.Line2D at 0x111cac610>]
np.mean(mcmc2.alpha), np.mean(mcmc2.beta)
(1.3363672429039526, 11.219127783512988)
plt.plot(mcmc2.alpha, mcmc2.beta)
[<matplotlib.lines.Line2D at 0x112281b50>]