The iterated filtering method,
called IF2
in mcstate
, is an alternative
method to infer the most likely parameters. Instead of conducting an
MCMC over particle filter likelihoods, as in the pmcmc
methods, the IF2
algorithm proceeds roughly as follows:
Over the course of the iterations, the parameters closer to the data tend to be selected in the resampling step. The noise allows them to vary and approach an optimum. The cooling schedule is similar to simulated annealing, allowing larger steps at the beginning, and finer steps as the optimum is approached.
This algorithm is described in the following paper:
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.
Based on the example in vignette("sir_models")
, we will
start by loading an SIR model and some simulated data. See that vignette
for more details on what is being done here
A data set to fit to
incidence <- read.csv(system.file("sir_incidence.csv", package = "mcstate"))
plot(cases ~ day, incidence,
type = "o", xlab = "Day", ylab = "New cases", pch = 19)
data <- mcstate::particle_filter_data(incidence, time = "day", rate = 4)
#> Warning in mcstate::particle_filter_data(incidence, time = "day", rate = 4):
#> 'initial_time' should be provided. I'm assuming '0' which is one time unit
#> before the first time in your data (1), but this might not be appropriate. This
#> will become an error in a future version of mcstate
A dust model generator object
A comparison function
compare <- function(state, observed, pars = NULL) {
if (is.na(observed$cases)) {
return(NULL)
}
exp_noise <- 1e6
incidence_modelled <- state[1, , drop = TRUE]
incidence_observed <- observed$cases
lambda <- incidence_modelled +
rexp(n = length(incidence_modelled), rate = exp_noise)
dpois(x = incidence_observed, lambda = lambda, log = TRUE)
}
An index function for filtering the run state
The differences for the IF2 method start here! We now need to specify the initial values and the search range for the two parameters β and γ:
pars <- mcstate::if2_parameters$new(
list(mcstate::if2_parameter("beta", 0.5, min = 0, max = 1),
mcstate::if2_parameter("gamma", 0.01, min = 0, max = 1)))
For this example, we start far away from the true values of β = 0.2 and γ = 0.1.
The IF2 algorithm needs a number of parameter sets to perturb, the
magnitude of the noise with which to perturb them, the amount to reduce
the perturbations by, and the number of iterations of which to perform
this reduction in magnitude. As with our mcstate::pmcmc
method, we organise these parameters into a single control object:
control <- mcstate::if2_control(
pars_sd = list("beta" = 0.02, "gamma" = 0.02),
iterations = n_par_sets,
n_par_sets = n_iterations,
cooling_target = 0.5,
progress = interactive())
Here we will iterate 10 parameter sets with
X ~ N(0, 0.02)
, with geometrically decreasing noise, until
after 30 iterations the magnitude of the noise has halved.
To start the IF2, we use a particle filter object in order to do some validation against the data. This may change in future and the particle filter is not actually run at this point in the process
Then, call mcstate::if2
, which will take a while to
run
We can plot the likelihood over the course of the iterations
plot(res$result$log_likelihood, main = "LL profile",
xlab = "IF iteration", ylab = "Log likelihood", type = "l")
and the values of β and γ (using a small helper function)
plot_if2_parameter <- function(name, obj) {
pars <- obj$result$pars
quantiles <- apply(pars[name, , ], 2, quantile, c(0.025, 0.975))
plot(apply(pars[name, , ], 2, mean), type = "l", lwd = 1, col = "#000000",
xlab = "IF iteration", ylab = name, ylim = range(quantiles))
matlines(t(quantiles), type = "l", lty = 2, lwd = 1, col = "#999999")
legend("topright", lwd = 1, legend = c("Mean", "95% quantile"), bty = "n")
}
plot_if2_parameter("beta", res)
Checking the mean values of each parameter at the final iteration, they are close to their true values:
pars_series <- res$result$pars
rowMeans(pars_series[, , dim(pars_series)[[3]]])
#> beta gamma
#> 0.3023458 0.1980262
We can also run particle filters at each final parameter estimate to get an estimate of the model likelihood, showing the mean value in red: