Run a SMC^2. This is experimental and subject to change. Use at your own risk.
smc2(pars, filter, control)A smc2_parameters object containing
information about parameters (ranges, priors, proposal kernel,
translation functions for use with the particle filter).
A particle_filter object
A smc2_control object to control the behaviour of the algorithm
A smc2_result object, with elements
pars: a matrix of sampled parameters (n_parameter_set long)
probabilities: a matrix of probabilities (log_prior, log_likelihood,
log_posterior and weight). The latter is the log posterior
normalised over all samples
statistics: interesting or useful statistics about your sample,
including the ess (effective sample size, over time),
acceptance_rate (where a regeneration step was done, the acceptance
rate), n_particles, n_parameter_sets and n_steps (inputs to
the simulation). The effort field is a rough calculation of the
number of particle-filter runs that this run was worth.
# We use an example from dust which implements an epidemiological SIR
# (Susceptible-Infected-Recovered) model; see vignette("sir_models")
# for more background, as this example follows from the pMCMC example
# there
# The key tuning here is the number of particles per filter and number
# of parameter sets to consider simultaneously. Ordinarily these would
# be set (much) higher with an increase in computing time
n_particles <- 42
n_parameter_sets <- 20
# Basic epidemiological (Susceptible-Infected-Recovered) model
sir <- dust::dust_example("sir")
# Pre-computed incidence data
incidence <- read.csv(system.file("sir_incidence.csv", package = "mcstate"))
# Annotate the data so that it is suitable for the particle filter to use
dat <- mcstate::particle_filter_data(incidence, "day", 4, 0)
# Subset the output during run
index <- function(info) {
list(run = 5L)
}
# The comparison function, used to compare simulated data with observe
# data, given the above subset
compare <- function(state, observed, pars) {
exp_noise <- 1e6
incidence_modelled <- state[1L, , 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)
}
# Finally, construct the particle filter:
filter <- mcstate::particle_filter$new(dat, sir, n_particles, compare,
index = index)
# To control the smc2 we need to specify the parameters to consider
pars <- mcstate::smc2_parameters$new(
list(
mcstate::smc2_parameter("beta",
function(n) runif(n, 0, 1),
function(x) dunif(x, 0, 1, log = TRUE),
min = 0, max = 1),
mcstate::smc2_parameter("gamma",
function(n) runif(n, 0, 1),
function(x) dunif(x, 0, 1, log = TRUE),
min = 0, max = 1)))
control <- mcstate::smc2_control(n_parameter_sets, progress = TRUE)
# Then we run the particle filter
res <- mcstate::smc2(pars, filter, control)
#> Finished 100 steps in 8 secs
# This returns quite a lot of information about the fit, and this will
# change in future versions
res
#> $pars
#> beta gamma
#> [1,] 0.3093155 0.18301115
#> [2,] 0.1911669 0.09043611
#> [3,] 0.3907869 0.22471620
#> [4,] 0.2153920 0.11965343
#> [5,] 0.2029480 0.08585487
#> [6,] 0.1386245 0.03068581
#> [7,] 0.4913324 0.30836149
#> [8,] 0.6128431 0.27551790
#> [9,] 0.1574880 0.02922703
#> [10,] 0.2165904 0.09988620
#> [11,] 0.2504913 0.08748742
#> [12,] 0.2349047 0.12025240
#> [13,] 0.4534552 0.14185847
#> [14,] 0.6092114 0.23627732
#> [15,] 0.5464832 0.21336094
#> [16,] 0.5354779 0.21819638
#> [17,] 0.3819656 0.11899958
#> [18,] 0.4627995 0.18108103
#> [19,] 0.4365939 0.18008579
#> [20,] 0.2299612 0.08552151
#>
#> $probabilities
#> log_prior log_likelihood log_posterior weight
#> [1,] 0 -252.6342 -252.6342 9.372691e-05
#> [2,] 0 -249.0119 -249.0119 3.507459e-03
#> [3,] 0 -255.6199 -255.6199 4.733299e-06
#> [4,] 0 -246.4818 -246.4818 4.403647e-02
#> [5,] 0 -247.3477 -247.3477 1.852518e-02
#> [6,] 0 -261.2342 -261.2342 1.725602e-08
#> [7,] 0 -257.6297 -257.6297 6.343559e-07
#> [8,] 0 -286.1875 -286.1875 2.510975e-19
#> [9,] 0 -282.5354 -282.5354 9.680940e-18
#> [10,] 0 -246.1817 -246.1817 5.944924e-02
#> [11,] 0 -266.1283 -266.1283 1.292600e-10
#> [12,] 0 -243.4933 -243.4933 8.743790e-01
#> [13,] 0 -358.7733 -358.7733 7.519999e-51
#> [14,] 0 -358.2318 -358.2318 1.292459e-50
#> [15,] 0 -294.4855 -294.4855 6.252750e-23
#> [16,] 0 -288.1854 -288.1854 3.405467e-20
#> [17,] 0 -320.4043 -320.4043 3.464707e-34
#> [18,] 0 -296.8766 -296.8766 5.722980e-24
#> [19,] 0 -277.0725 -277.0725 2.282691e-15
#> [20,] 0 -255.9242 -255.9242 3.491769e-06
#>
#> $statistics
#> $statistics$ess
#> [1] 15.816943 11.301591 9.580949 12.655357 10.057373 9.477349 12.149554
#> [8] 10.798700 8.500094 13.220614 11.093068 6.046357 10.082006 8.858798
#> [15] 7.945963 5.354276 3.833227 2.721528 3.411812 4.620651 2.272063
#> [22] 3.622429 1.321429 1.971645 3.284305 4.251124 3.151735 2.084010
#> [29] 1.816593 2.399074 3.860173 3.669836 8.469251 6.198018 15.683858
#> [36] 12.611303 11.676231 11.033467 10.192990 9.267872 18.603920 15.985515
#> [43] 17.221367 12.024333 11.282999 8.538440 16.909230 17.182555 16.076947
#> [50] 14.234458 13.000213 12.333011 11.462462 10.655559 9.411526 16.910551
#> [57] 17.166109 15.822495 15.190905 12.881715 12.671453 12.393773 11.654590
#> [64] 11.407807 10.708538 9.881198 15.493470 17.523067 16.390716 15.330937
#> [71] 13.233933 13.588038 12.926980 13.933442 13.238797 11.948426 10.325479
#> [78] 10.164376 8.878644 19.193786 16.739856 14.211743 16.487422 15.942664
#> [85] 14.696732 14.612732 13.887976 12.566439 11.443178 9.613109 5.505015
#> [92] 13.551766 14.985555 14.855052 15.349661 12.438938 11.623082 12.637786
#> [99] 9.934752 18.390931
#>
#> $statistics$acceptance_rate
#> [1] NA NA 0.75 NA NA 0.60 NA NA 0.50 NA NA 0.45 NA 0.55 0.45
#> [16] 0.50 0.70 0.30 0.50 0.30 0.35 0.30 0.45 0.40 0.35 0.45 0.45 0.50 0.55 0.50
#> [31] 0.60 0.55 0.50 0.60 NA NA NA NA NA 0.50 NA NA NA NA NA
#> [46] 0.30 NA NA NA NA NA NA NA NA 0.35 NA NA NA NA NA
#> [61] NA NA NA NA NA 0.05 NA NA NA NA NA NA NA NA NA
#> [76] NA NA NA 0.40 NA NA NA NA NA NA NA NA NA NA 0.35
#> [91] 0.25 NA NA NA NA NA NA NA 0.50 NA
#>
#> $statistics$n_particles
#> [1] 42
#>
#> $statistics$n_parameter_sets
#> [1] 20
#>
#> $statistics$n_steps
#> [1] 100
#>
#> $statistics$effort
#> [1] 240
#>
#>
#> attr(,"class")
#> [1] "smc2_result"
# Most useful is likely the predict method:
predict(res)
#> beta gamma
#> [1,] 0.2349047 0.1202524
#> [2,] 0.2349047 0.1202524
#> [3,] 0.2153920 0.1196534
#> [4,] 0.2165904 0.0998862
#> [5,] 0.2349047 0.1202524
#> [6,] 0.2349047 0.1202524
#> [7,] 0.2153920 0.1196534
#> [8,] 0.2349047 0.1202524
#> [9,] 0.2349047 0.1202524
#> [10,] 0.2349047 0.1202524
#> [11,] 0.2349047 0.1202524
#> [12,] 0.2349047 0.1202524
#> [13,] 0.2349047 0.1202524
#> [14,] 0.2349047 0.1202524
#> [15,] 0.2349047 0.1202524
#> [16,] 0.2349047 0.1202524
#> [17,] 0.2349047 0.1202524
#> [18,] 0.2349047 0.1202524
#> [19,] 0.2349047 0.1202524
#> [20,] 0.2349047 0.1202524