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_hier_standard <- parallel_h_fusion_exp_4_rcpp(N_schedule = rep(10000, 2),
mean = 0,
time_schedule = rep(1, 2),
start_beta = 1/4,
C_schedule = rep(2, 2),
L = 3,
base_samples = input_samples1,
seed = 21)
test1_hier_TA <- parallel_h_fusion_TA_exp_4_rcpp(N_schedule = rep(10000, 2),
mean = 0,
global_T = 1,
start_beta = 1/4,
C_schedule = rep(2, 2),
L = 3,
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_hier_standard$samples[[1]]), col = 'green')
lines(density(test1_hier_TA$samples[[1]]), col = 'blue')
# plot all levels
plot_levels_hier_exp_4(test1_hier_standard, from = -3, to = 3, plot_rows = 3, plot_columns = 1)
plot_levels_hier_exp_4(test1_hier_TA, from = -3, to = 3, plot_rows = 3, plot_columns = 1)
# acceptance rate comparisons
acceptance_rate_plots(hier1 = test1_hier_TA, hier2 = test1_hier_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_hier_standard <- parallel_h_fusion_exp_4_rcpp(N_schedule = rep(10000, 3),
mean = 0,
time_schedule = rep(1, 3),
start_beta = 1/8,
C_schedule = rep(2, 3),
L = 4,
base_samples = input_samples2,
seed = 21)
test2_hier_TA <- parallel_h_fusion_TA_exp_4_rcpp(N_schedule = rep(10000, 3),
mean = 0,
global_T = 1,
start_beta = 1/8,
C_schedule = rep(2, 3),
L = 4,
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_hier_standard$samples[[1]]), col = 'green')
lines(density(test2_hier_TA$samples[[1]]), col = 'blue')
# plot all levels
plot_levels_hier_exp_4(test2_hier_standard, from = -3, to = 3, plot_rows = 2, plot_columns = 2)
plot_levels_hier_exp_4(test2_hier_TA, from = -3, to = 3, plot_rows = 2, plot_columns = 2)
# acceptance rate comparisons
acceptance_rate_plots(hier1 = test2_hier_TA, hier2 = test2_hier_standard, time = 1)
#####
par(mfrow=c(1,2))
#####
# hier
#rho
plot(1:length(test2_hier_standard$rho_acc), test2_hier_standard$rho_acc, ylim = c(0,1), col = 'black',
ylab = expression(rho), xlab = 'Level', pch = 4, cex = 0.5)
lines(1:length(test2_hier_standard$rho_acc), test2_hier_standard$rho_acc, col = 'black')
points(1:length(test2_hier_TA$rho_acc), test2_hier_TA$rho_acc, col = 'orange', pch = 4, cex = 0.5)
lines(1:length(test2_hier_TA$rho_acc), test2_hier_TA$rho_acc, col = 'orange', lty = 2)
# Q
plot(1:length(test2_hier_standard$Q_acc), test2_hier_standard$Q_acc, ylim = c(0,1), col = 'black',
ylab = 'Q', xlab = 'Level', pch = 4, cex = 0.5)
lines(1:length(test2_hier_standard$Q_acc), test2_hier_standard$Q_acc, col = 'black')
points(1:length(test2_hier_TA$Q_acc), test2_hier_TA$Q_acc, col = 'orange', pch = 4, cex = 0.5)
lines(1:length(test2_hier_TA$Q_acc), test2_hier_TA$Q_acc, col = 'orange', lty = 2)
#####
# rho
round(test2_hier_standard$rho_acc, 3)
round(test2_hier_TA$rho_acc, 3)
# Q
round(test2_hier_standard$Q_acc, 3)
round(test2_hier_TA$Q_acc, 3)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.