R/reflectionmaximal4_pcn_synchronous_coupling.R

In jeremyhengjm/UnbiasedMultilevel: Unbiased estimators for discretized models

#' @rdname reflectionmaximal4_pcn_synchronous_coupling
#' @title Reflection maximal coupling of four pCN proposals
#' @description Generate pCN proposals for four chains via reflection maximal coupling of all four chains
#' @param chain_state_coarse1 a vector with coordinates of the first particle at the coarse level
#' @param chain_state_coarse2 a vector with coordinates of the second particle at the coarse level
#' @param chain_state_fine1 a vector with coordinates of the first particle at the fine level
#' @param chain_state_fine2 a vector with coordinates of the second particle at the fine level
#' @param identical_coarse a flag that is True if chains at the coarse level are identical (coincide) and False otherwise
#' @param identical_fine a flag that is True if chains at the fine level are identical (coincide) and False otherwise
#' @param tuning_coarse a list that contains parameters needed for genertaion of pCN proposal at the coarse level: standard devation and rho
#' @param tuning_coarse a list that contains parameters needed for genertaion of pCN proposal at the fine level: standard devation and rho
#'@return a list that contains state of the first chain at the coarse level, state of the second chain at the coarse level, state of the second chain, updated value of the flag "identical", cost of proposal generation
#'@export

reflectionmaximal4_pcn_synchronous_coupling <- function(chain_state_coarse1, chain_state_coarse2,
                                                        chain_state_fine1, chain_state_fine2,
                                                        identical_coarse, identical_fine,
                                                        tuning_coarse, tuning_fine){
  cost = 0  # cost of proposal generation
  # extract tuning parameters
  proposal_sd_coarse <- tuning_coarse$proposal_sd
  proposal_rho_coarse <- tuning_coarse$proposal_rho
  proposal_sd_factor_coarse <- sqrt(1-proposal_rho_coarse^2) * proposal_sd_coarse
  
  proposal_sd_fine <- tuning_fine$proposal_sd
  proposal_rho_fine <- tuning_fine$proposal_rho
  proposal_sd_factor_fine <- sqrt(1-proposal_rho_fine^2) * proposal_sd_fine
  
  ## sample 2-way maximal coupling for coarse chain
  # sample first proposal
  randn_coarse1 <- rnorm(dimension)
  coarse1 <- proposal_rho_coarse * chain_state_coarse1 + proposal_sd_factor_coarse * randn_coarse1
  
  # difference
  zdiff_coarse <- proposal_rho_coarse * (chain_state_coarse1 - chain_state_coarse2) / proposal_sd_factor_coarse
  
  # evaluate proposal transition densities at first proposal
  pdf1 <- sum(dnorm(randn_coarse1, log = TRUE))
  pdf2 <- sum(dnorm(randn_coarse1 + zdiff_coarse, log = TRUE))
  logacceptprob <- min(pdf1, pdf2) - pdf1
  loguniform <- log(runif(1))
  if (loguniform < logacceptprob){
    # return common proposal for both chains
    coarse2 <- coarse1
    identical_coarse <- TRUE
    
  } else {
    # perform reflection for second proposal
    evec_coarse <- zdiff_coarse / sqrt(sum(zdiff_coarse^2))
    randn_coarse2 <- randn_coarse1 - 2 * sum(evec_coarse * randn_coarse1) * evec_coarse
    coarse2 <- proposal_rho_coarse * chain_state_coarse2 + proposal_sd_factor_coarse * randn_coarse2
    identical_coarse <- FALSE
    
  }
  
  ## sample 2-way maximal coupling for fine chain
  # sample first proposal
  randn_fine1 <- randn_coarse1 # synchronous coupling
  fine1 <- proposal_rho_fine * chain_state_fine1 + proposal_sd_factor_fine * randn_fine1
  
  # difference
  zdiff_fine <- proposal_rho_fine * (chain_state_fine1 - chain_state_fine2) / proposal_sd_factor_fine
  
  # evaluate proposal transition densities at first proposal
  pdf2 <- sum(dnorm(randn_fine1 + zdiff_fine, log = TRUE))
  logacceptprob <- min(pdf1, pdf2) - pdf1
  
  if (loguniform < logacceptprob){
    # return common proposal for both chains
    fine2 <- fine1
    identical_fine <- TRUE
    
  } else {
    # perform reflection for second proposal
    evec_fine <- zdiff_fine / sqrt(sum(zdiff_fine^2))
    randn_fine2 <- randn_fine1 - 2 * sum(evec_fine * randn_fine1) * evec_fine
    fine2 <- proposal_rho_fine * chain_state_fine2 + proposal_sd_factor_fine * randn_fine2
    identical_fine <- FALSE
    
  }
  
  return(list(coarse1 = coarse1, coarse2 = coarse2, fine1 = fine1, fine2 = fine2,
              identical_coarse = identical_coarse, identical_fine = identical_fine, cost = cost))
  
}