Bugs in documentation reported by Rita

docs/quickstart.md

@@ -174,11 +174,11 @@ Setting up and simulating the model,
 # Define parameters, initial states and coordinates
 params={'alpha': np.array([0.05, 0.1, 0.2, 0.15]), 'gamma': 5, 'beta': 7}
 init_states = {'S': [606938, 1328733, 7352492, 2204478], 'S_v': 1e6, 'I_v': 2}
-coordinates={'age_groups': ['0-5','5-15', '15-65','65-120']}
+coordinates={'age_group': ['0-5','5-15', '15-65','65-120']}
 # Initialize model
-model = SIR_SI(states=init_states, parameters=params, coordinates=coordinates)
+model = ODE_SIR_SI(states=init_states, parameters=params, coordinates=coordinates)
 # Simulate the model
-model.sim(120)
+out = model.sim(120)
 print(out)
 ```
 

docs/workflow.md

@@ -19,7 +19,7 @@ I typically place all my dependencies together at the top of my script. However,
 ```
 import numpy as np
 import pandas as pd
-from matplotlib.pyplot import plt
+import matplotlib.pyplot as plt
 ```
 
 ### Load the dataset
@@ -33,16 +33,18 @@ n_{cases}(t) = \exp \Big( t *  \dfrac{\log 2}{t_d} \Big)
 We'll assume the first case was detected on December 1st, 2022 and data was collected on every weekday until December 21st, 2022. Then, we'll add observational noise to the synthetic data. For count based data, observational noise is typically the result of a poisson or negative binomial proces, depending on the occurence of overdispersion. For a poisson proces, the variance in the data is equal to the mean: {math}`\sigma^2 = \mu`, while for a negative binomial proces the mean-variance relationship is quadratic: {math}`\sigma^2 = \mu + \alpha \mu^2`. For this example we'll use the negative binomial distribution with a dispersion factor of `alpha=0.03`, which was representative for the data used during the COVID-19 pandemic in Belgium. Note that for {math}`\alpha=0`, the variance of the negative binomial distribution is equal to the variance of the poisson distribution.
 
 ```
-# Parameters
-alpha = 0.03 # Overdispersion
-t_d = 10 # Doubling time
+# Parameters 
+alpha = 0.03    # Overdispersion
+t_d = 10         # Exponential doubling time
 # Sample data
 dates = pd.date_range('2022-12-01','2023-01-21')
 t = np.linspace(start=0, stop=len(dates)-1, num=len(dates))
-y = np.random.negative_binomial(1/alpha, (1/alpha)/(np.exp(t*np.log(2)/td) + (1/alpha)))
+y = np.random.negative_binomial(1/alpha, (1/alpha)/(np.exp(t*np.log(2)/t_d) + (1/alpha)))
 # Place in a pd.Series
 d = pd.Series(index=dates, data=y, name='CASES')
-# Data collection only on weekdays only
+# Index name must be date for calibration to work
+d.index.name = 'date'
+# Data collection only on weekdays
 d = d[d.index.dayofweek < 5]
 ```
 
@@ -106,7 +108,7 @@ class ODE_SIR(ODE):
     """
 
     states = ['S','I','R']
-    parameter_names = ['beta','gamma']
+    parameters = ['beta','gamma']
 
     @staticmethod
     def integrate(t, S, I, R, beta, gamma):
@@ -165,7 +167,7 @@ from the above equations we can deduce that the SSEs's use is only appropriate w
 |------------------------------|--------------------------------------------|
 | Gaussian                     | {math}`\sigma^2 = c`                       |
 | Poisson                      | {math}`\sigma^2 = \mu`                     |
-| Quasi-Poison                 | {math}`\sigma^2 = \alpha * \mu`            |
+| Quasi-Poisson                | {math}`\sigma^2 = \alpha * \mu`            |
 | Negative Binomial            | {math}`\sigma^2 = \mu + \alpha * \mu^2`    |
 
 The following snippet performs the above procedure on our synthetic dataset. 
@@ -251,7 +253,7 @@ if __name__ == '__main__':
     # Update beta with the calibrated value
     model.parameters.update({'beta': theta[0]})
     # Simulate the model
-    out = model.sim([start_date, end_date])
+    out = model.sim([d.index.min(), d,index.max()])
     # Visualize result
     fig,ax=plt.subplots()
     ax.plot(out['date'], out['I'], color='red', label='Infectious')
@@ -301,7 +303,7 @@ if __name__ == '__main__':
     # Perturbate previously obtained estimate
     ndim, nwalkers, pos = perturbate_theta(theta, pert=[0.10,], multiplier=multiplier_mcmc, bounds=bounds)
     # Usefull settings to retain in the samples dictionary (no pd.Timestamps or np.arrays allowed!)
-    settings={'start_calibration': start_date.strftime("%Y-%m-%d"), 'end_calibration': end_date.strftime("%Y-%m-%d"),
+    settings={'start_calibration': d.index.min().strftime("%Y-%m-%d"), 'end_calibration': end_date.strftime("%Y-%m-%d"),
               'n_chains': nwalkers, 'starting_estimate': list(theta)}
     # Run the sampler
     sampler = run_EnsembleSampler(pos, n_mcmc, identifier, objective_function,
@@ -360,7 +362,7 @@ As demonstrated in the quickstart example, the `xarray` containing the model out
 
 ```
 # Simulate model
-out = model.sim([start_date, end_date+pd.Timedelta(days=28)], N=50, samples=samples_dict, draw_fcn=draw_fcn, processes=processes)
+out = model.sim([d.index.min(), d.index.max()+pd.Timedelta(days=28)], N=50, samples=samples_dict, draw_fcn=draw_fcn, processes=processes)
 # Add negative binomial observation noise
 out = add_negative_binomial_noise(out, dispersion)
 # Visualize result
@@ -398,18 +400,18 @@ This simple function reduces the parameter `param` with a factor two if the simu
 We will apply this function to the infectivity parameter `beta` by using the `time_dependent_parameters` argument of `sim()`.
 ```
 # Attach the additional arguments of the time-depenent function to the parameter dictionary
-params.update({'start_measures': '2023-01-21'})
+model.parameters.update({'start_measures': '2023-01-21'})
 
 # Initialize the model with the time dependent parameter funtion
-model_with = ODE_SIR(states=init_states, parameters=params,
+model_with = ODE_SIR(states=init_states, parameters=model.parameters,
                      time_dependent_parameters={'beta': lower_infectivity})
 ```
 
 Next, we compare the simulations with and without the use of our time-dependent model parameter function.
 
 ```
 # Simulate the model
-out_with = model_with.sim([start_date, end_date+pd.Timedelta(days=2*28)], N=50, samples=samples_dict, draw_fcn=draw_fcn, processes=processes)
+out_with = model_with.sim([d.index.min(), d.index.max()+pd.Timedelta(days=2*28)], N=50, samples=samples_dict, draw_fcn=draw_fcn, processes=processes)
 
 # Add negative binomial observation noise
 out_with = add_negative_binomial_noise(out_with, alpha)

tutorials/SIR/workflow_tutorial.py

@@ -28,13 +28,16 @@
 ## Generate a synthetic dataset with overdispersion ##
 ######################################################
 
-dates = pd.date_range('2022-12-01','2023-01-21')
-x = np.linspace(start=0, stop=len(dates)-1, num=len(dates))
+# Parameters 
 alpha = 0.03
-td = 10
-y = np.random.negative_binomial(1/alpha, (1/alpha)/(np.exp(x*np.log(2)/td) + (1/alpha)))
+t_d = 10
+# Sample data
+dates = pd.date_range('2022-12-01','2023-01-21')
+t = np.linspace(start=0, stop=len(dates)-1, num=len(dates))
+y = np.random.negative_binomial(1/alpha, (1/alpha)/(np.exp(t*np.log(2)/t_d) + (1/alpha)))
+# Place in a pd.Series
 d = pd.Series(index=dates, data=y, name='CASES')
-# Index name must be data for calibration to work
+# Index name must be date for calibration to work
 d.index.name = 'date'
 # Data collection only on weekdays
 d = d[d.index.dayofweek < 5]
@@ -144,7 +147,7 @@ def integrate(t, S, I, R, beta, gamma):
     fig_path = 'sampler_output/'
     identifier = 'username'
     # Perturbate previously obtained estimate
-    ndim, nwalkers, pos = perturbate_theta(theta, pert=0.10*np.ones(len(theta)), multiplier=multiplier_mcmc, bounds=bounds)
+    ndim, nwalkers, pos = perturbate_theta(theta, pert=[0.10,], multiplier=multiplier_mcmc, bounds=bounds)
     # Write some usefull settings to a pickle file (no pd.Timestamps or np.arrays allowed!)
     settings={'start_calibration': start_date.strftime("%Y-%m-%d"), 'end_calibration': end_date.strftime("%Y-%m-%d"),
               'n_chains': nwalkers, 'starting_estimate': list(theta)}