#' rho_j Importance Sampling Step
#'
#' rho_j Importance Sampling weighting for bivariate Gaussian distributions
#'
#' @param particle_set particles set prior to Q importance sampling step
#' @param m number of sub-posteriors to combine
#' @param time_mesh time mesh used in Bayesian Fusion
#' @param mean_vecs list of length m, where mean_vecs[[c]] is a vector of
#' length 2 for the mean of sub-posterior c
#' @param sd_vecs list of length m, where sd_vecs[[c]] is a vector of length 2
#' for the standard deviation of sub-posterior c
#' @param corrs vector of length m, where corrs[c] give the correlation value
#' between component 1 and component 2 for sub-posterior c
#' @param betas vector of length m, where betas[c] is the inverse temperature
#' value for c-th posterior
#' @param precondition_matrices list of length m, where precondition_matrices[[c]]
#' is the precondition matrix for sub-posterior c
#' @param inv_precondition_matrices list of length m, where inv_precondition_matrices[[c]]
#' is the inverse precondition matrix for sub-posterior c
#' @param Lambda inverse of the sum of the inverse precondition matrices (which
#' can be computed using inverse_sum_matrices(inv_precondition_matrices))
#' @param resampling_method method to be used in resampling, default is multinomial
#' resampling ('multi'). Other choices are stratified
#' resampling ('strat'), systematic resampling ('system'),
#' residual resampling ('resid')
#' @param ESS_threshold number between 0 and 1 defining the proportion of the
#' number of samples that ESS needs to be lower than for
#' resampling (i.e. resampling is carried out only when
#' ESS < N*ESS_threshold)
#' @param sub_posterior_means matrix with m rows and 2 columns, where sub_posterior_means[c,]
#' is the sub-posterior mean of sub-posterior c
#' @param adaptive_mesh logical value to indicate if an adaptive mesh is used
#' (default is FALSE)
#' @param adaptive_mesh_parameters list of parameters used for adaptive mesh
#' @param diffusion_estimator choice of unbiased estimator for the Exact Algorithm
#' between "Poisson" (default) for Poisson estimator
#' and "NB" for Negative Binomial estimator
#' @param beta_NB beta parameter for Negative Binomial estimator (default 10)
#' @param gamma_NB_n_points number of points used in the trapezoidal estimation
#' of the integral found in the mean of the negative
#' binomial estimator (default is 2)
#' @param seed seed number - default is NULL, meaning there is no seed
#' @param n_cores number of cores to use
#'
#' @return A list with components:
#' \describe{
#' \item{particle_set}{updated particle set after the iterative rho_j steps}
#' \item{proposed_samples}{proposal samples for the last time step}
#' \item{time}{elapsed time of each step of the algorithm}
#' \item{ESS}{effective sample size of the particles after each step}
#' \item{CESS}{conditional effective sample size of the particles after each step}
#' \item{resampled}{boolean value to indicate if particles were resampled
#' after each time step}
#' \item{E_nu_j}{approximation of the average variation of the trajectories
#' at each time step}
#' \item{chosen}{which term was chosen if using an adaptive mesh at each time step}
#' \item{mesh_terms}{the evaluated terms in deciding the mesh at each time step}
#' }
#'
#' @export
<- (particle_set,
m,
time_mesh,
mean_vecs,
sd_vecs,
corrs,
betas,
precondition_matrices,
inv_precondition_matrices,
Lambda,
resampling_method = 'multi',
ESS_threshold = 0.5,
sub_posterior_means = ,
adaptive_mesh = ,
adaptive_mesh_parameters = ,
diffusion_estimator = 'Poisson',
beta_NB = 10,
gamma_NB_n_points = 2,
seed = ,
n_cores = ::()) {
if (!("particle" %in% (particle_set))) {
("rho_j_biGaussian: particle_set must be a \"particle\" object")
} if (!(time_mesh)) {
("rho_j_biGaussian: time_mesh must be an ordered vector of length >= 2")
} if ((time_mesh) < 2) {
("rho_j_biGaussian: time_mesh must be an ordered vector of length >= 2")
} if (!(time_mesh, (time_mesh))) {
("rho_j_biGaussian: time_mesh must be an ordered vector of length >= 2")
} if (((1:m, () (!(mean_vecs[c]]) | ((mean_vecs[c]])!=2))))) {
("rho_j_biGaussian: mean_vecs[[c]] must be a vector of length 2 for each c")
} if (((1:m, () (!(sd_vecs[c]]) | ((sd_vecs[c]])!=2))))) {
("rho_j_biGaussian: sd_vecs[[c]] must be a vector of length 2 for each c")
} if (!(corrs) | ((corrs)!=m)) {
("rho_j_biGaussian: corrs must be a vector of length m")
} if (!(betas) | ((betas)!=m)) {
("rho_j_biGaussian: betas must be a vector of length m")
} if (!(precondition_matrices) | ((precondition_matrices)!=m)) {
("rho_j_biGaussian: precondition_matrices must be a list of length m")
} if (!(inv_precondition_matrices) | ((inv_precondition_matrices)!=m)) {
("rho_j_biGaussian: inv_precondition_matrices must be a list of length m")
} if ((ESS_threshold < 0) | (ESS_threshold > 1)) {
("rho_j_biGaussian: ESS_threshold must be between 0 and 1")
}
if (adaptive_mesh) {
if (!(sub_posterior_means)) {
("rho_j_biGaussian: if adaptive_mesh==TRUE, sub_posterior_means must be a (m x 2) matrix")
} if (((sub_posterior_means)!=(m,2))) {
("rho_j_biGaussian: if adaptive_mesh==TRUE, sub_posterior_means must be a (m x 2) matrix")
}
}
transform_matrices <- (1:m, () {
('to_Z' = expm::sqrtm(inv_precondition_matrices[c]]),
'to_X' = expm::sqrtm(precondition_matrices[c]]))
})
N <- particle_set$N
# ---------- creating parallel cluster
cl <- ::(n_cores, setup_strategy = "sequential")
::(cl, envir = (), varlist = ())
::(cl, varlist = ("package:DCFusion"))
::(cl, varlist = ("package:layeredBB"))
if (!(seed)) {
::(cl, iseed = seed)
}
max_samples_per_core <- (N/n_cores)
split_indices <- (1:N, ((1:N)/max_samples_per_core))
elapsed_time <- (, (time_mesh)-1)
ESS <- (particle_set$ESS[1], (, (time_mesh)-1))
CESS <- (particle_set$CESS[1], (, (time_mesh)-1))
resampled <- (, (time_mesh))
if (adaptive_mesh) {
E_nu_j <- (, (time_mesh))
chosen <- ("", (time_mesh))
mesh_terms <- (((,)), (time_mesh))
k4_choice <- (, (time_mesh))
} {
E_nu_j <-
chosen <-
mesh_terms <-
k4_choice <-
}
# iterative proposals
end_time <- time_mesh[length(time_mesh)]
j <- 1
while (time_mesh[j]!=end_time) {
pcm <- ()
j <- j+1
# ----------- resample particle_set (only resample if ESS < N*ESS_threshold)
if (particle_set$ESS < N*ESS_threshold) {
resampled[j-1] <-
particle_set <- resample_particle_x_samples(N = N,
particle_set = particle_set,
multivariate = ,
= j-1,
resampling_method = resampling_method,
seed = seed)
} {
resampled[j-1] <-
}
# ----------- if adaptive_mesh==TRUE, find mesh for jth iteration
if (adaptive_mesh) {
if (particle_set$number_of_steps < j) {
particle_set$number_of_steps <- j
particle_set$CESS[j] <-
particle_set$resampled[j] <-
}
tilde_Delta_j <- mesh_guidance_adaptive( = m,
d = 2,
data_size = adaptive_mesh_parameters$data_size,
b = adaptive_mesh_parameters$b,
threshold = adaptive_mesh_parameters$CESS_j_threshold,
particle_set = particle_set,
sub_posterior_means = sub_posterior_means,
inv_precondition_matrices = inv_precondition_matrices,
k3 = adaptive_mesh_parameters$k3,
k4 = adaptive_mesh_parameters$k4,
vanilla = adaptive_mesh_parameters$vanilla)
E_nu_j[j] <- tilde_Delta_j$E_nu_j
chosen[j] <- tilde_Delta_j$chosen
mesh_terms[j]] <- (tilde_Delta_j$T1, tilde_Delta_j$T2)
k4_choice[j] <- tilde_Delta_j$k4_choice
time_mesh[j] <- (end_time, time_mesh[j-1]+tilde_Delta_j$max_delta_j)
}
# split the x samples from the previous time marginal (and their means) into approximately equal lists
split_x_samples <- (split_indices, (indices) particle_set$x_samples[indices])
split_x_means <- (split_indices, (indices) particle_set$x_means[indices,, = ])
V <- construct_V(s = time_mesh[j-1],
= time_mesh[j],
end_time = end_time,
= m,
d = 2,
precondition_matrices = precondition_matrices,
Lambda = Lambda,
iteration = j)
rho_j_weighted_samples <- ::(cl, X = 1:(split_indices), fun = (core) {
split_N <- (split_indices[core]])
x_mean_j <- ( = , = split_N, = 2)
log_rho_j <- (0, split_N)
x_j <- (1:split_N, (i) {
M <- construct_M(s = time_mesh[j-1],
= time_mesh[j],
end_time = end_time,
= m,
d = 2,
sub_posterior_samples = split_x_samples[core]][i]],
sub_posterior_mean = split_x_means[core]][i,],
iteration = j)
if (time_mesh[j]!=end_time) {
(((N = 1, mu = M, Sigma = V), = m, = 2, byrow = ))
} {
((mvtnorm::rmvnorm(n = 1, = M, = V), = m, = 2, byrow = ))
}
})
for (i 1:split_N) {
x_mean_j[i,] <- ( = x_j[i]],
= inv_precondition_matrices,
inverse_sum_weights = Lambda)
log_rho_j[i] <- ((1:m, () {
(x0 = (split_x_samples[core]][i]][c,]),
y = (x_j[i]][c,]),
s = time_mesh[j-1],
= time_mesh[j],
mean_vec = mean_vecs[c]],
sd_vec = sd_vecs[c]],
corr = corrs[c],
= betas[c],
precondition_mat = precondition_matrices[c]],
transform_mats = transform_matrices[c]],
diffusion_estimator = diffusion_estimator,
beta_NB = beta_NB,
gamma_NB_n_points = gamma_NB_n_points,
logarithm = )}))
}
(('x_j' = x_j, 'x_mean_j' = x_mean_j, 'log_rho_j' = log_rho_j))
})
# ---------- update particle set
# update the weights and return updated particle set
particle_set$x_samples <- ((1:(split_indices), (i) {
rho_j_weighted_samples[i]]$x_j}), recursive = )
particle_set$x_means <- (, (1:(split_indices), (i) {
rho_j_weighted_samples[i]]$x_mean_j}))
# update weight and normalise
log_rho_j <- ((1:(split_indices), (i) {
rho_j_weighted_samples[i]]$log_rho_j}))
norm_weights <- particle_ESS(log_weights = particle_set$log_weights + log_rho_j)
particle_set$log_weights <- norm_weights$log_weights
particle_set$normalised_weights <- norm_weights$normalised_weights
particle_set$ESS <- norm_weights$ESS
ESS[j] <- particle_set$ESS
# calculate the conditional ESS (i.e. the 1/sum(inc_change^2))
# where inc_change is the incremental change in weight (= log_rho_j)
particle_set$CESS[j] <- particle_ESS(log_weights = log_rho_j)$ESS
CESS[j] <- particle_set$CESS[j]
elapsed_time[j-1] <- (()-pcm)'elapsed']
}
::(cl)
# set the y samples as the first element of each of the x_samples
proposed_samples <- ((1:N, (i) particle_set$x_samples[i]][1,]))
particle_set$y_samples <- proposed_samples
# ----------- resample particle_set (only resample if ESS < N*ESS_threshold)
if (particle_set$ESS < N*ESS_threshold) {
resampled[particle_set$number_of_steps] <-
particle_set <- resample_particle_y_samples(N = N,
particle_set = particle_set,
multivariate = ,
resampling_method = resampling_method,
seed = seed)
} {
resampled[particle_set$number_of_steps] <-
}
if (adaptive_mesh) {
CESS <- CESS[1:j]
ESS <- ESS[1:j]
resampled <- resampled[1:j]
particle_set$time_mesh <- time_mesh[1:j]
elapsed_time <- elapsed_time[1:(j-1)]
E_nu_j <- E_nu_j[1:j]
chosen <- chosen[1:j]
mesh_terms <- mesh_terms[1:j]
k4_choice <- k4_choice[1:j]
}
(('particle_set' = particle_set,
'proposed_samples' = proposed_samples,
'time' = elapsed_time,
'ESS' = ESS,
'CESS' = CESS,
'resampled' = resampled,
'E_nu_j' = E_nu_j,
'chosen' = chosen,
'mesh_terms' = mesh_terms,
'k4_choice' = k4_choice))
}
#' Generalised Bayesian Fusion [parallel]
#'
#' Generalised Bayesian Fusion with bivariate Gaussian target
#'
#' @param particles_to_fuse list of length m, where particles_to_fuse[[c]]
#' contains the particles for the c-th sub-posterior
#' (a list of particles to fuse can be initialised by
#' initialise_particle_sets() function)
#' @param N number of samples
#' @param m number of sub-posteriors to combine
#' @param time_mesh time mesh used in Bayesian Fusion
#' @param precondition_matrices list of length m, where precondition_matrices[[c]]
#' is the precondition matrix for sub-posterior c
#' @param mean_vecs list of length m, where mean_vecs[[c]] is a vector of
#' length 2 for the mean of sub-posterior c
#' @param sd_vecs list of length m, where sd_vecs[[c]] is a vector of length 2
#' for the standard deviation of sub-posterior c
#' @param corrs vector of length m, where corrs[c] give the correlation value
#' between component 1 and component 2 for sub-posterior c
#' @param betas vector of length m, where betas[c] is the inverse temperature
#' value for c-th posterior
#' @param resampling_method method to be used in resampling, default is multinomial
#' resampling ('multi'). Other choices are stratified
#' resampling ('strat'), systematic resampling ('system'),
#' residual resampling ('resid')
#' @param ESS_threshold number between 0 and 1 defining the proportion
#' of the number of samples that ESS needs to be
#' lower than for resampling (i.e. resampling is carried
#' out only when ESS < N*ESS_threshold)
#' @param sub_posterior_means matrix with m rows and 2 columns, where sub_posterior_means[c,]
#' is the sub-posterior mean of sub-posterior c
#' @param adaptive_mesh logical value to indicate if an adaptive mesh is used
#' (default is FALSE)
#' @param adaptive_mesh_parameters list of parameters used for adaptive mesh
#' @param diffusion_estimator choice of unbiased estimator for the Exact Algorithm
#' between "Poisson" (default) for Poisson estimator
#' and "NB" for Negative Binomial estimator
#' @param beta_NB beta parameter for Negative Binomial estimator (default 10)
#' @param gamma_NB_n_points number of points used in the trapezoidal estimation
#' of the integral found in the mean of the negative
#' binomial estimator (default is 2)
#' @param seed seed number - default is NULL, meaning there is no seed
#' @param n_cores number of cores to use
#'
#' @return A list with components:
#' \describe{
#' \item{particles}{particles returned from fusion sampler}
#' \item{proposed_samples}{proposal samples from fusion sampler}
#' \item{time}{run-time of fusion sampler}
#' \item{elapsed_time}{elapsed time of each step of the algorithm}
#' \item{ESS}{effective sample size of the particles after each step}
#' \item{CESS}{conditional effective sample size of the particles after each step}
#' \item{resampled}{boolean value to indicate if particles were resampled
#' after each time step}
#' \item{E_nu_j}{approximation of the average variation of the trajectories
#' at each time step}
#' \item{chosen}{which term was chosen if using an adaptive mesh at each time step}
#' \item{mesh_terms}{the evaluated terms in deciding the mesh at each time step}
#' \item{precondition_matrices}{list of length 2 where precondition_matrices[[2]]
#' are the pre-conditioning matrices that were used
#' and precondition_matrices[[1]] are the combined
#' precondition matrices}
#' \item{sub_posterior_means}{list of length 2, where sub_posterior_means[[2]]
#' are the sub-posterior means that were used and
#' sub_posterior_means[[1]] are the combined
#' sub-posterior means}
#' }
#'
#' @export
<- (particles_to_fuse,
N,
m,
time_mesh,
mean_vecs,
sd_vecs,
corrs,
betas,
precondition_matrices,
resampling_method = 'multi',
ESS_threshold = 0.5,
diffusion_estimator = 'Poisson',
sub_posterior_means = ,
adaptive_mesh = ,
adaptive_mesh_parameters = ,
beta_NB = 10,
gamma_NB_n_points = 2,
seed = ,
n_cores = ::()) {
if (!(particles_to_fuse) | ((particles_to_fuse)!=m)) {
("parallel_GBF_biGaussian: particles_to_fuse must be a list of length m")
} if (!((particles_to_fuse, (sub_posterior) ("particle" %in% (sub_posterior))))) {
("parallel_GBF_biGaussian: particles in particles_to_fuse must be \"particle\" objects")
} if (!((particles_to_fuse, (sub_posterior) (sub_posterior$y_samples)))) {
("parallel_GBF_biGaussian: the particles' samples for y should all be matrices")
} if (!((particles_to_fuse, (sub_posterior) (sub_posterior$y_samples)==2))) {
("parallel_GBF_biGaussian: the particles' samples for y should all be matrices with 2 columns")
} if (!(time_mesh)) {
("parallel_GBF_biGaussian: time_mesh must be an ordered vector of length >= 2")
} if ((time_mesh) < 2) {
("parallel_GBF_biGaussian: time_mesh must be an ordered vector of length >= 2")
} if (!(time_mesh, (time_mesh))) {
("parallel_GBF_biGaussian: time_mesh must be an ordered vector of length >= 2")
} if (((1:m, () (!(mean_vecs[c]]) | ((mean_vecs[c]])!=2))))) {
("parallel_GBF_biGaussian: mean_vecs[[c]] must be a vector of length 2 for each c")
} if (((1:m, () (!(sd_vecs[c]]) | ((sd_vecs[c]])!=2))))) {
("parallel_GBF_biGaussian: sd_vecs[[c]] must be a vector of length 2 for each c")
} if (!(corrs) | ((corrs)!=m)) {
("parallel_GBF_biGaussian: corrs must be a vector of length m")
} if (!(betas) | ((betas)!=m)) {
("parallel_GBF_biGaussian: betas must be a vector of length m")
} if (!(precondition_matrices) | ((precondition_matrices)!=m)) {
("parallel_GBF_biGaussian: precondition_matrices must be a list of length m")
} if ((ESS_threshold < 0) | (ESS_threshold > 1)) {
("parallel_GBF_biGaussian: ESS_threshold must be between 0 and 1")
}
# set a seed if one is supplied
if (!(seed)) {
(seed)
}
# start time recording
pcm <- ()
# ---------- first importance sampling step
# pre-calculating the inverse precondition matrices
inv_precondition_matrices <- (precondition_matrices, )
Lambda <- (inv_precondition_matrices)
pcm_rho_0 <- ()
particles <- (particles_to_fuse = particles_to_fuse,
= 2,
N = N,
m = m,
= time_mesh[length(time_mesh)],
inv_precondition_matrices = inv_precondition_matrices,
inverse_sum_inv_precondition_matrices = Lambda,
number_of_steps = (time_mesh),
time_mesh = time_mesh,
resampling_method = resampling_method,
seed = seed,
n_cores = n_cores)
elapsed_time_rho_0 <- (()-pcm_rho_0)'elapsed']
# ---------- iterative steps
rho_j <- (particle_set = particles,
m = m,
time_mesh = time_mesh,
mean_vecs = mean_vecs,
sd_vecs = sd_vecs,
corrs = corrs,
betas = betas,
precondition_matrices = precondition_matrices,
inv_precondition_matrices = inv_precondition_matrices,
Lambda = Lambda,
resampling_method = resampling_method,
ESS_threshold = ESS_threshold,
sub_posterior_means = sub_posterior_means,
adaptive_mesh = adaptive_mesh,
adaptive_mesh_parameters = adaptive_mesh_parameters,
diffusion_estimator = diffusion_estimator,
beta_NB = beta_NB,
gamma_NB_n_points = gamma_NB_n_points,
seed = seed,
n_cores = n_cores)
if ((precondition_matrices, (((1, 2)), m))) {
new_precondition_matrices <- ((1, 2), precondition_matrices)
} {
new_precondition_matrices <- ((inv_precondition_matrices),
precondition_matrices)
}
if (!(sub_posterior_means)) {
new_sub_posterior_means <- (( = sub_posterior_means,
= inv_precondition_matrices,
inverse_sum_weights = Lambda),
sub_posterior_means)
} {
new_sub_posterior_means <- (, sub_posterior_means)
}
(('particles' = rho_j$particle_set,
'proposed_samples' = rho_j$proposed_samples,
'time' = (()-pcm)'elapsed'],
'elapsed_time' = (elapsed_time_rho_0, rho_j$),
'ESS' = rho_j$ESS,
'CESS' = rho_j$CESS,
'resampled' = rho_j$resampled,
'E_nu_j' = rho_j$E_nu_j,
'chosen' = rho_j$chosen,
'mesh_terms' = rho_j$mesh_terms,
'k4_choice' = rho_j$k4_choice,
'precondition_matrices' = new_precondition_matrices,
'sub_posterior_means' = new_sub_posterior_means))
}
