scripts/study_time_exp_4.R

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

##### studying the best time T
source('h_fusion_exp_4.R')
source('seq_fusion_exp_4.R')

##### set seed
set.seed(26071996)

##### sampling from the tempered density to obtain samples for the base level using rejection sampling
base <- base_rejection_sampler(beta = 1/8, nsamples = 100000, proposal_mean = 0, 
                               proposal_sd = 1.5,dominating_M = 1.3)
curve(tempered_target_density(x, beta = 1/8), -4, 4)
for (i in 1:8) lines(density(base[[i]], adjust = 1.5), col = 'blue')

# times we would like to try out
testing_times <- c(0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.5)

######################################## hierarchical ########################################

h_fusion_standard <- lapply(testing_times, function(time_choice) {
  print(paste('standard - time_choice:', time_choice))
  parallel_h_fusion_exp_4(N_schedule = rep(10000, 3),
                          time_schedule = rep(time_choice, 3),
                          start_beta = 1/8,
                          m_schedule = rep(2, 3),
                          L = 4,
                          base_samples = base,
                          study = T)})

#####

h_fusion_TA <- list()
for (i in 1:length(testing_times)) {
  print(paste('time adapting - time_choice:', testing_times[i]))
  h_fusion_TA[[i]] <- parallel_h_fusion_TA_exp_4(N_schedule = rep(10000, 3),
                                                 global_T = testing_times[i],
                                                 start_beta = 1/8,
                                                 m_schedule = rep(2, 3),
                                                 L = 4,
                                                 base_samples = base,
                                                 study = T)
}

# ##### plots
for (i in 1:length(testing_times)) {
  acceptance_rate_plots(hier1_blue = h_fusion_TA[[i]], hier2_green = h_fusion_standard[[i]], time = testing_times[i])
}
par(mfrow=c(1,1))
plot(testing_times[1:length(h_fusion_TA)], sapply(h_fusion_TA, function(time_choice) sum(time_choice$time)), ylim = c(0,200),
     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)], sapply(h_fusion_TA, function(time_choice) sum(time_choice$time)), col = 'blue')
points(testing_times[1:length(h_fusion_standard)], sapply(h_fusion_standard, function(time_choice) sum(time_choice$time)), col = 'green', pch = 4)
lines(testing_times[1:length(h_fusion_standard)], sapply(h_fusion_standard, function(time_choice) sum(time_choice$time)), col = 'green')

##### checking that they all target the correct posterior
curve(target_density(x), -4, 4, ylim = c(0, 0.5))
lines(density(h_fusion_TA[[1]]$samples[[1]]))
for (i in 1:length(testing_times)) {
  lines(density(h_fusion_standard[[i]]$samples[[1]]), col = 'green')
  lines(density(h_fusion_TA[[i]]$samples[[1]]), col = 'blue')
}

######################################## sequential ########################################

seq_fusion_standard <- lapply(testing_times, function(time_choice) {
  print(paste('standard - time_choice:', time_choice))
  parallel_sequential_fusion_exp_4(N_schedule = rep(10000, 7), 
                                   time_schedule = rep(time_choice, 7), 
                                   start_beta = 1/8, 
                                   base_samples = base, 
                                   study = T)})

#####

seq_fusion_TA <- list()
for(i in 1:length(testing_times)) {
  print(paste('time adapting - time_choice:', testing_times[i]))
  seq_fusion_TA[[i]] <- parallel_sequential_fusion_TA_exp_4(N_schedule = rep(10000, 7), 
                                                            global_T = testing_times[i], 
                                                            start_beta = 1/8, 
                                                            base_samples = base, 
                                                            study = T)
}

##### plots
for (i in 1:length(testing_times)) {
  acceptance_rate_plots(hier1_blue = seq_fusion_TA[[i]], hier2_green = seq_fusion_standard[[i]], time = testing_times[i])
}
par(mfrow=c(1,1))
plot(testing_times[1:length(seq_fusion_TA)], sapply(seq_fusion_TA, function(time_choice) sum(time_choice$time)), ylim = c(0,500),
     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(seq_fusion_TA)], sapply(seq_fusion_TA, function(time_choice) sum(time_choice$time)), col = 'blue')
points(testing_times[1:length(seq_fusion_standard)], sapply(seq_fusion_standard, function(time_choice) sum(time_choice$time)), col = 'green', pch = 4)
lines(testing_times[1:length(seq_fusion_standard)], sapply(seq_fusion_standard, function(time_choice) sum(time_choice$time)), col = 'green')

##### checking that they all target the correct posterior
curve(target_density(x), -2.5, 2.5, ylim = c(0, 0.5))
for (i in 1:length(testing_times)) {
  lines(density(seq_fusion_standard[[i]]$samples[[1]]), col = 'green')
  lines(density(seq_fusion_TA[[i]]$samples[[1]]), col = 'blue')
}





##### save
save.image('study_time_exp_4.RData')