R/pmcmc.R

##' 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 pmcmc,ANY-method pmcmc,missing-method
##' @author Edward L. Ionides, Aaron A. King, Sebastian Funk
##' @seealso \link[=proposals]{MCMC proposals}
##' @family particle filter methods
##' @family \pkg{pomp} parameter estimation 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.
##'
##' @return An object of class \sQuote{pmcmcd_pomp}.
##'
##' @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{tol}, \code{max.fail}, and \code{verbose}, the defaults
##' of which are shown above).  If one does specify additional arguments, these
##' will override the defaults.
##'
##' @inheritSection pfilter Filtering failures
##'
##' @references
##' C. Andrieu, A. Doucet, and R. Holenstein (2010)
##' Particle Markov chain Monte Carlo methods.
##' Journal of the Royal Statistical Society, Series B, 72: 269–342.
##'
##' C. Andrieu and G.O. Roberts (2009)
##' The pseudo-marginal approach for computation
##' Annals of Statistics, 37:697-725.
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("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)
  }
)

##' @name pmcmc-data.frame
##' @aliases pmcmc,data.frame-method
##' @rdname pmcmc
##' @export
setMethod(
  "pmcmc",
  signature=signature(data="data.frame"),
  function (data,
    Nmcmc = 1, proposal,
    Np, tol = 1e-17, max.fail = Inf,
    params, rinit, rprocess, dmeasure, dprior,
    ..., verbose = getOption("verbose", FALSE)) {

    tryCatch(
      pmcmc.internal(
        data,
        Nmcmc=Nmcmc,
        proposal=proposal,
        Np=Np,
        tol=tol,
        max.fail=max.fail,
        params=params,
        rinit=rinit,
        rprocess=rprocess,
        dmeasure=dmeasure,
        dprior=dprior,
        ...,
        verbose=verbose
      ),
      error = function (e) pStop("pmcmc",conditionMessage(e))
    )

  }
)

##' @name pmcmc-pomp
##' @aliases pmcmc,pomp-method
##' @rdname pmcmc
##' @export
setMethod(
  "pmcmc",
  signature=signature(data="pomp"),
  function (data,
    Nmcmc = 1, proposal,
    Np, tol = 1e-17, max.fail = Inf,
    ..., verbose = getOption("verbose", FALSE)) {

    tryCatch(
      pmcmc.internal(
        data,
        Nmcmc=Nmcmc,
        proposal=proposal,
        Np=Np,
        tol=tol,
        max.fail=max.fail,
        ...,
        verbose=verbose
      ),
      error = function (e) pStop("pmcmc",conditionMessage(e))
    )

  }
)

##' @name pmcmc-pfilterd_pomp
##' @aliases pmcmc,pfilterd_pomp-method
##' @rdname pmcmc
##' @export
setMethod(
  "pmcmc",
  signature=signature(data="pfilterd_pomp"),
  function (data,
    Nmcmc = 1, proposal,
    Np, tol, max.fail = Inf,
    ..., verbose = getOption("verbose", FALSE)) {

    if (missing(Np)) Np <- data@Np
    if (missing(tol)) tol <- data@tol

    pmcmc(
      as(data,"pomp"),
      Nmcmc=Nmcmc,
      proposal=proposal,
      Np=Np,
      tol=tol,
      max.fail=max.fail,
      ...,
      verbose=verbose
    )

  }
)

##' @name pmcmc-pmcmcd_pomp
##' @aliases pmcmc,pmcmcd_pomp-method
##' @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
    )

  }
)

##' @name continue-pmcmcd_pomp
##' @aliases continue,pmcmcd_pomp-method
##' @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","variable")
    ft <- array(dim=replace(dim(obj@filter.traj),2L,ndone+Nmcmc),
      dimnames=replace(dimnames(obj@filter.traj),2L,
        list(seq_len(ndone+Nmcmc))))
    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, tol, max.fail, ...,
  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=0),
    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.")

  Nmcmc <- as.integer(Nmcmc)
  if (length(Nmcmc)!=1 || !is.finite(Nmcmc) || Nmcmc < 0)
    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+1,length(theta)+3),
    dimnames=list(
      iteration=seq(from=0,to=Nmcmc,by=1),
      variable=c("loglik","log.prior","nfail",names(theta))))

  if (.ndone==0L) { ## compute prior and likelihood on initial parameter vector

    pfp <- pfilter(object,params=theta,Np=Np,tol=tol,max.fail=max.fail,
      filter.traj=TRUE,.gnsi=gnsi,verbose=verbose)

    log.prior <- dprior(object,params=theta,log=TRUE,.gnsi=gnsi)

    gnsi <- FALSE

  } else { ## has been computed previously

    pfp <- .prev.pfp
    log.prior <- .prev.log.prior
    pfp@filter.traj <- pfp@filter.traj[,.ndone,,drop=FALSE]

  }

  traces[1,names(theta)] <- theta
  traces[1,c(1,2,3)] <- c(pfp@loglik,log.prior,pfp@nfail)

  filt.t <- array(data=numeric(0),dim=replace(dim(pfp@filter.traj),2L,Nmcmc),
    dimnames=replace(dimnames(pfp@filter.traj),2L,
      list(as.character(seq_len(Nmcmc)))))

  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,
        tol=tol,
        max.fail=max.fail,
        filter.traj=TRUE,
        verbose=verbose,
        .gnsi=gnsi
      )

      gnsi <- FALSE

      ## PMCMC update rule (OK because proposal is symmetric)
      alpha <- exp(pfp.prop@loglik+log.prior.prop-pfp@loglik-log.prior)
      if (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(1,2,3)] <- c(pfp@loglik,log.prior,pfp@nfail)

    if (verbose) cat("PMCMC iteration",n+.ndone,"of",Nmcmc+.ndone,
      "completed\nacceptance ratio:",
      round(.accepts/(n+.ndone),3),"\n")

  }

  pars <- apply(traces,2,function(x)diff(range(x))>0)
  pars <- setdiff(names(pars[pars]),c("loglik","log.prior","nfail"))

  new(
    "pmcmcd_pomp",
    pfp,
    params=theta,
    pars=pars,
    Nmcmc=Nmcmc,
    accepts=.accepts,
    proposal=proposal,
    traces=traces,
    log.prior=log.prior,
    filter.traj=filt.t
  )

}
kidusasfaw/pomp documentation built on May 20, 2019, 2:59 p.m.