scripts/h_fusion_examples.R

In rchan26/mixGaussTempering: What the Package Does (One Line, Title Case)

library(mixGaussTempering)
library(fusionAnalysis)

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

# ##### setting variables for mixture Normal distribution
w_test <- c(0.2, 0.35, 0.45)
m_test <- c(-8, -3, 18)
s_test <- c(1, 1, 1)
b_test <- 1/32

curve(dnorm_mix_tempered(x, w_test, m_test, s_test, b_test), -30, 40)

###### sampling from the tempered density to obtain samples for the base level
nsamples <- 500000
base <- hmc_base_sampler_mixG(w_test, m_test, s_test, b_test, nsamples, 32)

# check the base samples look okay
for (samples in base) {
  lines(density(samples), col = 'blue')
}

##### performing hierarchical fusion for when the number of samples to fuse at each stage is CONSTANT (m=2)
example_h_standard <- parallel_h_fusion_TA_mixG(N_schedule = rep(10000, 5), time_schedule = rep(0.5, 5),
                                                start_beta = b_test, C_schedule = rep(2, 5), L = 6, base_samples = base,
                                                weights = w_test, means = m_test, sds = s_test, study = T)

# plot results
plot_levels_hier_mixG(example_h_standard, weights = w_test, means = m_test, sds = s_test,
                      from = -30, to = 40, plot_rows = 3, plot_columns = 2)



##### performing hierarchical fusion for when the number of samples to fuse at each stage is *NOT* constant
example_h_non_constant <- parallel_h_fusion_mixG(N_schedule = rep(10000, 3), time_schedule = rep(0.5, 3),
                                                 start_beta = b_test, C_schedule = c(4,4,2), L = 4, base_samples = base,
                                                 weights = w_test, means = m_test, sds = s_test, study = T)

# plot results
plot_levels_hier_mixG(example_h_non_constant, weights = w_test, means = m_test, sds = s_test,
                      from = -30, to = 40, plot_rows = 2, plot_columns = 2)




##### time-adapting
example_h_TA <- parallel_h_fusion_TA_mixG(N_schedule = rep(10000, 5), global_T = 0.5, start_beta = b_test,
                                          C_schedule = rep(2, 5), L = 6, base_samples = base,
                                          weights = w_test, means = m_test, sds = s_test, study = T)

# plot results
plot_levels_hier_mixG(example_h_TA, weights = w_test, means = m_test, sds = s_test,
                      from = -30, to = 40, plot_rows = 3, plot_columns = 2)



##### time-adapting for when the number of samples to fuse at each stage is *NOT* constant
example_h_TA_non_constant <- parallel_h_fusion_TA_mixG(N_schedule = rep(10000, 3), global_T = 0.5, start_beta = b_test,
                                                       C_schedule = c(4,4,2), L = 4, base_samples = base,
                                                       weights = w_test, means = m_test, sds = s_test, study = T)

# plot results
plot_levels_hier_mixG(example_h_TA_non_constant, weights = w_test, means = m_test, sds = s_test,
                      from = -30, to = 40, plot_rows = 2, plot_columns = 2)


# acceptance rate plots
acceptance_rate_plots(hier1_blue = example_h_standard, hier2_green = example_h_TA, time = time_choice)
acceptance_rate_plots(hier1_blue = example_h_non_constant, hier2_green = example_h_TA_non_constant, time = time_choice)