PMCMC | R Documentation |
Runs particle Markov chain Monte Carlo (PMCMC) algorithm using a bootstrap particle filter for fitting infectious disease models to time series count data.
PMCMC(x, ...) ## S3 method for class 'PMCMC' PMCMC( x, niter = 1000, nprintsum = 100, adapt = TRUE, adaptmixprop = 0.05, nupdate = 100, ... ) ## Default S3 method: PMCMC( x, priors, func, u, npart = 100, iniPars = NA, fixpars = FALSE, niter = 1000, nprintsum = 100, adapt = TRUE, propVar = NA, adaptmixprop = 0.05, nupdate = 100, ... )
x |
A |
... |
Not used here. |
niter |
An integer specifying the number of iterations to run the MCMC. |
nprintsum |
Prints summary of MCMC to screen every |
adapt |
Logical determining whether to use adaptive proposal or not. |
adaptmixprop |
Mixing proportion for adaptive proposal. |
nupdate |
Controls when to start adaptive update. |
priors |
A |
func |
A
|
u |
A named vector of initial states. |
npart |
An integer specifying the number of particles for the bootstrap particle filter. |
iniPars |
A named vector of initial values for the parameters of the model. If left unspecified, then these are sampled from the prior distribution(s). |
fixpars |
A logical determining whether to fix the input parameters (useful for determining the variance of the marginal likelihood estimates). |
propVar |
A numeric (npars x npars) matrix with log (or logistic) covariances to use
as (initial) proposal matrix. If left unspecified then defaults to
|
Function runs a particle MCMC algorithm using a bootstrap particle filter for a given model.
If running with fixpars = TRUE
then this runs niter
simulations
using fixed parameter values. This can be used to optimise the number of
particles after a training run. Also has print()
, summary()
,
plot()
, predict()
and window()
methods.
If the code throws an error, then it returns a missing value (NA
). If
fixpars = TRUE
it returns a list of length 2 containing:
output
: a matrix with two columns. The first contains the simulated
log-likelihood, and the second is a binary indicator relating to whether the
simulation was 'skipped' or not (1 = skipped, 0 = not skipped);
pars
: a vector of parameters used for the simulations.
If fixpars = FALSE
, the routine returns a PMCMC
object, essentially a
list
containing:
pars
: an mcmc
object containing posterior samples for the parameters;
u
: a copy of the u
input;
accrate
: the cumulative acceptance rate;
npart
: the chosen number of particles;
time
: the time taken to run the routine (in seconds);
propVar
: the proposal covariance for the parameter updates;
data
: a copy of the x
input;
priors
: a copy of the priors
input;
func
: a copy of the func
input.
Andrieu C, Doucet A and Holenstein R (2010) <doi:10.1111/j.1467-9868.2009.00736.x>
print.PMCMC
, plot.PMCMC
, predict.PMCMC
, summary.PMCMC
window.PMCMC
## set up data to pass to PMCMC flu_dat <- data.frame( t = 1:14, Robs = c(3, 8, 26, 76, 225, 298, 258, 233, 189, 128, 68, 29, 14, 4) ) ## set up observation process obs <- data.frame( dataNames = "Robs", dist = "pois", p1 = "R + 1e-5", p2 = NA, stringsAsFactors = FALSE ) ## set up model (no need to specify tspan ## argument as it is set in PMCMC()) transitions <- c( "S -> beta * S * I / (S + I + R + R1) -> I", "I -> gamma * I -> R", "R -> gamma1 * R -> R1" ) compartments <- c("S", "I", "R", "R1") pars <- c("beta", "gamma", "gamma1") model <- mparseRcpp( transitions = transitions, compartments = compartments, pars = pars, obsProcess = obs ) ## set priors priors <- data.frame( parnames = c("beta", "gamma", "gamma1"), dist = rep("unif", 3), stringsAsFactors = FALSE) priors$p1 <- c(0, 0, 0) priors$p2 <- c(5, 5, 5) ## define initial states iniStates <- c(S = 762, I = 1, R = 0, R1 = 0) set.seed(50) ## run PMCMC algorithm post <- PMCMC( x = flu_dat, priors = priors, func = model, u = iniStates, npart = 25, niter = 5000, nprintsum = 1000 ) ## plot MCMC traces plot(post, "trace") ## continue for some more iterations post <- PMCMC(post, niter = 5000, nprintsum = 1000) ## plot traces and posteriors plot(post, "trace") plot(post) ## remove burn-in post <- window(post, start = 5000) ## summarise posteriors summary(post)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.