R/fusion_mixG.R

In rchan26/mixGaussTempering: What the Package Does (One Line, Title Case)

##########################
# these following functions allow user to perform fusion when the sub-posteriors tempered Mixture Gaussian distributions
### the function 'fusion_diff' is a function to implement fusion where the inverse temperatures may differ
# in this case, the parameter 'betas' are a vector of length m, corresponding to the inverse temperatures that we want to fuse
### the function 'fusion_TA' (TA = 'time adaptive') is a function to implement fusion where the inverse temperetures may differ AND time is scaled accordingly
##########################

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

#' Standard Fusion
#'
#' Monte Carlo Fusion for sub-posteriors which are tempered mixture Gaussians  with same weights, means, sds components
#'
#' @param N number of samples
#' @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 weights vector: weights of mixture Gaussian
#' @param means vector: means of mixture Gassuan
#' @param sds vector: st.devs of mixture Gaussian
#' @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
#' # setting variables
#' w_example <- c(0.35, 0.65)
#' m_example <- c(-3, 5)
#' s_example <- c(1, 1.5)
#' b_example <- 1/2
#'
#' # sampling from tempered density
#' nsamples <- 500000
#' base <- hmc_base_sampler_mixG(w_example, m_example, s_example, b_example, nsamples, 2)
#'
#' test <- fusion_diff_mixG(N = 10000, time = 1, C = 2, samples_to_fuse = base,
#'                          weights = w_example, means = m_example, sds = s_example,
#'                          betas = rep(b_example, 2), acceptance_rate = T, timed = T)
#' @export
fusion_diff_mixG <- function(N, time, C, samples_to_fuse, weights, means, sds, 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 the lower bound of the phi_function for the different betas
  K <- sapply(betas, function(beta) lower_bound_phi_function_mixG(weights, means, sds, beta))

  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_mixG(x0 = x[c], y = y, s = 0, t = time,
                                                                            K = K[c], accept_prob = TRUE,
                                                                            weights, means, sds, beta = betas[c])))
      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 which are tempered mixture Gaussians  with same weights, means, sds components
#'
#' @param N number of samples
#' @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 weights vector: weights of mixture Gaussian
#' @param means vector: means of mixture Gassuan
#' @param sds vector: st.devs of mixture Gaussian
#' @param betas vector of length C, where betas[c] is the beta for 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 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
#' # setting variables
#' w_example <- c(0.35, 0.65)
#' m_example <- c(-3, 5)
#' s_example <- c(1, 1.5)
#' b_example <- 1/2
#'
#' # sampling from tempered density
#' nsamples <- 500000
#' base <- hmc_base_sampler_mixG(w_example, m_example, s_example, b_example, nsamples, 2)
#'
#' test <- fusion_TA_mixG(N = 10000, time = 1, C = 2, samples_to_fuse = base,
#'                        weights = w_example, means = m_example, sds = s_example,
#'                        betas = rep(b_example, 2), acceptance_rate = T, timed = T)
#'
#'
#' @export
fusion_TA_mixG <- function(N, time, C, samples_to_fuse, weights, means, sds, 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 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 <- betas
  }

  # finding the lower bound of the phi_function for the different betas
  K <- sapply(betas, function(beta) lower_bound_phi_function_mixG(weights, means, sds, beta))

  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*((x - weighted_avg)^(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_mixG(x0 = x[c], y = y, s = T_star - new_times[c], t = T_star,
                                                                            K = K[c], accept_prob = TRUE,
                                                                            weights, means, sds, beta = betas[c])))
      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)
  }
}