Create an IF2 object for running and interacting with an IF2 inference.
if2(pars, filter, control)
if2_sample(obj, n_particles)An if2_parameters object, describing the parameters that will be varied in the simulation, and the method of transformation into model parameters.
A particle_filter object. We don't use
the particle filter directly (except for sampling with
mcstate::if2_sample) but this shares so much validation that
it's convenient. Be sure to set things like the seed and number
of threads here if you want to use anything other than the
default.
An if2_control() object
An object of class if2_fit, returned by mcstate::if2()
The number of particles to simulate, for each IF2 parameter set
An object of class if2_fit, which contains the sampled
parameters (over time) and their log-likelihoods
See: Ionides EL, Nguyen D, Atchadé Y, Stoev S, King AA (2015). "Inference for Dynamic and Latent Variable Models via Iterated, Perturbed Bayes Maps." PNAS, 112(3), 719–724. https://doi.org/10.1073/pnas.1410597112.
# A basic SIR model used in the particle filter example
gen <- dust::dust_example("sir")
# Some data that we will fit to, using 1 particle:
sir <- gen$new(pars = list(), time = 0, n_particles = 1)
dt <- 1 / 4
day <- seq(1, 100)
incidence <- rep(NA, length(day))
true_history <- array(NA_real_, c(5, 1, 101))
true_history[, 1, 1] <- sir$state()
for (i in day) {
state_start <- sir$state()
sir$run(i / dt)
state_end <- sir$state()
true_history[, 1, i + 1] <- state_end
# Reduction in S
incidence[i] <- state_start[1, 1] - state_end[1, 1]
}
# Convert this into our required format:
data_raw <- data.frame(day = day, incidence = incidence)
data <- particle_filter_data(data_raw, "day", 4, 0)
# A comparison function
compare <- function(state, observed, pars = NULL) {
if (is.null(pars$exp_noise)) {
exp_noise <- 1e6
} else {
exp_noise <- pars$exp_noise
}
incidence_modelled <- state[1,]
incidence_observed <- observed$incidence
lambda <- incidence_modelled +
rexp(length(incidence_modelled), exp_noise)
dpois(incidence_observed, lambda, log = TRUE)
}
# Range and initial values for model parameters
pars <- mcstate::if2_parameters$new(
list(mcstate::if2_parameter("beta", 0.15, min = 0, max = 1),
mcstate::if2_parameter("gamma", 0.05, min = 0, max = 1)))
# Set up of IF2 algorithm (the iterations and n_par_sets should be
# increased here for any real use)
control <- mcstate::if2_control(
pars_sd = list("beta" = 0.02, "gamma" = 0.02),
iterations = 10,
n_par_sets = 40,
cooling_target = 0.5,
progress = interactive())
# Create a particle filter object
filter <- mcstate::particle_filter$new(data, gen, 1L, compare)
# Then run the IF2
res <- mcstate::if2(pars, filter, control)
# Get log-likelihood estimates from running a particle filter at
# each final parameter estimate
ll_samples <- mcstate::if2_sample(res, 20)