R/rejection_sampler_exp_4.R

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

#' Target density exp(-((x-mean)^4)/2)
#'
#' Normalised target density
#'
#' @param x real value
#' @param mean mean value, defaults to 0
#'
#' @return real value
#'
#' @examples
#' curve(target_density_exp_4(x), from = -4, to = 4)
#'
#' @export
target_density_exp_4 <- function(x, mean = 0) {
  norm_constant <- integrate(function(y) exp(-((y-mean)^4)/2), lower = -Inf, upper = Inf)$value
  return(exp(-((x-mean)^4)/2) / norm_constant)
}

#' Tempered target density exp(-(((x-mean)^4)*beta)/2)
#'
#' Normalised tempered target density
#'
#' @param x real value
#' @param mean mean value, defaults to 0
#' @param beta real value between 0 and 1
#'
#' @return real value
#'
#' @examples
#' curve(tempered_target_density_exp_4(x, 1/2), from = -4, to = 4)
#'
#' @export
tempered_target_density_exp_4 <- function(x, mean = 0, beta) {
  norm_constant <- integrate(function(y) exp(-(((y-mean)^4)*beta)/2), lower = -Inf, upper = Inf)$value
  return(exp(-(((x-mean)^4)*beta)/2) / norm_constant)
}

#' Rejection sampler for tempered target exp(-(beta*(x-mean)^4)/2)
#'
#' Sample from tempered target using rejection sampling with normal proposals
#'
#' @param N number of samples
#' @param mean mean value, defaults to 0
#' @param proposal_mean mean of proposal density for rejection sampler
#' @param proposal_sd st.dev of proposal density for rejection sampler
#' @param dominiating_M constant M to bound the target density
#' @param beta real value between 0 and 1
#'
#' @return samples from tempered target density
#'
#' @examples
#' dominating_dnorm <- function(x) 1.3*dnorm(x, mean = 0, sd = 1)
#' curve(dominating_dnorm(x), -5, 5, col = 'red')
#' curve(dnorm(x, mean = 0, sd = 1), -5, 5, add = T, col = 'green')
#' curve(tempered_target_density_exp_4(x, 1/2), -5, 5, add = T, col = 'blue')
#' sample_from_fc(N = 100000, proposal_mean = 0, proposal_sd = 1, dominating_M = 1.3)
#'
#' @export
sample_from_fc <- function(N, mean = 0, proposal_mean, proposal_sd, dominating_M, beta) {
  # obtain samples from f_{c} by using a Standard Normal proposal density (g)
  # have to pass through the mean and sd for the proposal density with proposal_mean and proposal_sd
  # dominating_M is the constant M so that f(x) / M*g(x) for all x
  i <- 0; samples <- 1:N; iters <- 0

  dominating_function <- function(x) dominating_M*dnorm(x, mean = proposal_mean, sd = proposal_sd)
  target_fc <- function(x) exp(-(((x-mean)^4)*beta)/2) / integrate(function(y) exp(-(((y-mean)^4)*beta)/2), lower = -Inf, upper = Inf)$value

  while (i < N) {
    iters <- iters+1
    prop <- rnorm(1, mean = proposal_mean, sd = proposal_sd)
    if (runif(1,0,1) < (target_fc(prop) / dominating_function(prop))) {
      i <- i+1
      samples[i] <- prop
    }
  }

  # print out the acceptance rate
  print(paste('acceptance rate:', N, '/', iters, '=', N/iters))
  # return samples
  return(samples)
}

#' Rejection sampler for base level
#'
#' Sample for base level (samples of exp(-(beta*(x-mean)^4)/2))
#'
#' @param beta real value between 0 and 1
#' @param mean mean value, defaults to 0
#' @param nsamples number of samples per node
#' @param proposal_mean mean of proposal density for rejection sampler
#' @param proposal_sd st.dev of proposal density for rejection sampler
#' @param dominiating_M constant M to bound the target density
#'
#' @return return_name return description
#'
#' @examples
#' base_rejection_sampler_exp_4(beta = 1/2, nsamples = 100000, proposal_mean = 0, proposal_sd = 1, dominating_M = 1.3)
#'
#' @export
base_rejection_sampler_exp_4 <- function(beta, mean = 0, nsamples, proposal_mean, proposal_sd, dominating_M) {
  print(paste('sampling from tempered target density with beta =', 1, '/', (1/beta)))
  # rejection sampler for tempered target density
  base <- sample_from_fc(N = (1/beta)*nsamples,
                         mean = mean,
                         proposal_mean = proposal_mean,
                         proposal_sd = proposal_sd,
                         dominating_M = dominating_M,
                         beta = beta)
  base_samples <- split(base, ceiling((1:length(base))/nsamples))
  return(base_samples)
}