scripts/seq_fusion_examples.R

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

library(exp4Tempering)
library(fusionAnalysis)

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



######################################## beta = 1/4
dominating_dnorm <- function(x) 1.4*dnorm(x, mean = 0, sd = 1.5)
curve(dominating_dnorm(x), -5, 5, col = 'red')
curve(dnorm(x, mean = 0, sd = 1.5), -5, 5, add = T, col = 'green')
curve(tempered_target_density_exp_4(x, 1/4), -5, 5, add = T, col = 'blue')

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

test1_standard <- parallel_sequential_fusion_exp_4(N_schedule = rep(10000, 3), time_schedule = rep(1, 3),
                                                   start_beta = 1/4, base_samples = input_samples1, study = T)
test1_TA <- parallel_sequential_fusion_TA_exp_4(N_schedule = rep(10000, 3), global_T = 1,
                                                start_beta = 1/4, base_samples = input_samples1, study = T)

curve(target_density_exp_4(x), -2.5, 2.5, ylim = c(0,1))
lines(density(test_target_mc), col = 'black')
lines(density(test1_standard$samples[[1]]), col = 'green')
lines(density(test1_TA$samples[[1]]), col = 'blue')

# plot all levels
plot_levels_seq_exp_4(test1_standard, from = -3, to = 3, plot_rows = 2, plot_columns = 2)
plot_levels_seq_exp_4(test1_TA, from = -3, to = 3, plot_rows = 2, plot_columns = 2)

# acceptance rate comparisons
acceptance_rate_plots(hier1 = test1_TA, hier2 = test1_standard, time = 1)




##### beta = 1/8

# using Stan to sample from base level
input_samples2 <- hmc_base_sampler_exp_4(beta = 1/8, nsamples = 100000, nchains = 8)
curve(tempered_target_density_exp_4(x, beta = 1/8), -4, 4)
# check the samples look okay
for (samples in input_samples2) {
  lines(density(samples), col = 'blue')
}

test2_standard <- parallel_sequential_fusion_exp_4(N_schedule = rep(10000, 7), time_schedule = rep(1, 7),
                                                   start_beta = 1/8, base_samples = input_samples2, study = T)
test2_TA <- parallel_sequential_fusion_TA_exp_4(N_schedule = rep(10000, 7), global_T = 1,
                                                start_beta = 1/8, base_samples = input_samples2, study = T)

curve(target_density_exp_4(x), -2.5, 2.5)
lines(density(test_target_mc), col = 'black')
lines(density(test2_standard$samples[[1]]), col = 'green')
lines(density(test2_TA$samples[[1]]), col = 'blue')

# plot all levels
plot_levels_seq_exp_4(test2_standard, from = -3, to = 3, plot_rows = 2, plot_columns = 2)
plot_levels_seq_exp_4(test2_TA, from = -3, to = 3, plot_rows = 2, plot_columns = 2)

# acceptance rate comparisons
acceptance_rate_plots(hier1 = test2_TA, hier2 = test2_standard, time = 1)