R/fusion_exp_4.R

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

######################################## fusion for when fusing samples from different betas ########################################

#' Standard Fusion
#'
#' Monte Carlo Fusion for sub-posteriors of the form exp(-((x^4)*beta)/2)
#'
#' @param N number of samples
#' @param mean mean value
#' @param time time T for fusion algorithm
#' @param C number of sub-posteriors to combine
#' @param samples_to_fuse list of length C, where samples_to_fuse[c] containg the samples for the c-th sub-posterior
#' @param betas vector of length C, where betas[c] is the beta for c-th sub-posterior
#' @param level defaults to 1, used in hierarchical and sequential Monte Carlo Fusion
#' @param acceptance_rate boolean value: defaults to F, determines whether or not to return acceptance rates
#' @param timed boolen value: defaults to T, determines whether or not to return the time elapsed to run
#'
#' @return samples: fusion samples
#' @return iters_rho: number of iterations from the first accept/reject step (rho)
#' @return iters_Q: number of iterations from the second (diffusion) accept/reject step (Q)
#' @return time: run-time of fusion sampler
#'
#' @examples
#' input_samples <- base_rejection_sampler_exp_4(beta = 1/2, nsamples = 100000, proposal_mean = 0, proposal_sd = 1, dominating_M = 1.3)
#' test1_standard <- fusion_diff_exp_4(N = 10000, time = 1, C = 2, samples_to_fuse = input_samples, betas = rep(1/2, 2), timed = T)
#'
#' @export
fusion_diff_exp_4 <- function(N, mean = 0, time, C, samples_to_fuse, betas, level = 1, acceptance_rate = F, timed = F) {
  # obtain samples from fusion density by using Fusion
  if (is.list(samples_to_fuse) && length(samples_to_fuse)!=C) stop('length of samples_to_fuse must be equal to C')
  if (!is.list(samples_to_fuse)) stop('samples_to_fuse must be a list of length C')
  if (length(betas)!=C) {
    if (length(betas)==1) {
      # if only one is given, then the betas are the same
      betas <- rep(betas, C)
    } else {
      stop('betas must be a vector of length 1 (the same for all) or C (different)')
    }
  }

  # finding lower bound of the phi_function for given betas
  K <- sapply(betas, function(beta) optimise(function(x) phi_function_exp_4(x, mean = mean, beta = beta), interval = c(0, 100), maximum = FALSE)$objective)

  cat('Starting fusion for level', level, 'trying to get', N, 'samples\n', file = 'fusion_diff_progress.txt', append = T)
  f_samples <- rep(0, N); i <- 0
  iters_rho <- 0; iters_Q <- 0

  pcm <- proc.time()
  while (i < N) {
    # logging the number of iterations for rho probability
    iters_rho <- iters_rho+1
    # sampling from each of the C components
    x <- sapply(samples_to_fuse, function(node) sample(node, 1))
    # generate standard uniform random variable and calculating acceptance probability
    logrho <- -sum((x - mean(x))^2) / (2*time)

    if (log(runif(1,0,1)) < logrho) {
      # logging the number of iterations for the Q probability
      iters_Q <- iters_Q+1
      # simulate proposed value y from a Gaussian distribution
      y <- rnorm(1, mean = mean(x), sd = sqrt(time/C))
      # simulate C diffusions
      prob <- prod(sapply(1:C, function(c) simulate_langevin_diffusion_exp_4(x0 = x[c], y = y, s = 0, t = time, K = K[c],
                                                                             mean = mean, beta = betas[c], accept_prob = TRUE)))
      if (runif(1,0,1) < prob) {
        i <- i+1
        f_samples[i] <- y
        cat('Level:', level, '|| Iteration:', i, '/', N, '\n', file = 'fusion_diff_progress.txt', append = T)
      }
    }
  }
  final <- proc.time() - pcm

  # print completion
  cat('Completed fusion for level', level, '\n', file = 'fusion_diff_progress.txt', append = T)

  if (acceptance_rate && timed) {
    return(list('samples' = f_samples, 'iters_rho' = iters_rho, 'iters_Q' = iters_Q, 'time' = final['elapsed']))
  } else if (acceptance_rate && !timed) {
    return(list('samples' = f_samples, 'iters_rho' = iters_rho, 'iters_Q' = iters_Q))
  } else if (!acceptance_rate && timed) {
    return(list('samples' = f_samples, 'time' = final['elapsed']))
  } else {
    return(f_samples)
  }
}



######################################## fusion for when fusing samples from different betas and different times ########################################

#' Time-adapting Fusion
#'
#' Time-adapting Monte Carlo Fusion for sub-posteriors of the form exp(-((x^4)*beta)/2)
#'
#' @param N number of samples
#' @param mean mean value
#' @param time time T for time-adapting fusion algorithm
#' @param C number of sub-posteriors to combine
#' @param samples_to_fuse list of length C, where samples_to_fuse[c] containg the samples for the c-th sub-posterior
#' @param betas vector of length C, where betas[c] is the beta for c-th sub-posterior
#' @param sample_weights defaults to sqrt(betas) if missing: vector of length C, where sample_weights[c] is the sample wieght for the c-th sub-posterior
#' @param level defaults to 1, used in hierarchical and sequential Monte Carlo Fusion
#' @param acceptance_rate boolean value: defaults to F, determines whether or not to return acceptance rates
#' @param timed boolen value: defaults to T, determines whether or not to return the time elapsed to run
#'
#' @return samples: fusion samples
#' @return iters_rho: number of iterations from the first accept/reject step (rho)
#' @return iters_Q: number of iterations from the second (diffusion) accept/reject step (Q)
#' @return time: run-time of fusion sampler
#'
#' @examples
#' input_samples <- base_rejection_sampler_exp_4(beta = 1/2, nsamples = 100000, proposal_mean = 0, proposal_sd = 1, dominating_M = 1.3)
#' test1_standard <- fusion_diff_exp_4(N = 10000, time = 1, C = 2, samples_to_fuse = input_samples, betas = rep(1/2, 2), timed = T)
#'
#' @export
fusion_TA_exp_4 <- function(N, mean = 0, time, C, samples_to_fuse, betas, sample_weights, level = 1, acceptance_rate = F, timed = F) {
  # obtain samples from fusion density by using Fusion
  ### but now the time is scaled by the sqrt(betas) (i.e. the exponent) of the sub-posterior
  if (is.list(samples_to_fuse) && length(samples_to_fuse)!=C) stop('length of samples_to_fuse must be equal to C')
  if (!is.list(samples_to_fuse)) stop('samples_to_fuse must be a list of length C')
  if (length(betas)!=C) {
    if (length(betas)==1) {
      # if only one is given, then the betas are the same
      betas <- rep(betas, C)
    } else {
      stop('betas must be a vector of length 1 (the same for all) or C (different)')
    }
  }
  if (missing(sample_weights)) {
    sample_weights <- sqrt(betas)
  }
  if (length(sample_weights)==1) {
    # if only one is given, then the sample weights are the same
    sample_weights <- rep(sample_weights, C)
  }

  # finding lower bound of the phi_function for given betas
  K <- sapply(betas, function(beta) optimise(function(x) phi_function_exp_4(x, mean = mean, beta = beta), interval = c(0, 100), maximum = FALSE)$objective)

  cat('Starting fusion for level', level, 'trying to get', N, 'samples\n', file = 'fusion_TA_progress.txt', append = T)
  f_samples <- rep(0, N); i <- 0
  iters_rho <- 0; iters_Q <- 0

  pcm <- proc.time()
  while (i < N) {
    # logging the number of iterations for rho probability
    iters_rho <- iters_rho+1
    # sampling from each of the C components
    x <- sapply(samples_to_fuse, function(node) sample(node, 1))
    # calculating the weighted average of the samples
    weighted_avg <- weighted.mean(x, w = sample_weights)
    # generate standard uniform random variable and calculating acceptance probability
    logrho <- (-sum(sample_weights*((weighted_avg-x)^(2)))) / (2*time)

    if (log(runif(1,0,1)) < logrho) {
      # logging the number of iterations for the Q probability
      iters_Q <- iters_Q+1
      # finding T_{c}'s and storing them in new_times
      new_times <- time / sample_weights
      # setting T^{*} ('T-star') as the largest T_{c} = (time / w_{c}) for c = 1, ..., C
      T_star <- max(new_times)
      # simulate proposed value y from a Gaussian distribution
      y <- rnorm(1, mean = weighted_avg, sd = sqrt(time/sum(sample_weights)))
      # simulate each of the C diffusions
      prob <- prod(sapply(1:C, function(c) simulate_langevin_diffusion_exp_4(x0 = x[c], y = y, s = T_star - new_times[c],
                                                                             t = T_star, K = K[c], mean = mean, beta = betas[c],
                                                                             accept_prob = TRUE)))
      if (runif(1,0,1) < prob) {
        i <- i+1
        f_samples[i] <- y
        cat('Level:', level, '|| Iteration:', i, '/', N, '\n', file = 'fusion_TA_progress.txt', append = T)
      }
    }
  }
  final <- proc.time() - pcm

  # print completion
  cat('Completed fusion for level', level, '\n', file = 'fusion_TA_progress.txt', append = T)

  if (acceptance_rate && timed) {
    return(list('samples' = f_samples, 'iters_rho' = iters_rho, 'iters_Q' = iters_Q, 'time' = final['elapsed']))
  } else if (acceptance_rate && !timed) {
    return(list('samples' = f_samples, 'iters_rho' = iters_rho, 'iters_Q' = iters_Q))
  } else if (!acceptance_rate && timed) {
    return(list('samples' = f_samples, 'time' = final['elapsed']))
  } else {
    return(f_samples)
  }
}