scripts/seq_fusion_examples.R

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

library(mixGaussTempering)
library(fusionAnalysis)

######################################## example usage ########################################

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

# normal seq fusion
example_seq_standard <- parallel_sequential_fusion_mixG(N_schedule = rep(10000, 31), time_schedule = rep(0.5, 31),
                                                        start_beta = b_test, base_samples = base,
                                                        weights = w_test, means = m_test, sds = s_test, study = T)

# time-adapting fusion
example_seq_TA <- parallel_sequential_fusion_TA_mixG(N_schedule = rep(10000, 31), global_T = 0.5, start_beta = b_test,
                                                     base_samples = base, weights = w_test, means = m_test, sds = s_test,
                                                     study = T)

# plot the kde for the target (top level)
curve(dnorm_mix(x, w_test, m_test, s_test), -20, 30, col = 'red', ylab = 'pdf')
lines(density(example_seq_standard$samples[[1]], adjust = 0.1), col = 'blue')
lines(density(example_seq_TA$samples[[1]], adjust = 0.1), col = 'green')


# acceptance rate plots
acceptance_rate_plots(hier1_blue = example_seq_standard, hier2_green = example_seq_TA, time = 1)