R/hmc_exp_4.R

In rchan26/exp4Tempering: Hierarchical and Sequential Monte Carlo Fusion for exp(-(x^4)/2)

##########################
# function to use HMC (STAN) to sample from the tempered Gaussian density
# stan file 'mixture_gaussian.stan' must be in the same directory when trying to use the function
##########################

#' HMC sampler for tempered target exp(-(beta*(x-mean)^4)/2))
#'
#' Sample from tempered target using Stan
#'
#' @param beta temperature level
#' @param mean mean value, defaults to 0
#' @param iterations number of iterations per chain - note that burn in is 0.5 of this number
#' @param chains number of chains
#' @param output boolean value: defaults to T, determines whether or not to print output to console
#'
#' @return samples from the tempered target
#'
#' @examples
#' hmc_sample_tempered_exp_mean_4(beta = 1, iterations = 2*100000, chains = 1, output = F)
#'
#' @export
hmc_sample_tempered_exp_4 <- function(beta, mean = 0, iterations, chains, output = F) {
  # print output to console
  print('Sampling from tempered target density')

  # function to use Stan (HMC) to sample from the tempered mixture Gaussian distribution
  training_data <- list(beta = beta, mu = mean)

  if (output) {
    # print output to console
    model <- rstan::sampling(object = stanmodels$exp_4,
                             data = training_data,
                             iter = iterations,
                             chains = chains,
                             control = list(adapt_delta = 0.99, max_treedepth = 15))
  } else {
    # hide output from console
    model <- rstan::sampling(object = stanmodels$exp_4,
                             data = training_data,
                             iter = iterations,
                             chains = chains,
                             control = list(adapt_delta = 0.99, max_treedepth = 15),
                             verbose = FALSE,
                             refresh = 0)
  }

  # print completion
  print('Finished sampling from tempered target density')
  return(rstan::extract(model)$x)
}

#' HMC sampler for base level
#'
#' Sample for base level (samples of exp(-(beta*(x-mean)^4)/2))
#'
#' @param beta temperature level
#' @param mean mean value
#' @param nsamples number of samples per node
#' @param nchains number of nodes
#'
#' @return samples from tempered target
#'
#' @examples
#' hmc_base_sampler_exp_mean_4(beta = 1/2, nsamples = 100000, chains = 2, output = F)
#'
#' @export
hmc_base_sampler_exp_4 <- function(beta, mean = 0, nsamples, nchains) {
  base <- hmc_sample_tempered_exp_mean_4(beta = beta, mean = mean, iterations = 2*nsamples, chains = nchains, output = F)
  base_samples <- split(base, ceiling((1:length(base))/nsamples))
  return(base_samples)
}