R/PsMCMC.R

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

#' @name PsMCMC
#' @title Pseudo-Marginal MCMC
#' @description 
#' Generates dependent samples from target distribution
#' using a Pseudo-Marginal MCMC proposal scheme.
#' @param init initial value of epidemic parameters
#' @param epiModel epidemic model
#' @param obsFrame Generator function for observational model.
#' @param epiSample Observed epidemic data.
#' @param I0 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 N Number of particles to be used in Importance Sample Estimate of posterior. This
#'          can be optimised using `adapt_IS()`
#' @param noIts Number of iterations of MCMC sampler scheme to carry out.
#' @return 
#' @return 
#' A matrix of MCMC samples and proposal acceptance rate of the sample.
#' 
#' @export
PsMCMC <- function(init, epiModel, obsFrame, epiSample, I0, alpha, logPrior, lambda, V, N, noIts){
  # Set up functions
  IS_estimator <- ImportanceSampler(epiSample, epiModel, obsFrame) # Likelihood estimator
  k <- length(init)
  # Estimate Likelihood for initial parameters
  logLikeCurr <- -Inf
  while(is.infinite(logLikeCurr)){
    logLikeCurr <- IS_estimator(N, list(I0, init, alpha), func = log_mean_exp)$estimates
  }
  curr <- init
  accept <- 0
  draws <- matrix(ncol = k + 1, nrow = noIts)
  for(i in 1:noIts){
    # Propose new parameters
    prop <- abs(curr + mvnfast::rmvn(1, mu = rep(0, k), sigma = lambda*V))
    
    # Estimate Likelihood
    logLikeProp <- IS_estimator(N, list(I0, prop, alpha), log_mean_exp)$estimates
    
    if(!is.infinite(logLikeProp)){
      logAccProb <- (logLikeProp + logPrior(prop)) - (logLikeCurr + logPrior(curr)) 
      #print(logAccPRob)
      if(log(runif(1, 0, 1)) < logAccProb){
        curr <- prop
        logLikeCurr <- logLikeProp 
        accept <- accept + 1
      }
    }
    draws[i, ] <- c(curr, logLikeCurr)
  }
  return(list(draws = draws, acceptRate = accept/noIts))
}