scripts/seq_fusion_exp_4_example.R

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

library(exp4FusionRCPP)
library(fusionAnalysis)

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

set.seed(21)

# 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

# sampling from the tempered targets
# finding the right proposal distribution and constant M in rejection sampling
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, 0, 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_seq_standard <- parallel_sequential_fusion_exp_4_rcpp(N_schedule = rep(10000, 3),
                                                            mean = 0,
                                                            time_schedule = rep(1, 3),
                                                            start_beta = 1/4,
                                                            base_samples = input_samples1,
                                                            seed = 21)

test1_seq_TA <- parallel_sequential_fusion_TA_exp_4_rcpp(N_schedule = rep(10000, 3),
                                                         mean = 0,
                                                         global_T = 1,
                                                         start_beta = 1/4,
                                                         base_samples = input_samples1,
                                                         seed = 21)

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

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

# acceptance rate comparisons
acceptance_rate_plots(hier1 = test1_seq_TA, hier2 = test1_seq_standard, time = 1)




######################################## beta = 1/8

# 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.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, 0, 1/8), -5, 5, add = T, col = 'blue')

# using rejection sampling to obtain input samples
input_samples2 <- base_rejection_sampler_exp_4(beta = 1/8,
                                               nsamples = 100000,
                                               proposal_mean = 0,
                                               proposal_sd = 1.5,
                                               dominating_M = 1.3)
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_seq_standard <- parallel_sequential_fusion_exp_4_rcpp(N_schedule = rep(10000, 7),
                                                            mean = 0,
                                                            time_schedule = rep(1, 7),
                                                            start_beta = 1/8,
                                                            base_samples = input_samples2,
                                                            seed = 21)

test2_seq_TA <- parallel_sequential_fusion_TA_exp_4_rcpp(N_schedule = rep(10000, 7),
                                                         mean = 0,
                                                         global_T = 1,
                                                         start_beta = 1/8,
                                                         base_samples = input_samples2,
                                                         seed = 21)

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

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

# acceptance rate comparisons
acceptance_rate_plots(hier1 = test2_seq_TA, hier2 = test2_seq_standard, time = 1)

#####

par(mfrow=c(1,2))

#####
# seq

#rho
plot(1:length(test2_seq_standard$rho_acc), test2_seq_standard$rho_acc, ylim = c(0,1), col = 'black',
     ylab = expression(rho), xlab = 'Level', pch = 4, cex = 0.5)
lines(1:length(test2_seq_standard$rho_acc), test2_seq_standard$rho_acc, col = 'black')
points(1:length(test2_seq_TA$rho_acc), test2_seq_TA$rho_acc, col = 'orange', pch = 4, cex = 0.5)
lines(1:length(test2_seq_TA$rho_acc), test2_seq_TA$rho_acc, col = 'orange', lty = 2)

# Q
plot(1:length(test2_seq_standard$Q_acc), test2_seq_standard$Q_acc, ylim = c(0,1), col = 'black',
     ylab = 'Q', xlab = 'Level', pch = 4, cex = 0.5)
lines(1:length(test2_seq_standard$Q_acc), test2_seq_standard$Q_acc, col = 'black')
points(1:length(test2_seq_TA$Q_acc), test2_seq_TA$Q_acc, col = 'orange', pch = 4, cex = 0.5)
lines(1:length(test2_seq_TA$Q_acc), test2_seq_TA$Q_acc, col = 'orange', lty = 2)

#####

# rho
round(test2_seq_standard$rho_acc, 3)
round(test2_seq_TA$rho_acc, 3)

# Q
round(test2_seq_standard$Q_acc, 3)
round(test2_seq_TA$Q_acc, 3)