R/exact_algorithm_BLR.R

In rchan26/BayesLogitFusion: Monte Carlo Fusion for Bayesian Logistic Regression

#' @export
bound_phi_BLR <- function(initial_parameters,
                          y_labels,
                          X,
                          prior_means,
                          prior_variances,
                          C,
                          power,
                          lower,
                          upper,
                          precondition_mat) {
  # obtain lower bound using 'optim' function in R's stats package
  LB <- optim(par = initial_parameters,
              fn = function(beta) BayesLogitFusion:::phi_BLR(beta = beta,
                                                             y_labels = y_labels,
                                                             X = X,
                                                             prior_means = prior_means,
                                                             prior_variances = prior_variances,
                                                             C = C,
                                                             power = power,
                                                             precondition_mat = precondition_mat),
              method = "L-BFGS-B",
              lower = lower,
              upper = upper,
              control = list('fnscale' = 1,
                             'maxit' = 20))
  
  # obtain upper bound using 'optim' function in R's stats package
  UB <- optim(par = initial_parameters,
              fn = function(beta) BayesLogitFusion:::phi_BLR(beta = beta,
                                                             y_labels = y_labels,
                                                             X = X,
                                                             prior_means = prior_means,
                                                             prior_variances = prior_variances,
                                                             C = C,
                                                             power = power,
                                                             precondition_mat = precondition_mat),
              method = "L-BFGS-B",
              lower = lower,
              upper = upper,
              control = list('fnscale' = 1,
                             'maxit' = 20))
  
  return(list('LB' = LB$value,
              'UB' = UB$value))
}

#' @export
diffusion_probability_BLR <- function(dim,
                                      x0,
                                      y,
                                      s,
                                      t,
                                      K,
                                      initial_parameters,
                                      y_labels,
                                      X,
                                      prior_means,
                                      prior_variances,
                                      C,
                                      power, 
                                      precondition_mat) {
  # scale x0 and y by the pre-conditoning matrix
  x0 <- x0*diag(precondition_mat)
  y <- y*diag(precondition_mat)
  
  ##### simulate layer information to bound sample path
  # simulate layer information to bound sample path
  # create sequence vector
  a_seq <- seq(from = 0, to = ceiling(t-s), by = (t-s)/10)
  # simulate layer information (Bessel layer)
  bes_layers <- layeredBB::multi_bessel_layer_simulation(dim = dim,
                                                         x = x0,
                                                         y = y,
                                                         s = s,
                                                         t = t,
                                                         a = a_seq)
  
  ##### compute bounds for phi_BLR
  # set lower and upper bounds of the Brownian bridge
  lbound <- sapply(X = 1:dim, function(d) min(x0[d], y[d]) - bes_layers[[d]]$a[bes_layers[[d]]$l])
  ubound <- sapply(X = 1:dim, function(d) max(x0[d], y[d]) + bes_layers[[d]]$a[bes_layers[[d]]$l])
  
  # calculate upper and lower bounds of phi given the simulated sample layer information
  bounds <- bound_phi_BLR(initial_parameters = initial_parameters,
                          y_labels = y_labels,
                          X = X,
                          prior_means = prior_means,
                          prior_variances = prior_variances,
                          C = C,
                          power = power,
                          lower = lbound,
                          upper = ubound, 
                          precondition_mat = precondition_mat)
  LX <- bounds$LB
  UX <- bounds$UB
  
  ##### construct Poisson estimator
  # simulate skeletal points for Poisson thinning
  # simulate number of points to simulate from Poisson distribution
  kap <- rpois(1, (UX-LX)*(t-s))
  acc_prob <- 1
  
  # if kap > 0, then we simulate a layered Brownian bridge
  if (kap > 0) {
    # simulating Poisson times from a Uniform(0, (t-s)) distribution
    pois_times <- runif(kap, s, t)
    # simulating sample path at Poisson times
    layered_BB <- layeredBB::multi_layered_brownian_bridge(dim = dim,
                                                           x = x0,
                                                           y = y,
                                                           s = s,
                                                           t = t,
                                                           layers = bes_layers,
                                                           times = pois_times)
    
    # layered_bb includes points for times s and t - but we only want those simulated at pois_times
    pois_points <- as.matrix(layered_BB[1:(nrow(layered_BB)-1), 2:(ncol(layered_BB)-1)])
    
    # calculate acceptance probability
    for (j in 1:ncol(pois_points)) {
      phi_eval <- BayesLogitFusion:::phi_BLR(beta = pois_points[,j],
                                             y_labels = y_labels,
                                             X = X,
                                             prior_means = prior_means, 
                                             prior_variances = prior_variances, 
                                             C = C, 
                                             power = power,
                                             precondition_mat = precondition_mat)
      acc_prob <- acc_prob * ((UX-phi_eval)/(UX-LX))
    }
  }
  
  # return full acceptance probability
  return(acc_prob * exp(-(LX-K)*(t-s)))
}