Getting Started¶

Installation¶

There exists three different possiblities to run the models:

Clone the repository, with the latest release:

git clone --branch v0.1.8 https://github.com/Priesemann-Group/covid19_inference

Install the module via pip

pip install git+https://github.com/Priesemann-Group/covid19_inference.git@v0.1.8

3. Run the notebooks directly in Google Colab. At the top of the notebooks files there should be a symbol which opens them directly in a Google Colab instance.

First Steps¶

To get started, we recommend to look at one of the currently two example notebooks:

SIR model with one german state
This model is similar to the one discussed in our paper: Inferring COVID-19 spreading rates and potential change points for case number forecasts. The difference is that the delay between infection and report is now lognormal distributed and not fixed.
Hierarchical model of the German states
This builds a hierarchical bayesian model of the states of Germany. Caution, seems to be currently broken!

We can for example recommend the following articles about bayesian modeling:

As a introduction to Bayesian statistics and the python package (PyMC3) that we use: https://docs.pymc.io/notebooks/api_quickstart.html

This is a good post about hierarchical Bayesian models in general: https://statmodeling.stat.columbia.edu/2014/01/21/everything-need-know-bayesian-statistics-learned-eight-schools/

Examples¶

We supply a number of examples which can be found in the scripts folder of the github repository.

These examples are given as python files and interactive ipython notebooks. The python files get automatically converted into ipython notebooks for easier use with google colab. The conversion is done by a slightly modified version of the python2jupyter module, which can be found here.

For starters the most useful examples are the non hierarchical one bundesland example and the hierarchical analysis of the bundeslaender.

Disclaimer¶

We evaluate the data provided by the John Hopkins University link. We exclude any liability with regard to the quality and accuracy of the data used, and also with regard to the correctness of the statistical analysis. The evaluation of the different growth phases represents solely our personal opinion.

The number of cases reported may be significantly lower than the number of people actually infected. Also, we must point out that week-ends and changes in the test system may lead to fluctuations in reported cases that have no equivalent in actual case numbers.

Certainly, at this stage all statistical predictions are subject to great uncertainty because the general trends of the epidemic are not yet clear. In any case, the statistical trends that we interpret from the data are only suitable for predictions if the measures taken by the government and authorities to contain the pandemic remain in force and are being followed by the population. We must also point out that, even if the statistics indicate that the epidemic is under control, we may at any time see a resurgence of infection figures until the disease is eradicated worldwide.

Model ¶

If you are familiar with pymc3, then looking at the example below should explain how our model works. Otherwise, here is a quick overivew:

First, we have to create an instance of the base class (that is derived from pymc3s model class). It has some convenient properties to get the range of the data, simulation length and so forth.
We then add details that base model. They correspond to the actual (physical) model features, such as the change points, the reporting delay and the week modulation.
- Every feature has it’s own function that takes in arguments to set prior assumptions.
- Sometimes they also take in input (data, reported cases … ) but none of the function performs any actual modifactions on the data. They only tell pymc3 what it is supposed to do during the sampling.
- None of our functions actually modifies any data. They rather define ways how pymc3 should modify data during the sampling.
- Most of the feature functions add variables to the pymc3.trace, see the function arguments that start with name_.
in pymc3 it is common to use a context, as we also do in the example. everything within the block with cov19.model.Cov19Model(**params_model) as this_model: automagically applies to this_model. Alternatively, you could provide a keyword to each function model=this_model.

Example ¶

import datetime

import pymc3 as pm
import numpy as np
import covid19_inference as cov19

# limit the data range
bd = datetime.datetime(2020, 3, 2)

# download data
jhu = cov19.data_retrieval.JHU(auto_download=True)
new_cases = jhu.get_new(country="Germany", data_begin=bd)

# set model parameters
params_model = dict(
    new_cases_obs=new_cases,
    data_begin=bd,
    fcast_len=28,
    diff_data_sim=16,
    N_population=83e6,
)

# change points like in the paper
change_points = [
    dict(pr_mean_date_transient=datetime.datetime(2020, 3, 9)),
    dict(pr_mean_date_transient=datetime.datetime(2020, 3, 16)),
    dict(pr_mean_date_transient=datetime.datetime(2020, 3, 23)),
]

# create model instance and add details
with cov19.model.Cov19Model(**params_model) as this_model:
    # apply change points, lambda is in log scale
    lambda_t_log = cov19.model.lambda_t_with_sigmoids(
        pr_median_lambda_0=0.4,
        pr_sigma_lambda_0=0.5,
        change_points_list=change_points,
    )

    # prior for the recovery rate
    mu = pm.Lognormal(name="mu", mu=np.log(1 / 8), sigma=0.2)

    # new Infected day over day are determined from the SIR model
    new_I_t = cov19.model.SIR(lambda_t_log, mu)

    # model the reporting delay, our prior is ten days
    new_cases_inferred_raw = cov19.model.delay_cases(
        cases=new_I_t,
        pr_mean_of_median=10,
    )

    # apply a weekly modulation, fewer reports during weekends
    new_cases_inferred = cov19.model.week_modulation(new_cases_inferred_raw)

    # set the likeliehood
    cov19.model.student_t_likelihood(new_cases_inferred)

# run the sampling
trace = pm.sample(model=this_model, tune=50, draws=10, init="advi+adapt_diag")

Table of Contents

Model

Model Base Class ¶

class covid19_inference.model.Cov19Model(new_cases_obs, data_begin, fcast_len, diff_data_sim, N_population, name='', model=None, shifted_cases=True)[source]¶

Abstract base class for the dynamic model of covid-19 propagation. Derived from pymc3.Model.

Parameters below are passed to the constructor.

Attributes (Variables) are available after creation and can be accessed from every instance. Some background:

The simulation starts diff_data_sim days before the data.

The data has a certain length, on which the inference is based. This length is given by new_cases_obs.

After the inference, a forecast takes of length fcast_len takes place, starting on the day after the last data point in new_cases_obs.

In total, traces produced by a model run have the length sim_len = diff_data_sim + data_len + fcast_len

Date ranges include both boundaries. For example, if data_begin is March 1 and data_end is March 3 then data_len will be 3.

Parameters

new_cases_obs (1 or 2d array) – If the array is two-dimensional, an hierarchical model will be constructed. First dimension is then time, the second the region/country.
data_begin (datatime.datetime) – Date of the first data point
fcast_len (int) – Number of days the simulations runs longer than the data
diff_data_sim (int) – Number of days the simulation starts earlier than the data. Should be significantly longer than the delay between infection and report of cases.
N_population (number or 1d array) – Number of inhabitance in region, needed for the S(E)IR model. Is ideally 1 dimensional if new_cases_obs is 2 dimensional
name (string) – suffix appended to the name of random variables saved in the trace
model – specify a model, if this one should expand another
shifted_cases (bool) – when enabled (True), interprets short intervals of zero cases as days, where no reporting happens and adds model cases to next non-zero-case day

Variables

new_cases_obs (1 or 2d array) – as passed during construction
data_begin (datatime.datetime) – date of the first data point in the data
data_end (datatime.datetime) – date of the last data point in the data
sim_begin (datatime.datetime) – date at which the simulation begins
sim_end (datatime.datetime) – date at which the simulation ends (should match fcast_end)
fcast_begin (datatime.datetime) – date at which the forecast starts (should be one day after data_end)
fcast_end (datatime.datetime) – data at which the forecast ends
data_len (int) – total number of days in the data
sim_len (int) – total number of days in the simulation
fcast_len (int) – total number of days in the forecast
diff_data_sim (int) – difference in days between the simulation begin and the data begin. The simulation starting time is usually earlier than the data begin.

Example

with Cov19Model(**params) as model:
    # Define model here

property untransformed_freeRVs¶

Returns the names of all free parameters of the model, usefull for plotting!

Returns: list – all variable names

Compartmental models ¶

SIR — susceptible-infected-recovered¶

covid19_inference.model.SIR(lambda_t_log, mu, name_new_I_t='new_I_t', name_I_begin='I_begin', name_I_t='I_t', name_S_t='S_t', pr_I_begin=100, model=None, return_all=False)[source]¶

Implements the susceptible-infected-recovered model.

Parameters

lambda_t_log (TensorVariable) – time series of the logarithm of the spreading rate, 1 or 2-dimensional. If 2-dimensional the first dimension is time.
mu (TensorVariable) – the recovery rate $\mu$ , typically a random variable. Can be 0 or 1-dimensional. If 1-dimensional, the dimension are the different regions.
name_new_I_t (str, optional) – Name of the new_I_t variable
name_I_begin (str, optional) – Name of the I_begin variable
name_I_t (str, optional) – Name of the I_t variable, set to None to avoid adding as trace variable.
name_S_t (str, optional) – Name of the S_t variable, set to None to avoid adding as trace variable.
pr_I_begin (float or array_like or Variable) – Prior beta of the Half-Cauchy distribution of $I(0)$ . if type is tt.Constant, I_begin will not be inferred by pymc3
model (Cov19Model) – if none, it is retrieved from the context
return_all (bool) – if True, returns name_new_I_t, name_I_t, name_S_t otherwise returns only name_new_I_t

Returns

new_I_t (TensorVariable) – time series of the number daily newly infected persons.
I_t (TensorVariable) – time series of the infected (if return_all set to True)
S_t (TensorVariable) – time series of the susceptible (if return_all set to True)

More Details¶

$I_{new}(t) &= \lambda_t I(t-1) \frac{S(t-1)}{N} \\ S(t) &= S(t-1) - I_{new}(t) \\ I(t) &= I(t-1) + I_{new}(t) - \mu I(t)$

The prior distributions of the recovery rate $\mu$ and $I(0)$ are set to

$\mu &\sim \text{LogNormal}\left[ \log(\text{pr\_median\_mu}), \text{pr\_sigma\_mu} \right] \\ I(0) &\sim \text{HalfCauchy}\left[ \text{pr\_beta\_I\_begin} \right]$

SEIR-like — susceptible-exposed-infected-recovered¶

covid19_inference.model.SEIR(lambda_t_log, mu, name_new_I_t='new_I_t', name_new_E_t='new_E_t', name_I_t='I_t', name_S_t='S_t', name_I_begin='I_begin', name_new_E_begin='new_E_begin', name_median_incubation='median_incubation', pr_I_begin=100, pr_new_E_begin=50, pr_median_mu=0.125, pr_mean_median_incubation=4, pr_sigma_median_incubation=1, sigma_incubation=0.4, pr_sigma_mu=0.2, model=None, return_all=False)[source]¶

Implements a model similar to the susceptible-exposed-infected-recovered model. Instead of a exponential decaying incubation period, the length of the period is lognormal distributed.

Parameters

lambda_t_log (TensorVariable) – time series of the logarithm of the spreading rate, 1 or 2-dimensional. If 2-dimensional, the first dimension is time.
mu (TensorVariable) – the recovery rate $\mu$ , typically a random variable. Can be 0 or 1-dimensional. If 1-dimensional, the dimension are the different regions.
name_new_I_t (str, optional) – Name of the new_I_t variable
name_I_t (str, optional) – Name of the I_t variable
name_S_t (str, optional) – Name of the S_t variable
name_I_begin (str, optional) – Name of the I_begin variable
name_new_E_begin (str, optional) – Name of the new_E_begin variable
name_median_incubation (str) – The name under which the median incubation time is saved in the trace
pr_I_begin (float or array_like) – Prior beta of the HalfCauchy distribution of $I(0)$ . if type is tt.Variable, I_begin will be set to the provided prior as a constant.
pr_new_E_begin (float or array_like) – Prior beta of the HalfCauchy distribution of $E(0)$ .
pr_median_mu (float or array_like) – Prior for the median of the Lognormal distribution of the recovery rate $\mu$ .
pr_mean_median_incubation – Prior mean of the Normal distribution of the median incubation delay $d_{\text{incubation}}$ . Defaults to 4 days [Nishiura2020], which is the median serial interval (the important measure here is not exactly the incubation period, but the delay until a person becomes infectious which seems to be about 1 day earlier as showing symptoms).
pr_sigma_median_incubation (number or None) – Prior sigma of the Normal distribution of the median incubation delay $d_{\text{incubation}}$ . If None, the incubation time will be fixed to the value of pr_mean_median_incubation instead of a random variable Default is 1 day.
sigma_incubation – Scale parameter of the Lognormal distribution of the incubation time/ delay until infectiousness. The default is set to 0.4, which is about the scale found in [Nishiura2020], [Lauer2020].
pr_sigma_mu (float or array_like) – Prior for the sigma of the lognormal distribution of recovery rate $\mu$ .
model (Cov19Model) – if none, it is retrieved from the context
return_all (bool) – if True, returns name_new_I_t, name_new_E_t, name_I_t, name_S_t otherwise returns only name_new_I_t

Returns

name_new_I_t (TensorVariable) – time series of the number daily newly infected persons.
name_new_E_t (TensorVariable) – time series of the number daily newly exposed persons. (if return_all set to True)
name_I_t (TensorVariable) – time series of the infected (if return_all set to True)
name_S_t (TensorVariable) – time series of the susceptible (if return_all set to True)

More Details¶

$E_{\text{new}}(t) &= \lambda_t I(t-1) \frac{S(t)}{N} \\ S(t) &= S(t-1) - E_{\text{new}}(t) \\ I_\text{new}(t) &= \sum_{k=1}^{10} \beta(k) E_{\text{new}}(t-k) \\ I(t) &= I(t-1) + I_{\text{new}}(t) - \mu I(t) \\ \beta(k) & = P(k) \sim \text{LogNormal}\left[ \log(d_{\text{incubation}}), \text{sigma\_incubation} \right]$

The recovery rate $\mu$ and the incubation period is the same for all regions and follow respectively:

$P(\mu) &\sim \text{LogNormal}\left[ \text{log(pr\_median\_mu), pr\_sigma\_mu} \right] \\ P(d_{\text{incubation}}) &\sim \text{Normal}\left[ \text{pr\_mean\_median\_incubation, pr\_sigma\_median\_incubation} \right]$

The initial number of infected and newly exposed differ for each region and follow prior HalfCauchy distributions:

$E(t) &\sim \text{HalfCauchy}\left[ \text{pr\_beta\_E\_begin} \right] \:\: \text{for} \: t \in {-9, -8, ..., 0}\\ I(0) &\sim \text{HalfCauchy}\left[ \text{pr\_beta\_I\_begin} \right].$

References¶

Nishiura2020(1,2)

Nishiura, H.; Linton, N. M.; Akhmetzhanov, A. R. Serial Interval of Novel Coronavirus (COVID-19) Infections. Int. J. Infect. Dis. 2020, 93, 284–286. https://doi.org/10.1016/j.ijid.2020.02.060.
Lauer2020

Lauer, S. A.; Grantz, K. H.; Bi, Q.; Jones, F. K.; Zheng, Q.; Meredith, H. R.; Azman, A. S.; Reich, N. G.; Lessler, J. The Incubation Period of Coronavirus Disease 2019 (COVID-19) From Publicly Reported Confirmed Cases: Estimation and Application. Ann Intern Med 2020. https://doi.org/10.7326/M20-0504.

Likelihood ¶

covid19_inference.model.student_t_likelihood(cases, name_student_t='_new_cases_studentT', name_sigma_obs='sigma_obs', pr_beta_sigma_obs=30, nu=4, offset_sigma=1, model=None, data_obs=None)[source]¶

Set the likelihood to apply to the model observations (model.new_cases_obs) We assume a StudentT distribution because it is robust against outliers [Lange1989]. The likelihood follows:

$P(\text{data\_obs}) &\sim StudentT(\text{mu} = \text{new\_cases\_inferred}, sigma =\sigma, \text{nu} = \text{nu})\\ \sigma &= \sigma_r \sqrt{\text{new\_cases\_inferred} + \text{offset\_sigma}}$

The parameter $\sigma_r$ follows a HalfCauchy prior distribution with parameter beta set by pr_beta_sigma_obs. If the input is 2 dimensional, the parameter $\sigma_r$ is different for every region.

Parameters

cases (TensorVariable) – The daily new cases estimated by the model. Will be compared to the real world data data_obs. One or two dimensonal array. If 2 dimensional, the first dimension is time and the second are the regions/countries
name_student_t – The name under which the studentT distribution is saved in the trace.
name_sigma_obs – The name under which the distribution of the observable error is saved in the trace
pr_beta_sigma_obs (float) – The beta of the HalfCauchy prior distribution of $\sigma_r$ .
nu (float) – How flat the tail of the distribution is. Larger nu should make the model more robust to outliers. Defaults to 4 [Lange1989].
offset_sigma (float) – An offset added to the sigma, to make the inference procedure robust. Otherwise numbers of cases would lead to very small errors and diverging likelihoods. Defaults to 1.
model – The model on which we want to add the distribution
data_obs (array) – The data that is observed. By default it is model.new_cases_obs

Returns

None

References

Lange1989(1,2): Lange, K., Roderick J. A. Little, & Jeremy M. G. Taylor. (1989). Robust Statistical Modeling Using the t Distribution. Journal of the American Statistical Association, 84(408), 881-896. doi:10.2307/2290063

Spreading Rate ¶

covid19_inference.model.lambda_t_with_sigmoids(change_points_list, pr_median_lambda_0, pr_sigma_lambda_0=0.5, model=None, name_lambda_t='lambda_t', hierarchical=None, sigma_lambda_cp=None, sigma_lambda_week_cp=None, prefix_lambdas='', shape=None)[source]¶

Builds a time dependent spreading rate $\lambda_t$ with change points. The change points are marked by a transient with a sigmoidal shape, with at

Parameters

change_points_list –
pr_median_lambda_0 –
pr_sigma_lambda_0 –
model (Cov19Model) – if none, it is retrieved from the context

Returns

lambda_t_log

Todo

Documentation on this

Delay ¶

covid19_inference.model.delay_cases(cases, name_delay='delay', name_cases=None, name_width='delay-width', pr_mean_of_median=10, pr_sigma_of_median=0.2, pr_median_of_width=0.3, pr_sigma_of_width=None, model=None, len_input_arr=None, len_output_arr=None, diff_input_output=None, num_variants=None)[source]¶

Convolves the input by a lognormal distribution, in order to model a delay:

We have a kernel (a distribution) of delays, one realization of this kernel is applied to each pymc3 sample.
The kernel has a median delay D and a width that correspond to this one sample. Doing the ensemble average over all samples and the respective kernels, we get two distributions: one of the median delay D and one of the width.
The (normal) distribution of the median of D is specified using pr_mean_of_median and pr_sigma_of_median.
The (lognormal) distribution of the width of the kernel of D is specified using pr_median_of_width and pr_sigma_of_width. If pr_sigma_of_width is None, the width is fixed (skipping the second distribution).

Parameters

cases (TensorVariable) – The input, typically the number of newly infected cases from the output of SIR() or SEIR().
name_delay (str) – The name under which the delay is saved in the trace, suffixes and prefixes are added depending on which variable is saved. Default : “delay”
name_cases (str or None) – The name under which the delayed cases are saved in the trace. If None, no variable will be added to the trace. Default: “delayed_cases”
pr_mean_of_median (float) – The mean of the normal distribution which models the prior median of the LogNormal delay kernel. Default: 10.0 (days)
pr_sigma_of_median (float) – The standart devaiation of normal distribution which models the prior median of the LogNormal delay kernel. Default: 0.2
pr_median_of_width (float) – The scale (width) of the LogNormal delay kernel. Default: 0.3
pr_sigma_of_width (float or None) – Whether to put a prior distribution on the scale (width) of the distribution of the delays, too. If a number is provided, the scale of the delay kernel follows a prior LogNormal distribution, with median pr_median_scale_delay and scale pr_sigma_scale_delay. Default: None, and no distribution is applied.
model (Cov19Model or None) – The model to use. Default: None, model is retrieved automatically from the context

Other Parameters

len_input_arr – Length of new_I_t. By default equal to model.sim_len. Necessary because the shape of theano tensors are not defined at when the graph is built.
len_output_arr (int) – Length of the array returned. By default it set to the length of the cases_obs saved in the model plus the number of days of the forecast.
diff_input_output (int) – Number of days the returned array begins later then the input. Should be significantly larger than the median delay. By default it is set to the model.diff_data_sim.
num_variants (int) – If you are not using the hierachical model but still want to apply a delay to multidimensional casenumbers. This is the shape of your last dimension.

Returns

delayed_cases (TensorVariable) – The delayed input $y_\text{delayed}(t)$ , typically the daily number new cases that one expects to measure.

More Details¶

$y_\text{delayed}(t) &= \sum_{\tau=0}^T y_\text{input}(\tau) \text{LogNormal}\left[ \log(\text{delay}), \text{pr\_median\_scale\_delay} \right](t - \tau) \\ \log(\text{delay}) &= \text{Normal}\left[ \log(\text{pr\_sigma\_delay} ), \text{pr\_sigma\_delay} \right]$

The LogNormal distribution is a function evaluated at $t - \tau$ .

If the model is 2-dimensional (hierarchical), the $\log(\text{delay})$ is hierarchically modelled with the hierarchical_normal() function using the default parameters except that the prior sigma of delay_L2 is HalfNormal distributed (error_cauchy=False).

Week modulation ¶

covid19_inference.model.week_modulation(cases, name_cases=None, name_weekend_factor='weekend_factor', name_offset_modulation='offset_modulation', week_modulation_type='abs_sine', pr_mean_weekend_factor=0.3, pr_sigma_weekend_factor=0.5, weekend_days=(6, 7), model=None)[source]¶

Adds a weekly modulation of the number of new cases:

$\text{new\_cases} &= \text{new\_cases\_raw} \cdot (1-f(t))\,, \qquad\text{with}\\ f(t) &= f_w \cdot \left(1 - \left|\sin\left(\frac{\pi}{7} t- \frac{1}{2}\Phi_w\right)\right| \right),$

if week_modulation_type is "abs_sine" (the default). If week_modulation_type is "step", the new cases are simply multiplied by the weekend factor on the days set by weekend_days

The weekend factor $f_w$ follows a Lognormal distribution with median pr_mean_weekend_factor and sigma pr_sigma_weekend_factor. It is hierarchically constructed if the input is two-dimensional by the function hierarchical_normal() with default arguments.

The offset from Sunday $\Phi_w$ follows a flat VonMises distribution and is the same for all regions.

Parameters

cases (TensorVariable) – The input array of daily new cases, can be one- or two-dimensional
name_cases (str or None,) – The name under which to save the cases as a trace variable. Default: None, cases are not stored in the trace.
week_modulation_type (str) – The type of modulation, accepts "step" or "abs_sine (the default).
pr_mean_weekend_factor (float) – Sets the prior mean of the factor $f_w$ by which weekends are counted.
pr_sigma_weekend_factor (float) – Sets the prior sigma of the factor $f_w$ by which weekends are counted.
weekend_days (tuple of ints) – The days counted as weekend if week_modulation_type is "step"
model (Cov19Model) – if none, it is retrieved from the context

Returns

new_cases (TensorVariable)

Utility ¶

covid19_inference.model.utility.hierarchical_normal(pr_mean, pr_sigma, name_L1='delay_hc_L1', name_L2='delay_hc_L2', name_sigma='delay_hc_sigma', model=None, error_fact=2.0, error_cauchy=True, shape=None)[source]¶

Implements an hierarchical normal model:

$x_\text{L1} &= Normal(\text{pr\_mean}, \text{pr\_sigma})\\ y_{i, \text{L2}} &= Normal(x_\text{L1}, \sigma_\text{L2})\\ \sigma_\text{L2} &= HalfCauchy(\text{error\_fact} \cdot \text{pr\_sigma})$

It is however implemented in a non-centered way, that the second line is changed to:

$y_{i, \text{L2}} &= x_\text{L1} + Normal(0,1) \cdot \sigma_\text{L2}$

See for example https://arxiv.org/pdf/1312.0906.pdf

Parameters

name_L1 (str) – Name under which $x_\text{L1}$ is saved in the trace.
name_L2 (str) – Name under which $x_\text{L2}$ is saved in the trace. The non-centered distribution in addition saved with a suffix _raw added.
name_sigma (str) – Name under which $\sigma_\text{L2}$ is saved in the trace.
pr_mean (float) – Prior mean of $x_\text{L1}$
pr_sigma (float) – Prior sigma for $x_\text{L1}$ and (muliplied by error_fact) for $\sigma_\text{L2}$
len_L2 (int) – length of $y_\text{L2}$
error_fact (float) – Factor by which pr_sigma is multiplied as prior for sigma_text{L2}
error_cauchy (bool) – if False, a $HalfNormal$ distribution is used for $\sigma_\text{L2}$ instead of $HalfCauchy$

Returns

y (TensorVariable) – the random variable $y_\text{L2}$
x (TensorVariable) – the random variable $x_\text{L1}$

covid19_inference.model.utility.tt_lognormal(x, mu, sigma)[source]¶: Calculates a lognormal pdf for integer spaced x input.

covid19_inference.model.utility.tt_gamma(x, mu=None, sigma=None, alpha=None, beta=None)[source]¶: Calculates a gamma distribution pdf for integer spaced x input. Parametrized similarly to Gamma

Data Retrieval ¶

Table of Contents

Data Retrieval

Utility ¶

covid19_inference.data_retrieval.retrieval.set_data_dir(fname=None, permissions=None)[source]¶: Set the global variable _data_dir. New downloaded data is placed there. If no argument provided we try the default tmp directory. If permissions are not provided, uses defaults if fname is in user folder. If not in user folder, tries to set 777.

covid19_inference.data_retrieval.retrieval.backup_instances(trace=None, model=None, fname='latest_')[source]¶: helper to save or load trace and model instances. loads from fname if provided traces and model variables are None, else saves them there.

Johns Hops University ¶

class covid19_inference.data_retrieval.JHU(auto_download=False)[source]¶

This class can be used to retrieve and filter the dataset from the online repository of the coronavirus visual dashboard operated by the Johns Hopkins University.

Features

download all files from the online repository of the coronavirus visual dashboard operated by the Johns Hopkins University.
filter by deaths, confirmed cases and recovered cases
filter by country and state
filter by date

Example

jhu = cov19.data_retrieval.JHU()
jhu.download_all_available_data()

#Acess the data by
jhu.data
#or
jhu.get_new("confirmed","Italy")
jhu.get_total(filter)

__init__(auto_download=False)[source]¶

On init of this class the base Retrieval Class __init__ is called, with jhu specific arguments.

Parameters: auto_download (bool, optional) – Whether or not to automatically call the download_all_available_data() method. One should explicitly call this method for more configuration options (default: false)

download_all_available_data(force_local=False, force_download=False)[source]¶

Attempts to download from the main urls (self.url_csv) which was set on initialization of this class. If this fails it downloads from the fallbacks. It can also be specified to use the local files or to force the download. The download methods get inhereted from the base retrieval class.

Parameters

force_local (bool, optional) – If True forces to load the local files.
force_download (bool, optional) – If True forces the download of new files

get_total_confirmed_deaths_recovered(country: str = None, state: str = None, begin_date: datetime.datetime = None, end_date: datetime.datetime = None)[source]¶

Retrieves all confirmed, deaths and recovered cases from the Johns Hopkins University dataset as a DataFrame with datetime index. Can be filtered by country and state, if only a country is given all available states get summed up.

Parameters

country (str, optional) – name of the country (the “Country/Region” column), can be None if the whole summed up data is wanted (why would you do this?)
state (str, optional) – name of the state (the “Province/State” column), can be None if country is set or the whole summed up data is wanted
begin_date (datetime.datetime, optional) – intial date for the returned data, if no value is given the first date in the dataset is used
end_date (datetime.datetime, optional) – last date for the returned data, if no value is given the most recent date in the dataset is used

Returns

pandas.DataFrame

get_new(value='confirmed', country: str = None, state: str = None, data_begin: datetime.datetime = None, data_end: datetime.datetime = None)[source]¶

Retrieves all new cases from the Johns Hopkins University dataset as a DataFrame with datetime index. Can be filtered by value, country and state, if only a country is given all available states get summed up.

Parameters

value (str) – Which data to return, possible values are - “confirmed”, - “recovered”, - “deaths” (default: “confirmed”)
country (str, optional) – name of the country (the “Country/Region” column), can be None
state (str, optional) – name of the state (the “Province/State” column), can be None
begin_date (datetime.datetime, optional) – intial date for the returned data, if no value is given the first date in the dataset is used
end_date (datetime.datetime, optional) – last date for the returned data, if no value is given the most recent date in the dataset is used

Returns

pandas.DataFrame – table with new cases and the date as index

get_total(value='confirmed', country: str = None, state: str = None, data_begin: datetime.datetime = None, data_end: datetime.datetime = None)[source]¶

Retrieves all total/cumulative cases from the Johns Hopkins University dataset as a DataFrame with datetime index. Can be filtered by value, country and state, if only a country is given all available states get summed up.

Parameters

value (str) – Which data to return, possible values are - “confirmed”, - “recovered”, - “deaths” (default: “confirmed”)
country (str, optional) – name of the country (the “Country/Region” column), can be None
state (str, optional) – name of the state (the “Province/State” column), can be None
begin_date (datetime.datetime, optional) – intial date for the returned data, if no value is given the first date in the dataset is used
end_date (datetime.datetime, optional) – last date for the returned data, if no value is given the most recent date in the dataset is used

Returns

pandas.DataFrame – table with total/cumulative cases and the date as index

filter_date(df, begin_date: datetime.datetime = None, end_date: datetime.datetime = None)[source]¶

Returns give dataframe between begin and end date. Dataframe has to have a datetime index.

Parameters

begin_date (datetime.datetime, optional) – First day that should be filtered
end_date (datetime.datetime, optional) – Last day that should be filtered

Returns

pandas.DataFrame

get_possible_countries_states()[source]¶

Can be used to get a list with all possible states and coutries.

Returns: pandas.DataFrame in the format

Robert Koch Institute ¶

class covid19_inference.data_retrieval.RKI(auto_download=False)[source]¶

This class can be used to retrieve and filter the dataset from the Robert Koch Institute Robert Koch Institute. The data gets retrieved from the arcgis dashboard.

Features

download the full dataset
filter by date
filter by bundesland
filter by recovered, deaths and confirmed cases

Example

rki = cov19.data_retrieval.RKI()
rki.download_all_available_data()

#Acess the data by
rki.data
#or
rki.get_new("confirmed","Sachsen")
rki.get_total(filter)

__init__(auto_download=False)[source]¶

On init of this class the base Retrieval Class __init__ is called, with rki specific arguments.

Parameters: auto_download (bool, optional) – Whether or not to automatically call the download_all_available_data() method. One should explicitly call this method for more configuration options (default: false)

download_all_available_data(force_local=False, force_download=False)[source]¶

Attempts to download from the main url (self.url_csv) which was given on initialization. If this fails download from the fallbacks. It can also be specified to use the local files or to force the download. The download methods get inhereted from the base retrieval class.

Parameters

force_local (bool, optional) – If True forces to load the local files.
force_download (bool, optional) – If True forces the download of new files

get_total(value='confirmed', bundesland: str = None, landkreis: str = None, data_begin: datetime.datetime = None, data_end: datetime.datetime = None, date_type: str = 'date', age_group=None)[source]¶

Gets all total confirmed cases for a region as dataframe with date index. Can be filtered with multiple arguments.

Parameters

value (str) – Which data to return, possible values are - “confirmed”, - “recovered”, - “deaths” (default: “confirmed”)
bundesland (str, optional) – if no value is provided it will use the full summed up dataset for Germany
landkreis (str, optional) – if no value is provided it will use the full summed up dataset for the region (bundesland)
data_begin (datetime.datetime, optional) – initial date, if no value is provided it will use the first possible date
data_end (datetime.datetime, optional) – last date, if no value is provided it will use the most recent possible date
date_type (str, optional) – type of date to use: reported date ‘date’ (Meldedatum in the original dataset), or symptom date ‘date_ref’ (Refdatum in the original dataset)
age_group (str, optional) – Choosen age group. To get the possible combinations use possible_age_groups().

Returns

pandas.DataFrame

get_new(value='confirmed', bundesland: str = None, landkreis: str = None, data_begin: datetime.datetime = None, data_end: datetime.datetime = None, date_type: str = 'date', age_group=None)[source]¶

Retrieves all new cases from the Robert Koch Institute dataset as a DataFrame with datetime index. Can be filtered by value, bundesland and landkreis, if only a country is given all available states get summed up.

Parameters

value (str) – Which data to return, possible values are - “confirmed”, - “recovered”, - “deaths” (default: “confirmed”)
bundesland (str, optional) – if no value is provided it will use the full summed up dataset for Germany
landkreis (str, optional) – if no value is provided it will use the full summed up dataset for the region (bundesland)
data_begin (datetime.datetime, optional) – intial date for the returned data, if no value is given the first date in the dataset is used, if none is given could yield errors
data_end (datetime.datetime, optional) – last date for the returned data, if no value is given the most recent date in the dataset is used
age_group (str, optional) – Choosen age group. To get the possible combinations use possible_age_groups().

Returns

pandas.DataFrame – table with daily new confirmed and the date as index

filter(data_begin: datetime.datetime = None, data_end: datetime.datetime = None, variable='confirmed', date_type='date', level=None, value=None, age_group=None)[source]¶

Filters the obtained dataset for a given time period and returns an array ONLY containing only the desired variable.

Parameters

data_begin (datetime.datetime, optional) – initial date, if no value is provided it will use the first possible date
data_end (datetime.datetime, optional) – last date, if no value is provided it will use the most recent possible date
variable (str, optional) – type of variable to return possible types are: “confirmed” : cases (default) “AnzahlTodesfall” : deaths “AnzahlGenesen” : recovered
date_type (str, optional) – type of date to use: reported date ‘date’ (Meldedatum in the original dataset), or symptom date ‘date_ref’ (Refdatum in the original dataset)
level (str, optional) –

possible strings are:
”None” : return data from all Germany (default) “Bundesland” : a state “Landkreis” : a region
value (str, optional) – string of the state/region e.g. “Sachsen”
age_group (str, optional) – Choosen age group. To get the possible combinations use possible_age_groups().

Returns

pd.DataFrame – array with ONLY the requested variable, in the requested range. (one dimensional)

filter_all_bundesland(begin_date: datetime.datetime = None, end_date: datetime.datetime = None, variable='confirmed', date_type='date')[source]¶

Filters the full RKI dataset

Parameters

df (DataFrame) – RKI dataframe, from get_rki()
begin_date (datetime.datetime) – initial date to return
end_date (datetime.datetime) – last date to return
variable (str, optional) – type of variable to return: cases (“AnzahlFall”), deaths (“AnzahlTodesfall”), recovered (“AnzahlGenesen”)
date_type (str, optional) – type of date to use: reported date ‘date’ (Meldedatum in the original dataset), or symptom date ‘date_ref’ (Refdatum in the original dataset)

Returns

pd.DataFrame – DataFrame with datetime dates as index, and all German regions (bundesländer) as columns

possible_age_groups()[source]¶: Returns the valid age groups in the dataset.

Robert Koch Institute situation reports ¶

class covid19_inference.data_retrieval.RKIsituationreports(auto_download=False)[source]¶

As mentioned by Matthias Linden, the daily situation reports have more available data. This class retrieves this additional data from Matthias website and parses it into the format we use i.e. a datetime index.

Interesting new data is for example ICU cases, deaths and recorded symptoms. For now one can look at the data by running

Example

rki_si_re = cov19.data_retrieval.RKIsituationreports(True)
print(rki_si_re.data)

Todo

Filter functions for ICU, Symptoms and maybe even daily new cases for the respective categories.

__init__(auto_download=False)[source]¶

On init of this class the base Retrieval Class __init__ is called, with rki situation reports specific arguments.

Parameters: auto_download (bool, optional) – Whether or not to automatically call the download_all_available_data() method. One should explicitly call this method for more configuration options (default: false)

download_all_available_data(force_local=False, force_download=False)[source]¶

Attempts to download from the main url (self.url_csv) which was given on initialization. If this fails download from the fallbacks. It can also be specified to use the local files or to force the download. The download methods get inhereted from the base retrieval class.

Parameters

force_local (bool, optional) – If True forces to load the local files.
force_download (bool, optional) – If True forces the download of new files

Google ¶

class covid19_inference.data_retrieval.GOOGLE(auto_download=False)[source]¶

This class can be used to retrieve the mobility dataset from Google.

Example

gl = cov19.data_retrieval.GOOGLE()
gl.download_all_available_data()

#Acess the data by
gl.data
#or
gl.get_changes(filter)

__init__(auto_download=False)[source]¶

On init of this class the base Retrieval Class __init__ is called, with google specific arguments.

Parameters: auto_download (bool, optional) – Whether or not to automatically call the download_all_available_data() method. One should explicitly call this method for more configuration options (default: false)

download_all_available_data(force_local=False, force_download=False)[source]¶

Attempts to download from the main url (self.url_csv) which was given on initialization. If this fails download from the fallbacks. It can also be specified to use the local files or to force the download. The download methods get inhereted from the base retrieval class.

Parameters

force_local (bool, optional) – If True forces to load the local files.
force_download (bool, optional) – If True forces the download of new files

get_changes(country: str, state: str = None, region: str = None, data_begin: datetime.datetime = None, data_end: datetime.datetime = None)[source]¶

Returns a dataframe with the relative changes in mobility to a baseline, provided by google. They are separated into “retail and recreation”, “grocery and pharmacy”, “parks”, “transit”, “workplaces” and “residental”. Filterable for country, state and region and date.

Parameters

country (str) – Selected country for the mobility data.
state (str, optional) – State for the selected data if no value is selected the whole country is chosen
region (str, optional) – Region for the selected data if no value is selected the whole region/country is chosen
data_end (data_begin,) – Filter for the desired time period

Returns

pandas.DataFrame

get_possible_counties_states_regions()[source]¶

Can be used to obtain all different possible countries with there corresponding possible states and regions.

Returns: pandas.DataFrame

Our World in Data ¶

class covid19_inference.data_retrieval.OWD(auto_download=False)[source]¶

This class can be used to retrieve the testings dataset from Our World in Data.

Example

owd = cov19.data_retrieval.OWD()
owd.download_all_available_data()

__init__(auto_download=False)[source]¶

On init of this class the base Retrieval Class __init__ is called, with google specific arguments.

Parameters: auto_download (bool, optional) – Whether or not to automatically call the download_all_available_data() method. One should explicitly call this method for more configuration options (default: false)

download_all_available_data(force_local=False, force_download=False)[source]¶

Attempts to download from the main url (self.url_csv) which was given on initialization. If this fails download from the fallbacks. It can also be specified to use the local files or to force the download. The download methods get inhereted from the base retrieval class.

Parameters

force_local (bool, optional) – If True forces to load the local files.
force_download (bool, optional) – If True forces the download of new files

get_possible_countries()[source]¶

Can be used to obtain all different possible countries in the dataset.

Returns: pandas.DataFrame

get_total(value='tests', country=None, data_begin=None, data_end=None)[source]¶

Retrieves all new cases from the Our World in Data dataset as a DataFrame with datetime index. Can be filtered by value, country and state, if only a country is given all available states get summed up.

Parameters

value (str) – Which data to return, possible values are - “confirmed”, - “tests”, - “deaths”, - “vacination” (default: “confirmed”)
country (str) – name of the country
begin_date (datetime.datetime, optional) – intial date for the returned data, if no value is given the first date in the dataset is used
end_date (datetime.datetime, optional) – last date for the returned data, if no value is given the most recent date in the dataset is used

Returns

pandas.DataFrame – table with new cases and the date as index

get_new(value='tests', country=None, data_begin=None, data_end=None)[source]¶

Retrieves all new cases from the Our World in Data dataset as a DataFrame with datetime index. casesn be filtered by value, country and state, if only a country is given all available states get summed up.

Parameters

value (str) – Which data to return, possible values are - “confirmed”, - “tests”, - “deaths” (default: “confirmed”)
country (str) – name of the country
begin_date (datetime.datetime, optional) – intial date for the returned data, if no value is given the first date in the dataset is used
end_date (datetime.datetime, optional) – last date for the returned data, if no value is given the most recent date in the dataset is used

Returns

pandas.DataFrame – table with new cases and the date as index

Financial times ¶

class covid19_inference.data_retrieval.FINANCIAL_TIMES(auto_download=False)[source]¶

This class can be used to retrieve the excess mortality data from the Financial Times github repository.

Example

ft = cov19.data_retrieval.FINANCIAL_TIMES()
ft.download_all_available_data()

#Access the data by
ft.data
#or
ft.get(filter) #see below

__init__(auto_download=False)[source]¶

On init of this class the base Retrieval Class __init__ is called, with financial times specific arguments.

Parameters: auto_download (bool, optional) – Whether or not to automatically call the download_all_available_data() method. One should explicitly call this method for more configuration options (default: false)

download_all_available_data(force_local=False, force_download=False)[source]¶

Attempts to download from the main url (self.url_csv) which was given on initialization. If this fails download from the fallbacks. It can also be specified to use the local files or to force the download. The download methods get inhereted from the base retrieval class.

Parameters

force_local (bool, optional) – If True forces to load the local files.
force_download (bool, optional) – If True forces the download of new files

get(value='excess_deaths', country: str = 'Germany', state: str = None, data_begin: datetime.datetime = None, data_end: datetime.datetime = None)[source]¶

Retrieves specific data from the dataset, can be filtered by date, country and state.

Parameters

value (str, optional) – Which data to return, possible values are - “deaths”, - “expected_deaths”, - “excess_deaths”, - “excess_deaths_pct” (default: “excess_deaths”)
country (str, optional) –
state (str, optional) – Possible countries and states can be retrieved by the get_possible_countries_states() method.
begin_date (datetime.datetime, optional) – First day that should be filtered
end_date (datetime.datetime, optional) – Last day that should be filtered

get_possible_countries_states()[source]¶

Can be used to obtain all different possible countries with there corresponding possible states and regions.

Returns: pandas.DataFrame

Oxford COVID-19 Government Response Tracker ¶

class covid19_inference.data_retrieval.OxCGRT(auto_download=False)[source]¶

This class can be used to retrieve the datasset on goverment policies from the Oxford Covid-19 Government Response Tracker.

Example

gov_pol = cov19.data_retrieval.OxCGRT()
gov_pol.download_all_available_data()

__init__(auto_download=False)[source]¶

On init of this class the base Retrieval Class __init__ is called, with google specific arguments.

Parameters: auto_download (bool, optional) – Whether or not to automatically call the download_all_available_data() method. One should explicitly call this method for more configuration options (default: false)

download_all_available_data(force_local=False, force_download=False)[source]¶

Attempts to download from the main url (self.url_csv) which was given on initialization. If this fails download from the fallbacks. It can also be specified to use the local files or to force the download. The download methods get inhereted from the base retrieval class.

Parameters

force_local (bool, optional) – If True forces to load the local files.
force_download (bool, optional) – If True forces the download of new files

get_possible_countries()[source]¶

Can be used to obtain all different possible countries in the dataset.

Returns: pandas.DataFrame

get_possible_policies()[source]¶

Can be used to obtain all policies in there corresponding categories possible countries in the dataset.

Returns: dict

get_change_points(policies, country)[source]¶

Returns a list of change points, depending on the selected measure and country.

Parameters

policies (str, array of str) – The wanted policies. Can be an array of strings, use get_possible_policies() to get a dict of possible policies.
country (str) – Filter for country, use get_possible_countries() to get a list of possible ones.

Returns

array of dicts

get_time_data(policy, country, data_begin=None, data_end=None)[source]¶

Parameters

policy (str) – The wanted policy.
country (str) – Filter for country, use get_possible_countries() to get a list of possible ones.
data_begin (datetime.datetime, optional) – intial date for the returned data, if no value is given the first date in the dataset is used, if none is given could yield errors
data_end (datetime.datetime, optional) – last date for the returned data, if no value is given the most recent date in the dataset is used

Returns

Pandas dataframe with policy

Base Retrieval Class ¶

class covid19_inference.data_retrieval.retrieval.Retrieval(name, url_csv, fallbacks, update_interval=None, **kwargs)[source]¶

Each source class should inherit this base retrieval class, it streamlines alot of base functions. It manages downloads, multiple fallbacks and local backups via timestamp. At init of the parent class the Retrieval init should be called with the following arguments, these get saved as attributes.

An example for the usage can be seen in the _Google, _RKI and _JHU source files.

__init__(name, url_csv, fallbacks, update_interval=None, **kwargs)[source]¶

Parameters

name (str) – A name for the Parent class, mainly used for the local file backup.
url_csv (str) – The url to the main dataset as csv, if an empty string if supplied the fallback routines get used.
fallbacks (array) – Fallbacks can be filepaths to local or online sources or even methods defined in the parent class.
update_interval (datetime.timedelta) – If the local file is older than the update_interval it gets updated once the download all function is called.

_download_csv_from_source(filepath, **kwargs)[source]¶

Uses pandas read csv to download the csv file. The possible kwargs can be seen in the pandas documentation.

These kwargs can vary for the different parent classes and should be defined there!

filepathstr: Full path to the desired csv file

Returns: bool – True if the retrieval was a success, False if it failed

_fallback_handler()[source]¶: Recursivly iterate over all fallbacks and try to execute subroutines depending on the type of fallback.

_timestamp_local_old(force_local=False) → bool[source]¶

Get timestamp if it exists
compare with the date today
update if data is older than set intervall -> can be parent dependant

_save_to_local()[source]¶: Creates a local backup for the self.data pandas.DataFrame. And a timestamp for the source.

Sampling ¶

Table of Contents

Sampling

covid19_inference.sampling.robust_sample(model, tune, draws, tuning_chains, final_chains, return_tuning=False, args_start_points=None, tune_2nd=None, **kwargs)[source]¶

Samples the model by starting more chains than needed (tuning chains) and using only a reduced number final_chains for the final sampling. The final chains are randomly chosen (without replacement) weighted by their likelihood. :param model: The model :type model: Cov19Model :param tune: Number of tuning samples :type tune: int :param draws: Number of final samples :type draws: int :param tuning_chains: Number of tuning chains :type tuning_chains: int :param final_chains: Number of draw chains :type final_chains: int :param args_start_points: Arguments passed to get_start_points :type args_start_points: dict :param tune_2nd: If set, use different number of tuning samples for the second tuning :type tune_2nd: int :param **kwargs: Arguments passed to pm.sample

Returns

trace (trace as multitrace object)
trace_az (trace as arviz object)

covid19_inference.sampling.get_start_points(trace, trace_az, frames_start=None, SD_chain_logl=2.5)[source]¶

Returns the starting points such that the chains deviate at most SD_chain_logl standard deviations from the chain with the highest likelihood. :param trace: :type trace: multitrace object :param trace_az: :type trace_az: arviz trace object :param frames_start: Which frames to use for calculating the mean likelihood and its standard deviation.

By default it is set to the last third of the tuning samples

Parameters

SD_chain_logl (float) – The number of standard deviations. 2.5 as default

Returns

start_points – A list of starting points
logl_mean – The mean log-likelihood of the starting points

Plotting¶

We provide a lot of plotting functions which can be used to recreate our plots or create completely new visualizations. If you are familiar with matplotlib it should be no problem to use them extensively.

We provide three different types of functions here:

High level functions These can be used create figures similar to our paper Dehning et al. arXiv:2004.01105. The are neat little one liners which create a good looking plot from our model, but do not have a lot of customization options.
Low level functions These extend the normal matplotlib plotting functions and can be used to plot arbitrary data. They have a lot of customization options, it could take some time to get nicely looking plots with these functions though.
Helper functions These are mainly functions that manipulate data or retrieve data from our model. These do not have to be used most of the time and are only documented here for completeness.

If one just wants to recreate our figures with a different color. The easiest was is to change the default rc parameters.

covid19_inference.plot.get_rcparams_default()[source]¶: Get a Param (dict) of the default parameters. Here we set our default values. Assigned once to module variable rcParamsDefault on load.

covid19_inference.plot.set_rcparams(par)[source]¶

Sets the rcparameters used for plotting, provided instance of Param has to have the following keys (attributes):

Variables

locale (str) – region settings, passed to setlocale(). Default: “en_US”
date_format (str) – Format the date on the x axis of time-like data (see https://strftime.org/) example April 1 2020: “%m/%d” 04/01, “%-d. %B” 1. April Default “%b %-d”, becomes April 1
date_show_minor_ticks (bool) – whether to show the minor ticks (for every day). Default: True
rasterization_zorder (int or None) – Rasterizes plotted content below this value, set to None to keep everything a vector, Default: -1
draw_ci_95 (bool) – For timeseries plots, indicate 95% Confidence interval via fill between. Default: True
draw_ci_75 (bool) – For timeseries plots, indicate 75% Confidence interval via fill between. Default: False
draw_ci_50 (bool) – For timeseries plots, indicate 50% Confidence interval via fill between. Default: False
color_model (str) – Base color used for model plots, mpl compatible color code “C0”, “#303030” Default : “tab:green”
color_data (str) – Base color used for data Default : “tab:blue”
color_annot (str) – Color to use for annotations Default : “#646464”
color_prior (str) – Color to used for priors in distributions Default : “#708090”

Example

# Get default parameter
pars = cov.plot.get_rcparams_default()

# Change parameters
pars["locale"]="de_DE"
pars["color_data"]="tab:purple"

# Set parameters
cov.plot.set_rcparams(pars)

High level functions¶

covid19_inference.plot.timeseries_overview(model, trace, start=None, end=None, region=None, color=None, save_to=None, offset=0, annotate_constrained=True, annotate_watermark=True, axes=None, forecast_label='Forecast', forecast_heading='$\\bf Forecasts\\!:$', add_more_later=False)[source]¶

Create the time series overview similar to our paper. Dehning et al. arXiv:2004.01105 Contains $\lambda$ , new cases, and cumulative cases.

Parameters

model (Cov19Model) –
trace (trace instance) – needed for the data
offset (int) – offset that needs to be added to the (cumulative sum of) new cases at time model.data_begin to arrive at cumulative cases
start (datetime.datetime) – only used to set xrange in the end
end (datetime.datetime) – only used to set xrange in the end
color (str) – main color to use, default from rcParam
save_to (str or None) – path where to save the figures. default: None, not saving figures
annotate_constrained (bool) – show the unconstrained constrained annotation in lambda panel
annotate_watermark (bool) – show our watermark
axes (np.array of mpl axes) – provide an array of existing axes (from previously calling this function) to add more traces. Data will not be added again. Ideally call this first with add_more_later=True
forecast_label (str) – legend label for the forecast, default: “Forecast”
forecast_heading (str) – if add_more_later, how to label the forecast section. default: “$bf Forecasts!:$”,
add_more_later (bool) – set this to true if you plan to add multiple models to the plot. changes the layout (and the color of the fit to past data)

Returns

fig (mpl figure)
axes (np array of mpl axeses (insets not included))

Todo

Replace offset with an instance of data class that should yield the cumulative cases. we should not to calculations here.

Low level functions¶

covid19_inference.plot._timeseries(x, y, ax=None, what='data', draw_ci_95=None, draw_ci_75=None, draw_ci_50=None, date_format=True, **kwargs)[source]¶

low-level function to plot anything that has a date on the x-axis.

Parameters

x (array of datetime.datetime) – times for the x axis
y (array, 1d or 2d) – data to plot. if 2d, we plot the CI as fill_between (if CI enabled in rc params) if 2d, then first dim is realization and second dim is time matching x if 1d then first tim is time matching x
ax (mpl axes element, optional) – plot into an existing axes element. default: None
what (str, optional) – what type of data is provided in x. sets the style used for plotting: * data for data points * fcast for model forecast (prediction) * model for model reproduction of data (past)
date_format (bool, optional) – Automatic converting of index to dates default:True
kwargs (dict, optional) – directly passed to plotting mpl.

Returns

ax

covid19_inference.plot._distribution(model, trace, key, ax=None, color=None, draw_prior=True)[source]¶: Todo

documentation

Example¶

In this example we want to use the low level time series function to plot the new daily cases and deaths reported by the Robert Koch institute.

import datetime
import matplotlib.pyplot as plt
import covid19_inference as cov19

# Data retrieval i.e. download new data from RobertKochInstitue
rki = cov19.data_retrieval.RKI()
rki.download_all_available_data()

new_deaths = rki.get_new(
    value = "deaths",
    data_begin=datetime.datetime(2020,3,15), #arbitrary data
    data_end=datetime.datetime.today())

new_cases = rki.get_new(
    value = "confirmed",
    data_begin=datetime.datetime(2020,3,15),
    data_end=datetime.datetime.today())

# Create a multiplot
fig, axes = plt.subplots(2,1, figsize=(12,6))

# Plot the new cases onto axes[0]
cov19.plot._timeseries(
    x=new_cases.index,
    y=new_cases,
    ax=axes[0],
    what="model", #We define model here to get a line instead of data points
)

# Plot the new deaths onto axes[1]
cov19.plot._timeseries(
    x=new_deaths.index,
    y=new_deaths,
    ax=axes[1],
    what="model", #We define model here to get a line instead of data points
)

# Label the plots

axes[0].set_title("New cases")

axes[1].set_title("New deaths")

# Show the figure
fig.show()

Helper functions¶

covid19_inference.plot._get_array_from_trace_via_date(model, trace, var, start=None, end=None, dates=None)[source]¶

Parameters

model (model instance) –
trace (trace instance) –
var (str) – the variable name in the trace
start (datetime.datetime) – get all data for a range from start to end. (both boundary dates included)
end (datetime.datetime) –
dates (list of datetime.datetime objects, optional) – the dates for which to get the data. Default: None, will return all available data.

Returns

data (nd array, 3 dim) – the elements from the trace matching the dates. dimensions are as follows 0 samples, if no samples only one entry 1 data with time matching the returned dates (if compatible variable) 2 region, if no regions only one entry
dates (pandas DatetimeIndex) – the matching dates. this is essnetially an array of dates than can be passed to matplotlib

Example

import covid19_inference as cov
model, trace = cov.create_example_instance()
y, x = cov.plot._get_array_from_trace_via_date(
    model, trace, "lambda_t", model.data_begin, model.data_end
)
ax = cov.plot._timeseries(x, y[:,:,0], what="model")

covid19_inference.plot._new_cases_to_cum_cases(x, y, what, offset=0)[source]¶

so this conversion got ugly really quickly. need to check dimensionality of y

Parameters

x (pandas DatetimeIndex array) – will be padded accordingly
y (1d or 2d numpy array) – new cases matching dates in x. if 1d, we assume raw data (no samples) if 2d, we assume results from trace with 0th dim samples and 1st new cases matching x
what (str) – dirty workaround to differntiate between traces and raw data “data” or “trace”
offset (int or array like) – added to cum sum (should be the known cumulative case number at the first date of provided in x)

Returns

x_cum (pandas DatetimeIndex array) – dates of the cumulative cases
y_cum (nd array) – cumulative cases matching x_cum and the dimension of input y

Example

cum_dates, cum_cases = _new_cases_to_cum_cases(new_dates, new_cases)

covid19_inference.plot._label_for_varname(key)[source]¶

get the label for trace variable names (e.g. placed on top of distributions)

default for unknown keys is the key itself

Todo

add more parameters

covid19_inference.plot._math_for_varname(key)[source]¶

get the math string for trace variable name, e.g. used to print the median representation.

default for unknown keys is “$x$”

Todo

use regex

covid19_inference.plot._days_to_mpl_dates(days, origin)[source]¶

convert days as number to matplotlib compatible date numbers. this is not the same as pandas dateindices, but numpy operations work on them

Parameters

days (number, 1d array of numbers) – the day number to convert, e.g. integer values >= 0, one day per int
origin (datetime.datetime) – the date object corresponding to day 0

covid19_inference.plot._get_mpl_text_coordinates(text, ax)[source]¶

helper to get coordinates of a text object in the coordinates of the axes element [0,1]. used for the rectangle backdrop.

Returns: x_min, x_max, y_min, y_max

covid19_inference.plot._add_mpl_rect_around_text(text_list, ax, x_padding=0.05, y_padding=0.05, **kwargs)[source]¶

add a rectangle to the axes (behind the text)

provide a list of text elements and possible options passed to mpl.patches.Rectangle e.g. facecolor=”grey”, alpha=0.2, zorder=99,

covid19_inference.plot._rx_cp_id(key)[source]¶: get the change_point index from a compatible variable name

covid19_inference.plot._rx_hc_id(key)[source]¶: get the L1 / L2 value of hierarchical variable name

covid19_inference.plot._format_k(prec)[source]¶: format yaxis 10_000 as 10 k. _format_k(0)(1200, 1000.0) gives “1 k” _format_k(1)(1200, 1000.0) gives “1.2 k”

covid19_inference.plot._format_date_xticks(ax, minor=None)[source]¶

covid19_inference.plot._truncate_number(number, precision)[source]¶

covid19_inference.plot._string_median_CI(arr, prec=2)[source]¶

covid19_inference.plot._add_watermark(ax, mark='Dehning et al. 10.1126/science.abb9789')[source]¶: Add our arxive url to an axes as (upper right) title

class covid19_inference.plot.Param[source]¶

Paramters Base Class (a tweaked dict)

We inherit from dict and also provide keys as attributes, mapped to .get() of dict. This avoids the KeyError: if getting parameters via .the_parname, we return None when the param does not exist.

Avoid using keys that have the same name as class functions etc.

Example

foo = Param(lorem="ipsum")
print(foo.lorem)
>>> 'ipsum'
print(foo.does_not_exist is None)
>>> True

Variables saved in the trace¶

The trace by default contains the following parameters in the SIR/SEIR hierarchical model. XXX denotes a number.

Name in trace	Dimensions	Created by function
lambda_XXX_L1	samples	lambda_t_with_sigmoids/make_change_point_RVs
lambda_XXX_L2	samples x regions	lambda_t_with_sigmoids/make_change_point_RVs
sigma_lambda_XXX_L2	samples	lambda_t_with_sigmoids/make_change_point_RVs
transient_day_XXX_L1	samples	lambda_t_with_sigmoids/make_change_point_RVs
transient_day_XXX_L2	samples x regions	lambda_t_with_sigmoids/make_change_point_RVs
sigma_transient_day_XXX_L2	samples	lambda_t_with_sigmoids/make_change_point_RVs
transient_len_XXX_L1	samples	lambda_t_with_sigmoids/make_change_point_RVs
transient_len_XXX_L2	samples x regions	lambda_t_with_sigmoids/make_change_point_RVs
sigma_transient_len_XXX_L2	samples	lambda_t_with_sigmoids/make_change_point_RVs
delay_L1	samples	delay_cases
delay_L2	samples x regions	delay_cases
sigma_delay_L2	samples	delay_cases
weekend_factor_L1	samples	week_modulation
weekend_factor_L2	samples x regions	week_modulation
sigma_weekend_factor_L2	samples	week_modulation
offset_modulation	samples	week_modulation
new_cases_raw	samples x time x regions	week_modulation
mu	samples	SIR/SEIR
I_begin	samples x regions	SIR/SEIR
new_cases	samples x time x regions	SIR/SEIR
sigma_obs	samples x regions	SIR/SEIR
new_E_begin	samples x 11 x regions	SEIR
median_incubation_L1	samples	SEIR
median_incubation_L2	samples x regions	SEIR
sigma_median_incubation_L2	samples	SEIR

For the non-hierchical model, variables with _L2 suffixes are missing, and _L1 suffixes are removed from the name.

Contributing¶

We always welcome contributions. Here we gather some guidelines to make the process as smooth as possible.

Beginning¶

To see where help is needed, go to the issues page on Github. If you want to begin on an issue, make a comment below and begin a draft pull request: https://github.blog/2019-02-14-introducing-draft-pull-requests/ You can link the pull request on the right side of the commit to it.

When you have finished working on the issue, change it to a regular pull request. Check that there are no conflicts to the current master (https://www.digitalocean.com/community/tutorials/how-to-rebase-and-update-a-pull-request)

Code formatting¶

We use black https://github.com/psf/black as automatic code formatter. Please run your code through it before you open a pull request.

We do not check for formatting in the testing (travis) but have a config in the repository that uses black as a pre-commit hook.

This snippet should get you up and running:

conda install -c conda-forge black
conda install -c conda-forge pre-commit
pre-commit install

Try to stick to PEP 8. You can use type annotations if you want, but it is not necessary or encouraged.

Testing¶

We use travis and pytest. To check your changes locally:

python -m pytest --log-level=INFO --log-cli-level=INFO

It would be great if anything that is added to the code-base has an according test in the tests folder. We are not there yet, but it is on the todo. Be encouraged to add tests :)

Documentation¶

The documentation is built using Sphinx from the docstrings. To test it before submitting, navigate with a terminal to the docs/ directory. Install if necessary the packages listed in piprequirements.txt run make html. The documentation can then be accessed in docs/_build/html/index.html. As an example you can look at the documentation of covid19_inference.model.SIR()

Debugging¶

This is some pointer to help debugging models and sampling issues

General approach for nans/infs during sampling¶

The idea of this approach is to sample from the prior and then run the model. If the log likelihood is then -inf, there is a problem, and the output of the theano functions is inspected.

Sample from prior:

from pymc3.util import (
    get_untransformed_name,
    is_transformed_name)

varnames = list(map(str, model.vars))

for name in varnames:
    if is_transformed_name(name):
        varnames.append(get_untransformed_name(name))

with model:
    points = pm.sample_prior_predictive(var_names = varnames)
    points_list = []
    for i in range(len(next(iter(points.values())))):
        point_dict = {}
        for name, val in points.items():
                point_dict[name] = val[i]
        points_list.append(point_dict)

points_list is a list of the starting points for the model, sampled from the prior. Then to run the model and print the log-likelihood:

fn = model.fn(model.logpt)

for point in points_list[:]:
    print(fn(point))

To monitor the output and save it in a file (for use in ipython). Learned from: http://deeplearning.net/software/theano/tutorial/debug_faq.html#how-do-i-step-through-a-compiled-function

%%capture cap --no-stderr
def inspect_inputs(i, node, fn):
    print(i, node, "input(s) value(s):", [input[0] for input in fn.inputs],
          end='')

def inspect_outputs(i, node, fn):
    print(" output(s) value(s):", [output[0] for output in fn.outputs])

fn_monitor = model.fn(model.logpt,
                      mode=theano.compile.MonitorMode(
                           pre_func=inspect_inputs,
                           post_func=inspect_outputs).excluding(
                                 'local_elemwise_fusion', 'inplace'))

fn = model.fn(model.logpt)

for point in points_list[:]:
    if fn(point) < -1e10:
        print(fn_monitor(point))
        break

In a new cell:

with open('output.txt', 'w') as f:
    f.write(cap.stdout)

Then one can open output.txt in a text editor, and follow from where infs or nans come from by following the inputs and outputs up through the graph

Sampler: MCMC (Nuts)¶

Divergences¶

During sampling, a significant fraction of divergences are a sign that the sampler doesn’t sample the whole posterior. In this case the model should be reparametrized. See this tutorial for a typical example: https://docs.pymc.io/notebooks/Diagnosing_biased_Inference_with_Divergences.html

And these papers include some more details: https://pdfs.semanticscholar.org/7b85/fb48a077c679c325433fbe13b87560e12886.pdf https://arxiv.org/pdf/1312.0906.pdf

Bad initial energy¶

This typically occurs when some distribution in the model can’t be evaluated at the starting point of chain. Run this to see which distribution throws nans or infs:

for RV in model.basic_RVs:
    print(RV.name, RV.logp(model.test_point))

However, this is only evaluates the test_point. When PyMC3 starts sampling, it adds some jitter around this test_point, which then could lead to nans. Run this to add jitter and then evaluate the logp:

chains=4
for RV in model.basic_RVs:
    print(RV.name)

    for _ in range(chains):
        mean = {var: val.copy() for var, val in model.test_point.items()}
        for val in mean.values():
            val[...] += 2 * np.random.rand(*val.shape) - 1
        print(RV.logp(mean))

This code could potentially change in newer versions of PyMC3 (this is tested in 3.8). Read the source code, to know which random jitter PyMC3 currently adds at beginning.

Nans occur during sampling¶

Run the sampler with the debug mode of Theano.

from theano.compile.nanguardmode import NanGuardMode
mode = NanGuardMode(nan_is_error=True, inf_is_error=False, big_is_error=False,
                    optimizer='o1')
trace = pm.sample(mode=mode)

However this doesn’t lead to helpful messages if nans occur during gradient evaluations.

Sampler: Variational Inference¶

There exist some ways to track parameters during sampling. An example:

with model:
    advi = pm.ADVI()
    print(advi.approx.group)

    print(advi.approx.mean.eval())
    print(advi.approx.std.eval())

    tracker = pm.callbacks.Tracker(
        mean=advi.approx.mean.eval,  # callable that returns mean
        std=advi.approx.std.eval  # callable that returns std
    )

    approx = advi.fit(100000, callbacks=[tracker],
                      obj_optimizer=pm.adagrad_window(learning_rate=1e-3),)
                     #total_grad_norm_constraint=10) #constrains maximal gradient, could help


    print(approx.groups[0].bij.rmap(approx.params[0].eval()))

    plt.plot(tracker['mean'])
    plt.plot(tracker['std'])

For the tracker, the order of the parameters is saved in:

approx.ordering.by_name

and the indices encoded there in the slc field. To plot the mean value of a given parameter name, run:

plt.plot(np.array(tracker['mean'])[:, approx.ordering.by_name['name'].slc]

The debug mode is set with the following parameter:

from theano.compile.nanguardmode import NanGuardMode
mode = NanGuardMode(nan_is_error=True, inf_is_error=False, big_is_error=False,
                    optimizer='o1')
approx = advi.fit(100000, callbacks=[tracker],
             fn_kwargs={'mode':mode})

Getting Started¶

Installation¶

First Steps¶

Examples¶

Disclaimer¶

SIR — susceptible-infected-recovered¶

More Details¶

SEIR-like — susceptible-exposed-infected-recovered¶

More Details¶

References¶

More Details¶

Plotting¶

High level functions¶

Low level functions¶

Example¶

Helper functions¶

Variables saved in the trace¶

Contributing¶

Beginning¶

Code formatting¶

Testing¶

Documentation¶

Debugging¶

General approach for nans/infs during sampling¶

Sampler: MCMC (Nuts)¶

Divergences¶

Bad initial energy¶

Nans occur during sampling¶

Sampler: Variational Inference¶

Indices and tables¶