R/particleMCMC.R

In JMacDonaldPhD/REpi: Pseudo-Marginal Analysis of partially observed Epidemics

#' @name particleMCMC
#' @title Particle MCMC
#' @description 
#' Generates dependent samples from target distribution
#' using a Particle MCMC proposal scheme.
#' @param init initial value of epidemic parameters
#' @param epiModel epidemic model
#' @param obsFrame Generator function for observational model.
#' @param y Observed epidemic data.
#' @param X0 Initial state of the epidemic, which assumed to be known.
#' @param alpha Observational parameters to be passed on to the log-likelihood function
#'              generated by `obsFrame(X_sim)`
#' @param logPrior Functions which calculate log-density of the assumed priors of the 
#'                 epidemic parameters.
#' @param lambda Random Walk Metropolis (RWM) proposal scale parameter
#' @param V A square matrix containing the covariance values for RWM proposal. Shape
#'          of matrix should match the length of the init parameter.
#' @param K  Number of particles to be used in particle filter. This can be optimised using
#'           `adapt_BS_PF()`.
#' @param noIts Number of iterations of MCMC sampler scheme to carry out.
#' @return 
#' A matrix of MCMC samples and proposal acceptance rate of the sample.
#' 
#' @export
particleMCMC <- function(init, epiModel, obsFrame, y, X0, alpha, logPrior, lambda, V, K, noIts){
  start <- as.numeric(Sys.time())
  # Set up functions
  particleFilter <- VBS_PF(y, X0, obs = obsFrame, epiModel = epiModel)
  
  # Estimate Likelihood for initial parameters
  logLikeCurr <- -Inf
  while(is.infinite(logLikeCurr)){
    PF_curr <- particleFilter(K, init, alpha)
    logLikeCurr <- PF_curr$logLikeEst
  }
  curr <- init
  accept <- 0
  k <- length(init)
  draws <- matrix(ncol = k + 1, nrow = noIts)
  
  ESS_store <- matrix(nrow = noIts, ncol = length(PF_curr$ESS))
  for(i in 1:noIts){
    # Propose new parameters
    prop <- abs(curr + mvtnorm::rmvnorm(1, mean = rep(0, k), sigma = (lambda^2)*V))
    
    # Estimate Likelihood
    PF_prop <- particleFilter(K, prop, alpha)
    logLikeProp <- PF_prop$logLikeEst

    if(!is.infinite(logLikeProp)){
      logAccProb <- (logLikeProp + logPrior(prop)) - (logLikeCurr + logPrior(curr)) 
      #print(logAccPRob)
      if(log(runif(1, 0, 1)) < logAccProb){
        curr <- prop
        logLikeCurr <- logLikeProp 
        PF_curr <- PF_prop
        accept <- accept + 1
      }
    }
    ESS_store[i, ] <- PF_curr$ESS
    draws[i, ] <- c(curr, logLikeCurr)
  }
  return(list(draws = draws, acceptRate = accept/noIts,
              curr = curr, ESS_store = ESS_store, time = as.numeric(Sys.time()) - start))
}