R/exact_algorithm_exp_4.R

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

######################################## 'phi function' for the Exact Algorithm for tempered target density ########################################

#' 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)
}

#' phi-function for exp(-(((x-mean)^4)*beta)/2)
#'
#' phi-function for the Exact Algorithm for exp(-(((x-mean)^4)*beta)/2)
#'
#' @param x real value
#' @param mean mean value, defaults to 0
#' @param beta real value between 0 and 1
#'
#' @return real value
#'
#' @examples
#' curve(phi_function_exp_4(x, 1), from = -4, to = 4)
#'
#' @export
phi_function_exp_4 <- function(x, mean = 0, beta) {
  return(2*(beta^2)*((x-mean)^6) - 3*beta*((x-mean)^2))
}

#' Obtain bounds for phi function
#'
#' Finds the lower and upper bounds of the phi function between an interval
#'
#' @param beta real value between 0 and 1
#' @param mean mean value, defaults to 0
#' @param lower lower end of interval
#' @param upper upper end of interval
#'
#' @return vector of lower and upper bound of phi-function between [lower, upper]
#'
#' @examples
#' bounds_phi_function_exp_4(beta = 1/2, interval = c(-2, -1))
#'
#' @export
bounds_phi_function_exp_4 <- function(beta, mean = 0, lower, upper) {
  # finding maxima and minimma in the interval
  maxi <- optimise(function(x) phi_function_exp_4(x, mean, beta), interval = c(lower, upper), maximum = TRUE)$objective
  mini <- optimise(function(x) phi_function_exp_4(x, mean, beta), interval = c(lower, upper), maximum = FALSE)$objective

  # checking end points
  if (max(phi_function_exp_4(lower, mean, beta), phi_function_exp_4(upper, mean, beta)) > maxi) {
    up_bound <- max(phi_function_exp_4(lower, mean, beta), phi_function_exp_4(upper, mean, beta))
  } else {
    up_bound <- maxi
  }

  # checking end points
  if (min(phi_function_exp_4(lower, mean, beta), phi_function_exp_4(upper, mean, beta)) < mini) {
    low_bound <- min(phi_function_exp_4(lower, mean, beta), phi_function_exp_4(upper, mean, beta))
  } else {
    low_bound <- mini
  }

  return(list('fun_ubound' = up_bound, 'fun_lbound' = low_bound))
}

#' 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) {
  # 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)
}

######################################## exact algorithm to simulate langevin diffusion ########################################

#' Exact Algorithm for langevin diffusion with pi = exp(-(beta*(x-mean)^4)/2)
#'
#' Simulate langevin diffusion using the Exact Algorithm with pi = exp(-(beta*(x-mean)^4)/2)
#'
#' @param x0 start value
#' @param y end value
#' @param s start time
#' @param t end time
#' @param K lower bound of the phi function
#' @param mean mean value
#' @param beta real value between 0 and 1
#' @param accept_prob boolean value: defaults to F, determined whether or not to return acceptance prob (T) or a skeleton of the diffusion (F)
#'
#' @return accept_prob real value if accept_prob = TRUE
#' @return layered_bb matrix holding skeleton of the diffusion
#'
#' @examples
#' lower_bound <- optimise(function(x) phi_function_exp_4(x, beta = 1/2), interval = c(0, 100), maximum = FALSE)$objective
#' test <- simulate_langevin_diffusion_exp_4(x0 = 0, y = 0.3, s = 0, t = 1, K = lower_bound, beta = 1/2, accept_prob = FALSE)
#' plot(x = test['time',], y = test['X',], type = 'l')
#'
#' @export
simulate_langevin_diffusion_exp_4 <- function (x0, y, s, t, K, mean = 0, beta, accept_prob = FALSE) {
  # simulates langevin diffusion for tempered target
  # if accept_prob = TRUE, function returns the acceptance probability only

  repeat {
    # first simulate layer information to bound sample path
    a_seq <- seq(from = 0, to = ceiling(t-s), by = 0.5)
    bes_layer <- layeredBB::bessel_layer_simulation(x0, y, s, t, a_seq)
    lbound <- min(x0,y) - bes_layer$a[bes_layer$l]
    ubound <- max(x0,y) + bes_layer$a[bes_layer$l]

    # calculate upper and lower bounds of phi given the simulated sample layer information
    bounds <- bounds_phi_function_exp_4(beta, mean, lower = lbound, upper = ubound)
    LX <- bounds$fun_lbound
    UX <- bounds$fun_ubound

    if (!(accept_prob)) {
      if (runif(1,0,1) < (1 - exp(-(LX-K)*(t-s)))) {
        # go back to step 1: start again
        next
      }
    }

    # simulating skeletal points via Poisson thinning
    # simulating number of points from Poisson distribution
    kap <- rpois(1, (UX - LX)*(t-s))
    # simulating Poisson times from a Uniform(0, (t-s)) distribution
    pois_times <- sort(runif(kap, s, t))

    # simulating sample path at Poisson times
    layered_bb <- layeredBB::layered_brownian_bridge(x0, y, s, t, bes_layer$a, bes_layer$l, pois_times)

    # layered_bb includes points for times s,t,tau, so only keeping those from pois_times
    pois_points <- layered_bb[1, (layered_bb[2,] %in% pois_times)]

    # calculating acceptance probability
    acc_prob <- ifelse(test = length(pois_times)==0,
                       yes = 1,
                       no = prod((UX - phi_function_exp_4(x = pois_points, mean, beta = beta))/(UX - LX)))

    if (accept_prob) {
      # the full acceptance probability is acc_prob multiplied by the first acceptance prob
      return(acc_prob * (exp(-(LX-K)*(t-s))))
    }

    if (runif(1,0,1) < acc_prob) {
      return(layered_bb) # accept sample path and return the entire path
    }
  }
}