scripts/multivariate_Gaussian/bivG_verify/bivariate_Gaussian_bf_example_multi.R

In rchan26/hierarchicalFusion: Divide-and-Conquer Monte Carlo Fusion

library(DCFusion)
library(MASS)

seed <- 1983
set.seed(seed)
nsamples <- 10000
fusion_time <- 1
mean <- c(1, 2)
sd <- c(sqrt(0.5), sqrt(2))
C <- 4
cov_mat <- C*matrix(c(sd[1]^2, 0.9, 0.9, sd[2]^2), nrow = 2, ncol = 2, byrow = T)
diffusion_estimator <- 'NB'

# sampling from the sub-posteriors (target at inverse temperature 1/4)
input_samples <- lapply(1:C, function(i) mvrnormArma(1000000, mu = mean, Sigma = cov_mat))

# sample from true target density
true_samples <- mvrnormArma(1000000, mu = mean, Sigma = cov_mat/C)
true_kde <- MASS::kde2d(true_samples[,1], true_samples[,2], n = 50)
image(true_kde)
contour(true_kde, add = T)

##### Monte Carlo Fusion #####
input_particles_MCF <- initialise_particle_sets(samples_to_fuse = input_samples,
                                                multivariate = TRUE)

print('### performing standard fusion')
MCF_standard <- parallel_fusion_SMC_multiGaussian(particles_to_fuse = input_particles_MCF,
                                                  N = nsamples,
                                                  m = C,
                                                  time = fusion_time,
                                                  dim = 2,
                                                  mu = mean,
                                                  Sigma = cov_mat,
                                                  betas = rep(1, C),
                                                  precondition_matrices = rep(list(diag(1,2)), C),
                                                  ESS_threshold = 0.5,
                                                  diffusion_estimator = diffusion_estimator,
                                                  seed = seed)
print('ESS:'); print(MCF_standard$ESS)
print('CESS:'); print(MCF_standard$CESS)

print('### performing fusion with a preconditioning matrix')
MCF_generalised <- parallel_fusion_SMC_multiGaussian(particles_to_fuse = input_particles_MCF,
                                                     N = nsamples,
                                                     m = C,
                                                     time = fusion_time,
                                                     dim = 2,
                                                     mu = mean,
                                                     Sigma = cov_mat,
                                                     betas = rep(1, C),
                                                     precondition_matrices = lapply(input_samples, cov),
                                                     ESS_threshold = 0.5,
                                                     diffusion_estimator = diffusion_estimator,
                                                     seed = seed)
print('ESS:'); print(MCF_generalised$ESS)
print('CESS:'); print(MCF_generalised$CESS)

# plots
MCF_standard$particles <- resample_particle_y_samples(particle_set = MCF_standard$particles,
                                                      multivariate = TRUE,
                                                      resampling_method = 'resid',
                                                      seed = seed)
MCF_generalised$particles <- resample_particle_y_samples(particle_set = MCF_generalised$particles,
                                                         multivariate = TRUE,
                                                         resampling_method = 'resid',
                                                         seed = seed)
compare_samples_bivariate(list(true_samples,
                               MCF_standard$particles$y_samples,
                               MCF_generalised$particles$y_samples),
                          c('black', 'red', 'blue'),
                          c(-4, 4))

##### Bayesian Fusion (with n=1, so equal to Monte Carlo Fusion) #####
time_mesh_BF_n1 <- seq(0, fusion_time, 1)
input_particles_BF_n1 <- initialise_particle_sets(samples_to_fuse = input_samples,
                                                  multivariate = TRUE,
                                                  number_of_steps = length(time_mesh_BF_n1))

print('### performing standard Bayesian Fusion (with n=1)')
BF_standard_n1 <- parallel_GBF_multiGaussian(particles_to_fuse = input_particles_BF_n1,
                                             N = nsamples,
                                             m = C,
                                             time_mesh = time_mesh_BF_n1,
                                             dim = 2,
                                             mean_vecs = rep(list(mean), C),
                                             Sigmas = rep(list(cov_mat), C),
                                             precondition_matrices = rep(list(diag(1,2)), C),
                                             ESS_threshold = 0.5,
                                             diffusion_estimator = diffusion_estimator,
                                             seed = seed)
print('ESS:'); print(BF_standard_n1$ESS)
print('CESS:'); print(BF_standard_n1$CESS)

print('### performing Bayesian Fusion (with n=1) with a preconditioning matrix')
BF_generalised_n1 <- parallel_GBF_multiGaussian(particles_to_fuse = input_particles_BF_n1,
                                                N = nsamples,
                                                m = C,
                                                time_mesh = time_mesh_BF_n1,
                                                dim = 2,
                                                mean_vecs = rep(list(mean), C),
                                                Sigmas = rep(list(cov_mat), C),
                                                precondition_matrices = lapply(input_samples, cov),
                                                ESS_threshold = 0.5,
                                                diffusion_estimator = diffusion_estimator,
                                                seed = seed)
print('ESS:'); print(BF_generalised_n1$ESS)
print('CESS:'); print(BF_generalised_n1$CESS)

# plots
BF_standard_n1$particles <- resample_particle_y_samples(particle_set = BF_standard_n1$particles,
                                                        multivariate = TRUE,
                                                        resampling_method = 'resid',
                                                        seed = seed)
BF_generalised_n1$particles <- resample_particle_y_samples(particle_set = BF_generalised_n1$particles,
                                                           multivariate = TRUE,
                                                           resampling_method = 'resid',
                                                           seed = seed)
compare_samples_bivariate(list(true_samples,
                               BF_standard_n1$particles$y_samples,
                               BF_generalised_n1$particles$y_samples),
                          c('black', 'red', 'blue'),
                          c(-4, 4))

##### Bayesian Fusion (with n=20) #####
time_mesh_BF_n20 <- seq(0, fusion_time, fusion_time/20)
input_particles_BF_n20 <- initialise_particle_sets(samples_to_fuse = input_samples,
                                                   multivariate = TRUE,
                                                   number_of_steps = length(time_mesh_BF_n20))

print('### performing standard Bayesian Fusion (with n=20)')
BF_standard_n20 <- parallel_GBF_multiGaussian(particles_to_fuse = input_particles_BF_n20,
                                              N = nsamples,
                                              m = C,
                                              time_mesh = time_mesh_BF_n20,
                                              dim = 2,
                                              mean_vecs = rep(list(mean), C),
                                              Sigmas = rep(list(cov_mat), C),
                                              precondition_matrices = rep(list(diag(1,2)), C),
                                              ESS_threshold = 0.5,
                                              diffusion_estimator = diffusion_estimator,
                                              seed = seed)
print('ESS:'); print(BF_standard_n20$ESS)
print('CESS:'); print(BF_standard_n20$CESS)

print('### performing Bayesian Fusion (with n=20) with a preconditioning matrix')
BF_generalised_n20 <- parallel_GBF_multiGaussian(particles_to_fuse = input_particles_BF_n20,
                                                 N = nsamples,
                                                 m = C,
                                                 time_mesh = time_mesh_BF_n20,
                                                 dim = 2,
                                                 mean_vecs = rep(list(mean), C),
                                                 Sigmas = rep(list(cov_mat), C),
                                                 precondition_matrices = lapply(input_samples, cov),
                                                 ESS_threshold = 0.5,
                                                 diffusion_estimator = diffusion_estimator,
                                                 seed = seed)
print('ESS:'); print(BF_generalised_n20$ESS)
print('CESS:'); print(BF_generalised_n20$CESS)

# plots
# proposals
compare_samples_bivariate(list(true_samples,
                               BF_standard_n20$proposed_samples,
                               BF_generalised_n20$proposed_samples),
                          c('black', 'red', 'blue'),
                          c(-4, 4))
# resampled
BF_standard_n20$particles <- resample_particle_y_samples(particle_set = BF_standard_n20$particles,
                                                         multivariate = TRUE,
                                                         resampling_method = 'resid',
                                                         seed = seed)
BF_generalised_n20$particles <- resample_particle_y_samples(particle_set = BF_generalised_n20$particles,
                                                            multivariate = TRUE,
                                                            resampling_method = 'resid',
                                                            seed = seed)
compare_samples_bivariate(list(true_samples,
                               BF_standard_n20$particles$y_samples,
                               BF_generalised_n20$particles$y_samples),
                          c('black', 'red', 'blue'),
                          c(-4, 4))

save.image('bivariate_Gaussian_bf_example_multi.RData')