Package manual

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.

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 complete system state over time.
• 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.

SingleModelResults

The results of a single-model inference analysis.

Fields

• model – model names.
• ibis – primary analysis IBIS results.
• mcmc – validation analysis MCMC results.

ModelComparisonResults

The results of a model comparison, based on the Bayesian model evidence (BME.)

Fields

• names – model names.
• bme – matrix of BME estimates.
• mu – vector of BME estimate means (by model.)
• sigma – vector of BME estimate standard deviations.
• n_runs – number of independent analyses for each model.
• run_time – total application run time.

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.

Named parameters

• freq_dep – epidemiological models only, set to true for frequency-dependent contact rates.
• obs_error – average observation error (default = 2.)
• t0_index – index of the parameter that represents the initial time. 0 if fixed at 0.0.

model_name options

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

Examples

generate_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_smc2_analysis(model, obs_data; ... )

Run an SMC^2 (i.e. particle filter IBIS) analysis based on model and obs_data of type Observations.

Parameters

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

Optional

• np – number of [outer, i.e. theta] particles (default = 2000.)
• npf – number of [inner] particles (default = 200.)
• ess_rs_crit – resampling criteria (default = 0.5.)
• ind_prop – true for independent theta proposals (default = false.)
• alpha – user-defined, increase for lower acceptance rate targeted (default = 1.002.)

Example

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

run_mbp_ibis_analysis(model, obs_data; ... )

Run an MBP-IBIS analysis based on model, and obs_data of type Observations.

Parameters

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

Optional

• np – number of particles (default = 4000.)
• ess_rs_crit – resampling criteria (default = 0.5.)
• n_props – MBP mutations per step (default = 3.)
• ind_prop – true for independent theta proposals (default = false.)
• alpha – user-defined, increase for lower acceptance rate targeted (default = 1.002.)

Example

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

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_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.

resample(sample; n = 10000 [, discard = 0])

Resample a results type of some kind, e.g. MCMCSample, ImportanceSample or SingleModelResults.

Parameters

• sample – parameter inference results of some kind.
• n – the index of the model parameter to be plotted on the x axis.
• discard – number of initial samples to be discarded - applicable to sample of type RejectionSamples only.

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(results; [null_index = 1], [display=true])

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. The display flag can be set =false to return a single DataFrame for valid types.

plot_trajectory(x; plot_index=[:])

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

The only input parameter required is x of type SimResults, i.e. from a call to gillespie_sim. All system states are plotted by default, but a subset can be specified by passing an integer array to the plot_index option, which contains the indices of the desired subset. E.g. [1,2] for the first two 'compartments' only.

plot_observations(x; plot_index=1)

Plot the trajectory of observation values for one or more [simulated or real] data sets using UnicodePlots.jl.

The only input parameter required is x of type SimResults, i.e. from a call to gillespie_sim. The first observation value is plotted by default, but another can be specified by passing an integer to the plot_index option.

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)

Plot the Bayesian model evidence (BME) from a [multi-] model inference workflow, using UnicodePlots.jl.

Parameters

• results – ModelComparisonResults, i.e. from a call to run_model_comparison_analysis.