scripts/study_time_variance_exp_4.R

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

library(exp4Tempering)
library(fusionAnalysis)

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

variance <- var(unlist(base))

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

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

h_fusion_standard_var <- list()
for (i in 1:length(testing_times)) {
  print(paste('standard - time_choice:', testing_times[i]))
  h_fusion_standard_var[[i]] <- parallel_h_fusion_exp_4(N_schedule = rep(10000, 3),
                                                        time_schedule = rep(testing_times[i], 3),
                                                        start_beta = 1/8,
                                                        C_schedule = rep(2, 3),
                                                        L = 4,
                                                        base_samples = base,
                                                        study = T)
}

#####

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

# ##### plots
for (i in 1:length(testing_times)) {
  acceptance_rate_plots(hier1 = h_fusion_TA_var[[i]], hier2 = h_fusion_standard_var[[i]], time = testing_times[i])
}
par(mfrow=c(1,1))
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')
points(testing_times[1:length(h_fusion_standard_var)], sapply(h_fusion_standard_var, function(time_choice) sum(time_choice$time)), col = 'green', pch = 4)
lines(testing_times[1:length(h_fusion_standard_var)], sapply(h_fusion_standard_var, 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_var[[1]]$samples[[1]]))
for (i in 1:length(testing_times)) {
  lines(density(h_fusion_standard_var[[i]]$samples[[1]]), col = 'green')
  lines(density(h_fusion_TA_var[[i]]$samples[[1]]), col = 'blue')
}

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

seq_fusion_standard_var <- list()
for (i in 1:length(testing_times)) {
  print(paste('standard - time_choice:', testing_times[i]))
  seq_fusion_standard_var[[i]] <- parallel_sequential_fusion_exp_4(N_schedule = rep(10000, 7),
                                                                   time_schedule = rep(testing_times[i], 7),
                                                                   start_beta = 1/8,
                                                                   base_samples = base,
                                                                   study = T)
}

#####

seq_fusion_TA_var <- list()
for(i in 1:length(testing_times)) {
  print(paste('time adapting - time_choice:', testing_times[i]))
  seq_fusion_TA_var[[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 = seq_fusion_TA_var[[i]], hier2 = seq_fusion_standard_var[[i]], time = testing_times[i])
}
par(mfrow=c(1,1))
plot(testing_times[1:length(seq_fusion_TA_var)], sapply(seq_fusion_TA_var, function(time_choice) sum(time_choice$time)), ylim = c(0,500),
     col = 'blue', ylab = 'Computing Time', xlab = 'T', pch = 4, main = 'Sequential Fusion', xaxt="n")
axis(1, at = testing_times, las=1)
lines(testing_times[1:length(seq_fusion_TA_var)], sapply(seq_fusion_TA_var, function(time_choice) sum(time_choice$time)), col = 'blue')
points(testing_times[1:length(seq_fusion_standard_var)], sapply(seq_fusion_standard_var, function(time_choice) sum(time_choice$time)), col = 'green', pch = 4)
lines(testing_times[1:length(seq_fusion_standard_var)], sapply(seq_fusion_standard_var, 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_var[[i]]$samples[[1]]), col = 'green')
  lines(density(seq_fusion_TA_var[[i]]$samples[[1]]), col = 'blue')
}



dev.new();
plot(x=1:7, y=seq_fusion_TA[[2]]$rho_acc, ylab = 'Rho', xlab = 'Level', pch = 4, cex = 0.5, ylim = c(0,1))
lines(x=1:7, y=seq_fusion_TA[[2]]$rho_acc)
plot(x=1:7, y=seq_fusion_TA[[2]]$Q_acc, ylab = 'Q', xlab = 'Level', pch = 4, cex = 0.5, ylim = c(0,1))
lines(x=1:7, y=seq_fusion_TA[[2]]$Q_acc)
plot(x=1:7, y=seq_fusion_TA[[2]]$rhoQ_acc, ylab = 'rhoQ', xlab = 'Level', pch = 4, cex = 0.5, ylim = c(0,1))
lines(x=1:7, y=seq_fusion_TA[[2]]$rhoQ_acc)


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



######################################## line ########################################

line_fusion_standard_var <- list()
for (i in 1:length(testing_times)) {
  print(paste('standard - time_choice:', testing_times[i]))
  line_fusion_standard_var[[i]] <- parallel_line_fusion_exp_4(N_schedule = rep(10000, 3),
                                                              time_schedule = rep(testing_times[i], 3),
                                                              start_beta = 1/8,
                                                              C_schedule = rep(2, 3),
                                                              L = 4,
                                                              base_samples = base[[1]],
                                                              study = T)
}

#####

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

##### plots comparing standard hierarchical and line
for (i in 1:length(testing_times)) {
  acceptance_rate_plots(hier1 = h_fusion_standard_var[[i]], hier2 = line_fusion_standard_var[[i]], colours = c('blue', 'green'), time = testing_times[i])
}
par(mfrow=c(1,1))
plot(testing_times[1:length(h_fusion_standard_var)], sapply(h_fusion_standard_var, function(time_choice) sum(time_choice$time)), ylim = c(0,500),
     col = 'blue', ylab = 'Computing Time', xlab = 'T', pch = 4, main = 'Hierarchical (Blue) vs. Line (Green) Standard Fusion', xaxt="n")
axis(1, at = testing_times, las=1)
lines(testing_times[1:length(h_fusion_standard_var)], sapply(h_fusion_standard_var, function(time_choice) sum(time_choice$time)), col = 'blue')
points(testing_times[1:length(line_fusion_standard_var)], sapply(line_fusion_standard_var, function(time_choice) sum(time_choice$time)), col = 'green', pch = 4)
lines(testing_times[1:length(line_fusion_standard_var)], sapply(line_fusion_standard_var, function(time_choice) sum(time_choice$time)), col = 'green')

##### plots comparing TA hierarchical and line
for (i in 1:length(testing_times)) {
  acceptance_rate_plots(hier1 = h_fusion_TA_var[[i]], hier2 = line_fusion_TA_var[[i]], colours = c('blue', 'green'), time = testing_times[i])
}
par(mfrow=c(1,1))
plot(testing_times[1:length(h_fusion_TA_var)], sapply(h_fusion_TA_var, function(time_choice) sum(time_choice$time)), ylim = c(0,500),
     col = 'blue', ylab = 'Computing Time', xlab = 'T', pch = 4, main = 'Hierarchical (Blue) vs. Line (Green) TA 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')
points(testing_times[1:length(line_fusion_TA_var)], sapply(line_fusion_TA_var, function(time_choice) sum(time_choice$time)), col = 'green', pch = 4)
lines(testing_times[1:length(line_fusion_TA_var)], sapply(line_fusion_TA_var, 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_var[[i]]$samples[[1]]), col = 'green')
  lines(density(seq_fusion_TA_var[[i]]$samples[[1]]), col = 'blue')
}

# get total variation
for (i in 1:6) {
  h_fusion_standard_var[[i]]$tv <- total_variation_true_exp_4(h_fusion_standard_var[[i]]$samples[[1]])
  line_fusion_standard_var[[i]]$tv <- total_variation_true_exp_4(line_fusion_standard_var[[i]]$samples[[1]])
  seq_fusion_standard_var[[i]]$tv <- total_variation_true_exp_4(seq_fusion_standard_var[[i]]$samples[[1]])
}

for (i in 1:7) {
  h_fusion_TA_var[[i]]$tv <- total_variation_true_exp_4(h_fusion_TA_var[[i]]$samples[[1]])
  line_fusion_TA_var[[i]]$tv <- total_variation_true_exp_4(line_fusion_TA_var[[i]]$samples[[1]])
  seq_fusion_TA_var[[i]]$tv <- total_variation_true_exp_4(seq_fusion_TA_var[[i]]$samples[[1]])
}