Revise timeseries slides

models/exercise/plot_two_strain_ts.py

@@ -23,8 +23,8 @@
 epsilon = 0.1
 gamma = np.array([1, 1])/7.0
 mu = 1/(10*365.0)
-alpha = np.array([1., 1.])
-a = np.array([1., 1.])
+alpha = np.array([0.5, 1.])
+a = np.array([1., 2.5])
 omega = 2*np.pi/365.
 obs_sd = 0.01
 

models/exercise/plot_two_strain_ts.py~

@@ -0,0 +1,89 @@
+###################################################################																   #
+# Basic plot for two-strain SIR model:
+# Time series given some initial conditions
+####################################################################
+
+import sys
+import csv
+import numpy as np
+import matplotlib as mpl
+mpl.use('TkAgg')
+from matplotlib.font_manager import FontProperties 
+import matplotlib.pyplot as plt
+from two_strain import *
+
+# Run parameters
+run_num = 1  # sys.argv[1]
+end_time = 100*365
+output_interval = 1.0
+step_size = 0.1
+
+# Strain parameters, including initial conditions
+beta = np.array([5, 5])/7.0
+epsilon = 0.1
+gamma = np.array([1, 1])/7.0
+mu = 1/(10*365.0)
+alpha = np.array([0.5, 1.])
+a = np.array([1., 1.5])
+omega = 2*np.pi/365.
+obs_sd = 0.01
+
+NSS = 0.2
+NIS = 1e-3
+NRS = 0.02
+NRI = 0.0
+NSI = 1e-3
+NSR = 0.02
+NIR = 0.0
+
+# Organize and run simulation
+params = np.array([gamma, mu, alpha, a, omega, beta, epsilon])
+SI = np.array([NSS, NIS, NRS, NRI, NSI, NSR, NIR])
+ic = np.array([NSS, NIS, NRS, NRI, NSI, NSR, NIR, 1-np.sum(SI)])
+output = run_two_strain(end_time, output_interval, step_size, params, ic)
+
+# Save output (NIS+NIR, NSI+NRI) to csv and plot
+infecteds = np.asarray([output[:, 1] + output[:, 6], output[:, 3] + output[:, 4]])
+times = np.arange(0,infecteds.shape[1])
+infecteds_t = np.vstack((times, infecteds))
+filename = 'infecteds_' + str(run_num) + '.csv'
+with open(filename, 'w') as csvfile:
+    writer = csv.writer(csvfile)
+    writer.writerow(['times', 'I1', 'I2'])
+    writer.writerows(infecteds_t.T)
+
+# Add observation error if present
+if obs_sd > 0:
+    errors = np.random.normal(1, obs_sd, infecteds.shape)
+    infecteds_obs = infecteds*errors
+    filename = 'infecteds_obs_' + str(run_num) + '.csv'
+    with open(filename, 'w') as csvfile:
+        writer = csv.writer(csvfile)
+        writer.writerow(['times', 'I1', 'I2'])
+        writer.writerows(infecteds_t.T)
+
+plt.subplot(3, 1, 1)
+plt.plot(output[:, 0], 'b-', label=r'$N_{SS}$')
+plt.plot(output[:, 2], 'g-', label=r'$N_{RS}$')
+plt.plot(output[:, 5], 'r-', label=r'$N_{SR}$')
+plt.plot(output[:, 7], 'c-', label=r'$N_{RR}$')
+plt.xlabel('Time')
+plt.ylabel('Uninfected')
+plt.legend(loc=1, prop=FontProperties(size='smaller'))
+plt.subplot(3, 1, 2)
+plt.plot(output[:, 1], 'b-', label=r'$N_{IS}$')
+plt.plot(output[:, 6], 'g-', label=r'$N_{IR}$')
+plt.plot((output[:, 1]+a[0]*output[:, 6]), 'r-', label=r'$I_1$')
+plt.xlabel('Time')
+plt.ylabel('Infected 1')
+plt.legend(loc=1, prop=FontProperties(size='smaller'))
+plt.subplot(3, 1, 3)
+plt.plot(output[:, 4], 'b-', label=r'$N_{SI}$')
+plt.plot(output[:, 3], 'g-', label=r'$N_{RI}$')
+plt.plot((output[:, 4]+a[1]*output[:, 3]), 'r-', label=r'$I_2$')
+plt.xlabel('Time')
+plt.ylabel('Infected 2')
+plt.legend(loc=1, prop=FontProperties(size='smaller'))
+plt.savefig("time_series_" + str(run_num) + ".png")
+plt.show()
+plt.close()

timeseries/slides.html

@@ -47,9 +47,9 @@ <h2 class="title">Fitting mechanistic models</h2>
 
 				<section>
 					<h3>A general approach to model fitting</h3>
-	                <p>Pick a model and some parameters
+	                <p>Pick a model, including parameter values
 					<p>Evaluate how likely the observed data are, given the model
-					<p>Tweak the parameters/model to make the observations more likely
+					<p>Tweak the model to make the observations more likely
 					<p>Is this model superior to other models?
 				</section>
 
@@ -87,7 +87,7 @@ <h3>Maximum likelihood (ML) inference</h3>
 
 				<section>
 					<h3>Likelihood in timeseries models</h3>
-					<p>Observed trajectory $D=(t_0,...,t_n)$ depends on unknown parameter $\theta$
+					<p>Observed trajectory $D=(t_0,...,t_n)$ depends on unknown parameter(s) $\theta$
 					<p>Probability density function of $D$ is $f_{\theta}$
 					<p>$$L_D(\theta)=f_{\theta}(D)$$
                     <p>Problem: These data aren't independent.
@@ -217,8 +217,7 @@ <h3>Methods for Bayesian integration</h3>
 				</section>
 
 				<section>
-					<h3><a href="http://cobeylab.uchicago.edu/mcmc/intro/">Introduction to Monte Carlo integration</a></h3>
-                    <h3><a href="https://chi-feng.github.io/mcmc-demo/app.html">(and an even better one)</a></h3>
+                    <h3><a href="https://chi-feng.github.io/mcmc-demo/app.html">Introduction to Monte Carlo integration</a></h3>
 				</section>
 
 				<section>