From 782e24cd80a766e378da647460048e4b5a5c591f Mon Sep 17 00:00:00 2001 From: Tijs Alleman Date: Mon, 11 Dec 2023 13:18:44 +0100 Subject: [PATCH 1/6] first error --- docs/quickstart.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/quickstart.md b/docs/quickstart.md index e3bf752..c90f232 100644 --- a/docs/quickstart.md +++ b/docs/quickstart.md @@ -174,9 +174,9 @@ 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) print(out) From 3e181bb7d03e1449b8172a073aa04a04ae5adc47 Mon Sep 17 00:00:00 2001 From: Tijs Alleman Date: Mon, 11 Dec 2023 13:20:31 +0100 Subject: [PATCH 2/6] second mistake in docs --- docs/quickstart.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/quickstart.md b/docs/quickstart.md index c90f232..0424e7f 100644 --- a/docs/quickstart.md +++ b/docs/quickstart.md @@ -178,7 +178,7 @@ coordinates={'age_group': ['0-5','5-15', '15-65','65-120']} # Initialize model model = ODE_SIR_SI(states=init_states, parameters=params, coordinates=coordinates) # Simulate the model -model.sim(120) +out = model.sim(120) print(out) ``` From 58ec8b2f7450576a11a93487901ffc39cec5b8e2 Mon Sep 17 00:00:00 2001 From: Tijs Alleman Date: Mon, 11 Dec 2023 13:27:14 +0100 Subject: [PATCH 3/6] bugs in workflow tutorial --- docs/workflow.md | 14 ++++++++------ tutorials/SIR/workflow_tutorial.py | 13 ++++++++----- 2 files changed, 16 insertions(+), 11 deletions(-) diff --git a/docs/workflow.md b/docs/workflow.md index e54ce3e..6549f15 100644 --- a/docs/workflow.md +++ b/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] ``` diff --git a/tutorials/SIR/workflow_tutorial.py b/tutorials/SIR/workflow_tutorial.py index 1d64c22..65b9b58 100644 --- a/tutorials/SIR/workflow_tutorial.py +++ b/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] From b7829cfc6ae243700968965f4cb483306e92a088 Mon Sep 17 00:00:00 2001 From: Tijs Alleman Date: Mon, 11 Dec 2023 13:31:16 +0100 Subject: [PATCH 4/6] type in name 'Quasi-Poison' --> 'Quasi-Poisson' --- docs/workflow.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/workflow.md b/docs/workflow.md index 6549f15..e6d7b16 100644 --- a/docs/workflow.md +++ b/docs/workflow.md @@ -167,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. From d251de8fac81bc83e43b2ba1ece62ecee72472df Mon Sep 17 00:00:00 2001 From: Tijs Alleman Date: Mon, 11 Dec 2023 13:31:29 +0100 Subject: [PATCH 5/6] 'parameter_names' --> 'parameters' --- docs/workflow.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/workflow.md b/docs/workflow.md index e6d7b16..5458a76 100644 --- a/docs/workflow.md +++ b/docs/workflow.md @@ -108,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): From ceeb855ae9245d7f1e2eb7660e62178a21e09029 Mon Sep 17 00:00:00 2001 From: Tijs Alleman Date: Mon, 11 Dec 2023 13:42:29 +0100 Subject: [PATCH 6/6] other bugs in workflow tutorial --- docs/workflow.md | 12 ++++++------ tutorials/SIR/workflow_tutorial.py | 2 +- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/docs/workflow.md b/docs/workflow.md index 5458a76..73ae15f 100644 --- a/docs/workflow.md +++ b/docs/workflow.md @@ -253,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') @@ -303,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, @@ -362,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 @@ -400,10 +400,10 @@ 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}) ``` @@ -411,7 +411,7 @@ Next, we compare the simulations with and without the use of our time-dependent ``` # 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) diff --git a/tutorials/SIR/workflow_tutorial.py b/tutorials/SIR/workflow_tutorial.py index 65b9b58..ed6717d 100644 --- a/tutorials/SIR/workflow_tutorial.py +++ b/tutorials/SIR/workflow_tutorial.py @@ -147,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)}