scripts/fusion_examples.R

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

library(exp4Tempering)
library(fusionAnalysis)

CMC <- function(N) {
  x <- matrix(sample_from_fc(N=4*N, proposal_mean = 0, proposal_sd = 1.5, dominating_M = 1.4, 1/4), ncol = 4)
  weights <- apply(x, 2, var)
  samples <- sapply(1:N, function(i) weighted.mean(x[i,], weights))

  return(samples)
}

######################################## examples ########################################

set.seed(41)

# rejection sampling target
dominating_dnorm <- function(x) 1.35*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(target_density_exp_4(x), -5, 5, add = T, col = 'blue')
test_target_mc <- sample_from_fc(N = 10000, proposal_mean = 0, proposal_sd = 1, dominating_M = 1.35, beta = 1)



######################################## fusion example for fusing pi^(1/2) twice to get pi

# sampling from the tempered targets
# finding the right proposal distribution and constant M in rejection sampling
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')

# using rejection sampling to obtain input samples
input_samples <- base_rejection_sampler_exp_4(beta = 1/2, nsamples = 100000, proposal_mean = 0, proposal_sd = 1, dominating_M = 1.3)
curve(tempered_target_density_exp_4(x, beta = 1/2), -4, 4)
# check the samples look okay
for (samples in input_samples) {
  lines(density(samples), col = 'blue')
}

# normal fusion
test1_standard <- fusion_diff_exp_4(N = 10000, time = 1, C = 2, samples_to_fuse = input_samples, betas = rep(1/2, 2), timed = T)
# time-adapting fusion
test1_TA <- fusion_TA_exp_4(N = 10000, time = 1, C = 2, samples_to_fuse = input_samples, betas = rep(1/2, 2),
                            sample_weights = 1, timed = T)

# plot kdes for fusion samples
curve(target_density_exp_4(x), -2.5, 2.5)
lines(density(test_target_mc), col = 'black')
lines(density(test1_standard$samples), col = 'green')
lines(density(test1_TA$samples), col = 'blue')



######################################## fusion example for fusing pi^(1/4) four times to get pi

# sampling from the tempered targets
input_samples <- hmc_base_sampler_exp_4(beta = 1/4, nsamples = 100000, nchains = 4)
# check the samples look okay
curve(tempered_target_density_exp_4(x, beta = 1/4), -4, 4)
for (samples in input_samples) {
  lines(density(samples), col = 'blue')
}
# normal fusion
test2_standard <- fusion_diff_exp_4(N = 10000, time = 1, C = 4, samples_to_fuse = input_samples, betas = rep(1/4, 4), timed = T)
# time-adapting fusion
test2_TA <-fusion_TA_exp_4(N = 10000, time = 0.25, C = 4, samples_to_fuse = input_samples, betas = rep(1/4, 4), timed = T)
# Consensus Monte Carlo
samples_CMC <- CMC(N = 10000)
# plot kdes for fusion samples
curve(target_density_exp_4(x), -2.5, 2.5, lwd = 2)
plot(density(test_target_mc, adjust = 2), col = 'black', lwd = 2, main = '', xlab = 'x', ylim = c(0,1))
lines(density(test2_standard$samples, adjust = 2), col = 'red', lwd = 2)
lines(density(samples_CMC, adjust = 2), col = 'darkgrey', lwd = 2)
lines(density(test2_TA$samples), col = 'blue')



######################################## fusion example for fusing pi^(1/3) and pi^(1/4) times to get pi^(7/12)

# sampling from the tempered targets
nsamples <- 100000
tempered_samples_1 <- hmc_sample_tempered_exp_4(beta = 1/3, iterations = 2*nsamples, chains = 1, output = F)
tempered_samples_2 <- hmc_sample_tempered_exp_4(beta = 1/4, iterations = 2*nsamples, chains = 1, output = F)
input_samples <- list(tempered_samples_1, tempered_samples_2)
# check the samples look okay
curve(tempered_target_density_exp_4(x, beta = 1/3), -4, 4)
lines(density(tempered_samples_1), col = 'blue')
curve(tempered_target_density_exp_4(x, beta = 1/4), -4, 4)
lines(density(tempered_samples_2), col = 'blue')
# normal fusion
test3_standard <- fusion_diff_exp_4(N = 10000, time = 1, C = 2, samples_to_fuse = input_samples, betas = c(1/3, 1/4), timed = T)
# time-adapting fusion
test3_TA <- fusion_diff_exp_4(N = 10000, time = 1/3, C = 2, samples_to_fuse = input_samples, betas = c(1/3, 1/4), timed = T)
# plot kdes for fusion samples
curve(tempered_target_density_exp_4(x, beta = 7/12), -2.5, 2.5)
lines(density(test3_standard$samples), col = 'green')
lines(density(test3_TA$samples), col = 'blue')






########################################################################################################################

library(exp4Tempering)

nsamples <- 100000
target_nsamples <- 10000
test <- list()
samples_per_core <- rep(floor(target_nsamples/n_cores), n_cores)
if (sum(samples_per_core) != target_nsamples) {
  remainder <- target_nsamples %% n_cores
  samples_per_core[1:remainder] <- samples_per_core[1:remainder] + 1
}

for (denominator in 2:10) {
  print(denominator)
  input_samples <- hmc_base_sampler_exp_4(beta = 1/denominator, nsamples = nsamples, nchains = denominator)
  # check samples look okay
  curve(tempered_target_density_exp_4(x, beta = 1/denominator), -4, 4)
  for (samples in input_samples) {
    lines(density(samples), col = 'blue')
  }

  # creating parallel cluster
  n_cores <- parallel::detectCores()
  cl <- parallel::makeCluster(n_cores)
  # creating variable and functions list to pass into cluster using clusterExport
  varlist <- list("phi_function_exp_4", "bounds_phi_function_exp_4",
                  "simulate_langevin_diffusion_exp_4", "fusion_diff_exp_4", "input_samples",
                  "samples_per_core", "denominator")
  parallel::clusterExport(cl, envir = environment(), varlist = varlist)
  # exporting functions from layeredBB package to simulate layered Brownian bridges
  parallel::clusterExport(cl, varlist = ls("package:layeredBB"))

  pcm <- proc.time()
  fused <- parallel::parLapplyLB(cl, X = 1:length(samples_per_core), fun = function(i) {
    fusion_diff_exp_4(N = samples_per_core[i],
                      time = 1,
                      C = denominator,
                      samples_to_fuse = input_samples,
                      betas = 1/denominator)
  })
  final <- proc.time() - pcm

  test[[denominator-1]] <- list('samples' = unlist(fused), 'time' = final['elapsed'])

  parallel::stopCluster(cl)
}

curve(target_density_exp_4(x), -2.5, 2.5)
for (i in 1:9) {
  lines(density(test[[i]]$samples))
  print(paste(i, ':', test[[i]]$time))
}

plot(testing_times[1:length(h_fusion_TA_var)], sapply(h_fusion_TA_var, function(time_choice) sum(time_choice$time)), ylim = c(0,2000),
     col = 'blue', ylab = 'Computing Time', xlab = 'T', pch = 4, main = 'Hierarchical Fusion', xaxt="n")
axis(1, at = testing_times, las=1)
lines(testing_times[1:length(h_fusion_TA_var)], sapply(h_fusion_TA_var, function(time_choice) sum(time_choice$time)), col = 'blue')

plot(x = 1:7, y = log(sapply(1:7, function(i) test[[i]]$time)), col = 'blue', ylab = 'Computing Time in log(seconds)', xlab = 'Number of sub-posteriors (C)',
     pch = 4, xaxt='n', ylim = c(1.8, 5))
axis(1, at = 1:7, labels = 2:8)

plot(x = 1:7, y = sapply(1:7, function(i) test[[i]]$time), col = 'blue', ylab = 'Computing Time in seconds', xlab = 'Number of sub-posteriors (C)',
     pch = 4, xaxt='n', ylim = c(0, 100))
axis(1, at = 1:7, labels = 2:8)