Package manual

DPOMPModel

A mutable struct which represents a DSSCT model (see Models for further details).

Fields

• model_name – string, e,g, "SIR".
• rate_function – event rate function.
• initial_condition – initial condition.
• m_transition – transition matrix.
• `obs_function – observation function, use this to add 'noise' to simulated observations.
• obs_model – observation model likelihood function.
• prior – prior [multivariate] Distributions.Distribution.
• t0_index – index of the parameter that represents the initial time. 0 if fixed at 0.0.

Particle

E.g. the main results of a simulation including the initial and final conditions, but not the full state trajectory.

Fields

• theta – e.g. simulation parameters.
• initial_condition – initial system state.
• final_condition – final system state.
• trajectory – the event history.
• log_like – trajectory log likelihood, mainly for internal use.

Event

Requires no explanation.

Fields

• time – the time of the event.
• event_type – indexes the rate function and transition matrix.

Observation

A single observation. Note that by default val has the same size as the model state space. However that is not necessary - it need only be compatible with the observation model.

Fields

• time – similar to Event.time, the time of the observation.
• obs_id – <1 if not a resampling step.
• prop – optional information for the observation model.
• val – the observation value.

ARQModel

Contains the PDF (or estimate, or approximation of the target density - a function) and parameter density, which together specify the model.

Fields

• pdf – prior density function.
• sample_interval – An array specifying the (fixed or fuzzy) interval between samples.

SimResults

The results of a simulation, including the full state trajectory.

Fields

• model_name – string, e,g, "SIR".
• particle – the 'trajectory' variable, of type Particle.
• population – records the final system state.
• observations – simulated observations data (an Array of Observation types.)

ImportanceSample

The results of an importance sampling analysis, such as iterative batch importance sampling algorithms.

Fields

• mu – weighted sample mean.
• cv – weighted covariance.
• theta – two dimensional array of samples, e.g. parameter; iteration.
• weight – sample weights.
• run_time – application run time.
• bme – Estimate (or approximation) of the Bayesian model evidence.

RejectionSample

Essentially, the main results of an MCMC analysis, consisting of samples, mean, and covariance matrix.

Fields

• samples – three dimensional array of samples, e.g. parameter; iteration; Markov chain.
• mu – sample mean.
• cv – sample covariance matrix.

MCMCSample

The results of an MCMC analysis, mainly consisting of a RejectionSample.

Fields

• samples – samples of type RejectionSample.
• adapt_period – adaptation (i.e. 'burn in') period.
• sre – scale reduction factor estimate, i.e. Gelman diagnostic.
• run_time – application run time.

ARQMCMCSample

The results of an ARQ MCMC analysis including the ImportanceSample and resampled RejectionSample.

The sre scale factor reduction estimates relate the rejection (re)samples to the underlying importance sample.

Fields

• imp_sample – main results, i.e. ImportanceSample.
• samples – resamples, of type RejectionSample.
• adapt_period – adaptation (i.e. 'burn in') period.
• sample_dispersal – number of distinct [possible] sample values along each dimension in the unit cube.
• sample_limit – maximum number of samples per theta tupple.
• grid_range – bounds of the parameter space.
• sre – scale reduction factor estimate, i.e. Gelman diagnostic. NB. only valid for resamples.
• run_time – application run time.
• sample_cache – a link to the underlying likelihood cache - can be reused.

generate_model(model_name, initial_condition; freq_dep = false, obs_error = 2.0)

Generates an DPOMPModel instance. Observation models are generated using the partial_gaussian_obs_model function, with `σ = obs_error (see that functions entry for further details.)

Parameters

• model_name – the model, e.g. "SI"; "SIR"; "SEIR"; etc
• initial_condition – initial condition.

Optional parameters

• freq_dep – epidemiological models only, set to true for frequency-dependent contact rates.
• obs_error – average observation error (default = 2.)

model_name options

• "SI"
• "SIR"
• "SIS"
• "SEI"
• "SEIR"
• "SEIS"
• "SEIRS"
• "PREDPREY"
• "ROSSMAC"

Examples

generate_model("SIS", [100,1])

generate_custom_model(model_name, rate_function, initial_condition, m_transition; ... )

Generates an DPOMPModel instance. Observation models are generated using the partial_gaussian_obs_model function, with `σ = obs_error (see that functions entry for further details.)

Parameters

• model_name – the model, e.g. "SIR", "SEIR-custom", etc.
• rate_function – event rate function.
• initial_condition – initial condition
• m_transition – transition matrix.

Optional parameters

• `observation_function – observation function, use this to add 'noise' to simulated observations.
• obs_error – average observation error (default = 2.)
• obs_model – use this option to manuallu specify the observation model likelihood function.
• prior_density – prior density function.
• t0_index – index of the parameter that represents the initial time. 0 if fixed at 0.0.

Examples

generate_custom_model("SIS", [100,1])

partial_gaussian_obs_model(σ = 2.0; seq = 2, y_seq = seq)

Generate a simple Gaussian observation model. So-called because the accuracy of observations is 'known' and [assumed to be] normally distributed according to~N(0, σ), where observation error σ can be specified by the user.

Parameters

• σ – observation error.
• seq – the indexing sequence of the observed state, e.g. 2 for that state only, 3:4 for the third and fourth, etc.
• y_seq – as above, the corresponding [indexing] values for the observations data, seq unless otherwise specified.

test latex eqn:

\[rac{n!}{k!(n - k)!} = inom{n}{k}\]

Examples

p = partial_gaussian_obs_model(1.0, seq = 2)

gillespie_sim(model, parameters; tmax = 100.0, num_obs = 5)

Run a Doob-Gillespie (DGA) simulation based on model.

Returns a SimResults type containing the trajectory and observations data, or an array of the same if n_sims > 1.

Parameters

• model – DPOMPModel (see [DCTMPs.jl models]@ref).
• parameters – model parameters.

Optional

• tmax – maximum time (default: 100.)
• n_obs – the number of observations to draw (default: 5.)
• n_sims – number of simulations to draw (default: 1.)

Example

using DiscretePOMP
m = generate_model("SIR", [50, 1, 0])
x = DiscretePOMP.gillespie_sim(model, [0.005, 0.12])
println(DiscretePOMP.plot_trajectory(x))

run_mcmc_analysis(model, obs_data; ... )

Run an n_chains-MCMC analysis using the designated algorithm (MBP-MCMC by default.)

The initial_parameters are sampled from the prior distribution unless otherwise specified by the user. A Gelman-Rubin convergence diagnostic is automatically carried out (for n_chains > 1) and included in the [multi-chain] analysis results.

Parameters

• model – DPOMPModel (see [DCTMPs.jl models]@ref).
• obs_data – Observations data.

Optional

• n_chains – number of Markov chains (default: 3.)
• initial_parameters – 2d array of initial model parameters. Each column vector correspondes to a single model parameter.
• steps – number of iterations.
• adapt_period – number of discarded samples.
• mbp – model based proposals (MBP). Set mbp = false for standard proposals.
• ppp – the proportion of parameter (vs. trajectory) proposals in Gibbs sampler. Default: 30%. NB. not required for MBP.
• fin_adapt – finite adaptive algorithm. The default is false, i.e. [fully] adaptive.
• mvp – increase for a higher proportion of 'move' proposals. NB. not applicable if MBP = true (default: 2.)

Example

y = x.observations                          # some simulated data
model = generate_model("SIR", [50, 1, 0])   # a model
results = run_mcmc_analysis(model, y; fin_adapt = true) # finite-adaptive MCMC
tabulate_results(results)                   # optionally, show the results

run_ibis_analysis(model, obs_data; ... )

Run an iterative-batch-importance-sampling (IBIS) analysis based on model and obs_data of type Observations.

The default algorithm is SMC^2 (i.e. particle filter IBIS), the other option is model-based-proposal IBIS (use algorithm = "MBPI".) Note that the default value of the optional parameters below is sometimes affected by choice of algorithm. However these are overridden when specified by the user.

Parameters

• model – DPOMPModel (see [DCTMPs.jl models]@ref).
• obs_data – Observations data.

Optional parameters

• algorithm – String, per above.
• np – number of [outer, i.e. theta] particles.

Algorithm options

• ess_rs_crit – resampling criteria.
• ind_prop – true for independent theta proposals.
• alpha – increase to lower the targeted acceptance rate (default = 1.002.)
• npf – number of [inner] particles for SMC^2 (default = 200.)
• n_props – MBP mutations per step (default: 3.)

Example

# NB. using 'y' and 'model' as above
results = run_ibis_analysis(model, y)   # importance sample
tabulate_results(results)               # show the results

run_arq_mcmc_analysis(model, obs_data, theta_range; ... )

Run ARQ-MCMC analysis with n_chains Markov chains.

The Gelman-Rubin convergence diagnostic is computed automatically.

Parameters

• model – DPOMPModel (see docs.)
• obs_data – Observations data.
• sample_interval – An array specifying the (fixed or fuzzy) interval between samples.

Optional

• sample_dispersal – i.e. the length of each dimension in the importance sample.
• sample_limit – sample limit, should be increased when the variance of model.pdf is high (default: 1.)
• n_chains – number of Markov chains (default: 3.)
• steps – number of iterations.
• burnin – number of discarded samples.
• tgt_ar – acceptance rate (default: 0.33.)
• np – number of SMC particles in PF (default: 200.)
• ess_crit – acceptance rate (default: 0.33.)
• sample_cache – the underlying model likelihood cache - can be retained and reused for future analyses.

run_model_comparison_analysis(model, obs_data; ... )

Run n_runs independent analyses for each DPOMPModel element in models, and compute [estimate] the Bayesian model evidence (BME.)

Returns an object of type ModelComparisonResults, which includes the mean and standard deviation of the estimates obtained.

Parameters

• models – An Array of DPOMPModel.
• obs_data – An array of type Observation.

Optional parameters

• n_runs – number of independent analyses used to estimate the BME for each model.
• algorithm – String representing the inference method used for the analysis, "SMC2" for SMC^2 (default); or "MBPI" for MBP-IBIS.

Optional inference algorithm parameters

• np – number of [outer, i.e. theta] particles used in IBIS procedures (doesn't apply to ARQ-MCMC.)
• ess_rs_crit – Effective sample size (ESS) resampling criteria.
• npf – number of particles used in particle filter (doesn't apply to MBP-IBIS.)
• n_props – see the docs for run_mbp_ibis_analysis.

Example

# NB. first define some models to compare, e.g. as m1, m2, etc
models = [m1, m2, m3]
results = run_model_comparison_analysis(models, y; n_runs = 10)
tabulate_results(results)   # show the results (optional)

get_observations(source)

Return an array of type Observation, based on a two-dimensional array, DataFrame or file location (i.e. String.)

Note that a observation times must be in the first column of the input variable.

tabulate_results

Display the results of an inference analysis.

The main parameter is results – a data structure of type MCMCSample, ImportanceSample, ARQMCMCSample or ModelComparisonResults. When invoking the latter, the named parameter null_index = 1 by default but can be overridden. This determines the 'null' model, used to compute the Bayes factor.

print_results

Print the results of a parameter inference analysis to file.

Parameters

• samples – a data structure of type SimResults, ImportanceSample, MCMCSample or ARQMCMCSample.
• dpath – the directory where the results will be saved.

get_particle_filter_lpdf(model, obs_data; ... )

Generate a SMC [log] likelihood estimating function for a DPOMPModel – a.k.a. a particle filter.

Parameters

• model – DPOMPModel (see [DCTMPs.jl models]@ref).
• obs_data – Observations data.

Named parameters

• np – number of particles.
• rs_type – resampling method: 2 for stratified, 3 for multinomial, or 1 for systematic (the default.)
• ess_rs_crit – effective-sample-size resampling criteria.

```

plot_trajectory(x)

Plot the trajectory of a a DGA simulation using UnicodePlots.jl.

The only input parameter required is x of type SimResults, i.e. from a call to gillespie_sim.

plot_parameter_trace(mcmc, [parameter::Int64])

Produce a trace plot of samples using UnicodePlots.jl.

The mcmc input is of type MCMCSample, ARQMCMCSample or RejectionSample. The parameter index can be optionally specified, else all parameters are plotted and returned as an Array of unicode plots.

plot_parameter_marginal(sample, parameter)

Plot the marginal distribution of samples from an MCMC analysis for a given model parameter using UnicodePlots.jl.

Parameters

• results – Results object, e.g. of type MCMCSample.
• parameter – the index of the model parameter to be plotted.
• adapt_period– Adaptation period to be discarded, only required for RejectionSample.

Optional

• use_is – Resample IS rather than using MCMC [re]samples (ARQMCMCSample results only.)

plot_parameter_heatmap(mcmc, x_parameter, y_parameter)

Plot the marginal distribution of samples from an MCMC analysis for two model parameters using UnicodePlots.jl.

Parameters

• mcmc – MCMCResults, e.g. from a call to run_met_hastings_mcmc.
• x_parameter – the index of the model parameter to be plotted on the x axis.
• y_parameter – the index of the model parameter to be plotted on the y axis.

plot_model_comparison(results; boxplot = true)

Plot the Bayesian model evidence (BME) from a model comparison analysis, using UnicodePlots.jl.

Parameters

• results – ModelComparisonResults, i.e. from a call to run_model_comparison_analysis.
• boxplot – true for a series of boxplots, else a simple UnicodePlots.barplot showing only the average BME for each model variant (default.)

run_custom_mcmc_analysis(model, obs_data, trajectory_prop, [x0_prop], ... )

Run an n_chains data-augmented MCMC analysis, based on the Gibbs sampler with a user defined proposal function.

A function for conditionally sampling the initial trajectory variable can optionally be specified, use the Doob-Gillespie algorithm is used by default.

Elsewise, this function is equivalent to calling runmcmcanalysis with mbp = false, which invokes the standard Gibbs sampler.

Parameters

• model – DPOMPModel (see [DCTMPs.jl models]@ref).
• obs_data – Observations data.
• trajectory_prop – .
• x0_prop – Initial state variable sampler, DGA by default.

Optional

• n_chains – number of Markov chains (optional, default: 3.)
• initial_parameters – 2d array of initial model parameters. Each column vector correspondes to a single model parameter.
• steps – number of iterations.
• adapt_period – number of discarded samples.
• ppp – the proportion of parameter (vs. trajectory) proposals in Gibbs sampler. Default: 0.3, or 30%.
• fin_adapt – finite adaptive algorithm. The default is false, i.e. [fully] adaptive.

generate_custom_particle(model, obs_data, trajectory_prop, [x0_prop], ... )

For use with run_custom_mcmc_analysis(). Initialises a Particle based on an array of type Event. Also evaluates the likelihood function.

Parameters

• model – DPOMPModel (see [DCTMPs.jl models]@ref).
• obs_data – Observations data.
• trajectory – Array of Event types.

Optional

• theta – model parameters, sampled from prior unless otherwise specified.