##' The particle Markov chain Metropolis-Hastings algorithm
##'
##' The Particle MCMC algorithm for estimating the parameters of a
##' partially-observed Markov process. Running \code{pmcmc} causes a particle
##' random-walk Metropolis-Hastings Markov chain algorithm to run for the
##' specified number of proposals.
##'
##' @name pmcmc
##' @rdname pmcmc
##' @include pfilter.R proposals.R load.R continue.R
##' @aliases pmcmc,ANY-method pmcmc,missing-method
##' @author Edward L. Ionides, Aaron A. King, Sebastian Funk
##' @family estimation methods
##' @family particle filter methods
##' @family full-information methods
##' @family MCMC methods
##' @family Bayesian methods
##'
##' @importFrom stats runif
##' @inheritParams pomp
##' @inheritParams pfilter
##' @param Nmcmc The number of PMCMC iterations to perform.
##' @param proposal optional function that draws from the proposal
##' distribution. Currently, the proposal distribution must be symmetric for
##' proper inference: it is the user's responsibility to ensure that it is.
##' Several functions that construct appropriate proposal function are
##' provided: see \link[=proposals]{MCMC proposals} for more information.
##'
##' @section Methods:
##' The following can be applied to the output of a \code{pmcmc} operation:
##' \describe{
##' \item{\code{pmcmc}}{repeats the calculation, beginning with the last state}
##' \item{\code{\link{continue}}}{continues the \code{pmcmc} calculation}
##' \item{\code{plot}}{produces a series of diagnostic plots}
##' \item{\code{\link{filter_traj}}}{extracts a random sample from the smoothing distribution}
##' \item{\code{\link{traces}}}{produces an \code{\link[coda]{mcmc}} object, to which the various \pkg{coda} convergence diagnostics can be applied}
##' }
##'
##' @return An object of class \sQuote{pmcmcd_pomp}.
##'
##' @inheritSection pomp Note for Windows users
##'
##' @section Re-running PMCMC Iterations:
##' To re-run a sequence of PMCMC
##' iterations, one can use the \code{pmcmc} method on a \sQuote{pmcmc} object.
##' By default, the same parameters used for the original PMCMC run are re-used
##' (except for \code{verbose}, the default of which is shown above). If one
##' does specify additional arguments, these will override the defaults.
##'
##' @references
##'
##' \Andrieu2010
##'
NULL
setClass(
"pmcmcd_pomp",
contains="pfilterd_pomp",
slots=c(
pars = "character",
Nmcmc = "integer",
accepts = "integer",
proposal = "function",
traces = "array",
log.prior = "numeric"
),
prototype=prototype(
pars = character(0),
Nmcmc = 0L,
accepts = 0L,
proposal = function (...)
pStop(who="pmcmc","proposal not specified."),
traces=array(dim=c(0,0)),
log.prior=numeric(0)
)
)
setGeneric(
"pmcmc",
function (data, ...)
standardGeneric("pmcmc")
)
setMethod(
"pmcmc",
signature=signature(data="missing"),
definition=function (...) {
reqd_arg("pmcmc","data")
}
)
setMethod(
"pmcmc",
signature=signature(data="ANY"),
definition=function (data, ...) {
undef_method("pmcmc",data)
}
)
##' @rdname pmcmc
##' @export
setMethod(
"pmcmc",
signature=signature(data="data.frame"),
function (data,
Nmcmc = 1, proposal,
Np,
params, rinit, rprocess, dmeasure, dprior,
..., verbose = getOption("verbose", FALSE)) {
tryCatch(
pmcmc_internal(
data,
Nmcmc=Nmcmc,
proposal=proposal,
Np=Np,
params=params,
rinit=rinit,
rprocess=rprocess,
dmeasure=dmeasure,
dprior=dprior,
...,
verbose=verbose
),
error = function (e) pStop(who="pmcmc",conditionMessage(e))
)
}
)
##' @rdname pmcmc
##' @export
setMethod(
"pmcmc",
signature=signature(data="pomp"),
function (data,
Nmcmc = 1, proposal,
Np, ...,
verbose = getOption("verbose", FALSE)) {
tryCatch(
pmcmc_internal(
data,
Nmcmc=Nmcmc,
proposal=proposal,
Np=Np,
...,
verbose=verbose
),
error = function (e) pStop(who="pmcmc",conditionMessage(e))
)
}
)
##' @rdname pmcmc
##' @export
setMethod(
"pmcmc",
signature=signature(data="pfilterd_pomp"),
function (data,
Nmcmc = 1, proposal,
Np, ...,
verbose = getOption("verbose", FALSE)) {
if (missing(Np)) Np <- data@Np
pmcmc(
as(data,"pomp"),
Nmcmc=Nmcmc,
proposal=proposal,
Np=Np,
...,
verbose=verbose
)
}
)
##' @rdname pmcmc
##' @export
setMethod(
"pmcmc",
signature=signature(data="pmcmcd_pomp"),
function (data,
Nmcmc, proposal,
..., verbose = getOption("verbose", FALSE)) {
if (missing(Nmcmc)) Nmcmc <- data@Nmcmc
if (missing(proposal)) proposal <- data@proposal
pmcmc(
as(data,"pfilterd_pomp"),
Nmcmc=Nmcmc,
proposal=proposal,
...,
verbose=verbose
)
}
)
##' @rdname continue
##' @param Nmcmc positive integer; number of additional PMCMC iterations to perform
##' @export
setMethod(
"continue",
signature=signature(object="pmcmcd_pomp"),
function (object, Nmcmc = 1, ...) {
ndone <- object@Nmcmc
accepts <- object@accepts
obj <- pmcmc(
object,
Nmcmc=Nmcmc,
...,
.ndone=ndone,
.accepts=accepts,
.prev.pfp=as(object,"pfilterd_pomp"),
.prev.log.prior=object@log.prior
)
obj@traces <- rbind(
object@traces[,colnames(obj@traces)],
obj@traces[-1,]
)
names(dimnames(obj@traces)) <- c("iteration","name")
ft <- array(
dim=replace(dim(obj@filter.traj),2L,ndone+Nmcmc),
dimnames=list(
name=rownames(obj@filter.traj),
rep=NULL,
time=NULL
)
)
ft[,seq_len(ndone),] <- object@filter.traj
ft[,ndone+seq_len(Nmcmc),] <- obj@filter.traj
obj@filter.traj <- ft
obj@Nmcmc <- as.integer(ndone+Nmcmc)
obj@accepts <- as.integer(accepts+obj@accepts)
obj
}
)
pmcmc_internal <- function (object, Nmcmc, proposal, Np, ...,
verbose, .ndone = 0L, .accepts = 0L, .prev.pfp = NULL, .prev.log.prior = NULL,
.gnsi = TRUE) {
verbose <- as.logical(verbose)
object <- pomp(object,...,verbose=verbose)
gnsi <- as.logical(.gnsi)
.ndone <- as.integer(.ndone)
.accepts <- as.integer(.accepts)
if (missing(proposal) || is.null(proposal))
pStop_(sQuote("proposal")," must be specified")
proposal <- tryCatch(
match.fun(proposal),
error = function (e) {
pStop_(sQuote("proposal")," must be a function: ",conditionMessage(e))
}
)
pompLoad(object,verbose=verbose)
on.exit(pompUnload(object,verbose=verbose))
start <- coef(object)
## test proposal distribution
theta <- tryCatch(
proposal(start,.n=0L),
error = function (e) {
pStop_("error in proposal function: ",conditionMessage(e))
}
)
if (is.null(names(theta)) || !is.numeric(theta) || any(names(theta)==""))
pStop_(sQuote("proposal")," must return a named numeric vector.")
theta <- start
Nmcmc <- as.integer(Nmcmc)
if (length(Nmcmc)!=1L || !is.finite(Nmcmc) || Nmcmc < 0L)
pStop_(sQuote("Nmcmc")," must be a positive integer")
if (verbose)
cat("performing",Nmcmc,"PMCMC iteration(s) using",Np[1L],"particles\n")
traces <- array(
data=NA_real_,
dim=c(Nmcmc+1L,length(theta)+2L),
dimnames=list(
iteration=NULL,
name=c("loglik","log.prior",names(theta))
)
)
if (.ndone==0L) { ## compute prior and likelihood on initial parameter vector
log.prior <- dprior(object,params=theta,log=TRUE,.gnsi=gnsi)
if (!is.finite(log.prior))
pStop_("non-finite log prior at starting parameters.")
pfp <- pfilter(
object,
params=theta,
Np=Np,
filter.traj=TRUE,
.gnsi=gnsi,
verbose=verbose
)
if (!is.finite(pfp@loglik))
pStop_("non-finite log likelihood at starting parameters.")
} else { ## has been computed previously
pfp <- .prev.pfp
log.prior <- .prev.log.prior
pfp@filter.traj <- pfp@filter.traj[,.ndone,,drop=FALSE]
}
traces[1L,names(theta)] <- theta
traces[1L,c(1L,2L)] <- c(pfp@loglik,log.prior)
filt.t <- array(
data=NA_real_,
dim=replace(dim(pfp@filter.traj),2L,Nmcmc),
dimnames=list(
name=rownames(pfp@filter.traj),
rep=NULL,
time=NULL
)
)
for (n in seq_len(Nmcmc)) { # main loop
theta.prop <- tryCatch(
proposal(theta,.n=n+.ndone,.accepts=.accepts,verbose=verbose),
error = function (e) {
pStop_("error in proposal function: ",conditionMessage(e))
}
)
## compute log prior
log.prior.prop <- dprior(object,params=theta.prop,log=TRUE,.gnsi=gnsi)
if (is.finite(log.prior.prop)) {
## run the particle filter on the proposed new parameter values
pfp.prop <- pfilter(
pfp,
params=theta.prop,
Np=Np,
filter.traj=TRUE,
verbose=verbose,
.gnsi=gnsi
)
## PMCMC update rule (OK because proposal is symmetric)
alpha <- exp(pfp.prop@loglik+log.prior.prop-pfp@loglik-log.prior)
if (!is.finite(alpha))
pWarn_("non-finite log likelihood or log prior encountered.") #nocov
if (is.finite(alpha) && runif(1) < alpha) {
pfp <- pfp.prop
theta <- theta.prop
log.prior <- log.prior.prop
.accepts <- .accepts+1L
}
}
## add filtered trajectory to the store
filt.t[,n,] <- pfp@filter.traj[,1L,]
## store a record of this iteration
traces[n+1,names(theta)] <- theta
traces[n+1,c(1L,2L)] <- c(pfp@loglik,log.prior)
if (verbose) cat("PMCMC iteration",n+.ndone,"of",Nmcmc+.ndone,
"completed\nacceptance ratio:",
round(.accepts/(n+.ndone),3),"\n")
gnsi <- FALSE
}
pars <- apply(traces,2L,\(x)diff(range(x))>0)
pars <- setdiff(names(pars[pars]),c("loglik","log.prior"))
new(
"pmcmcd_pomp",
pfp,
params=theta,
pars=pars,
Nmcmc=Nmcmc,
accepts=.accepts,
proposal=proposal,
traces=traces,
log.prior=log.prior,
filter.traj=filt.t
)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.