R/c_chains_to_dataframe.R

In pierrejacob/debiasedmcmc: Unbiased MCMC with couplings

#'@rdname c_chains_to_dataframe
#'@title Obtain data frame representation of measure from list of coupled chains
#'@description From coupled chains,
#' presumably generated via \code{\link{sample_coupled_chains}},
#' and a choice of integers k and m, the function constructs
#' a data frame representation of an empirical signed measure, i.e.
#'
#'  \deqn{\hat{\pi}(dx) = \sum_{n=1}^N \omega_n \delta_{Z_n}(dx)}
#'
#' The function returns a data frame with first column "rep" indicating index of coupled chain,
#' second column "MCMC" indicating whether atom is part of the "MCMC" part of the signed measure (1) or the bias correction part (0)
#' third column "weight" indicating weights,
#' remaining columns "atom.1", "atom.2", etc containing components of the atoms
#'
#'@param c_chains A list containing coupled chains generated by \code{\link{sample_coupled_chains}}.
#'@param k An integer at which to start computing the signed measure; should be less than m
#'@param m A time horizon: should be less than the length of the chains; typically the same
#'m that was used in the call to \code{\link{sample_coupled_chains}}, or a smaller value
#'@param dopar Boolean (default to FALSE) indicating whether to perform calculation using registered parallel cores
#'@return A data frame
#'@export
c_chains_to_dataframe <- function(c_chains, k, m, dopar = FALSE, prune = TRUE){
  if ("meetingtime" %in% names(c_chains)){
    # the user supplied a coupled chain, rather than a list of coupled chain
    c_chains <- list(c_chains)
  }
  # number of coupled chains
  nsamples <- length(c_chains)
  approximation <- NULL
  if (dopar){
    # get signed measure representation of each coupled chain
    approximation <-  foreach(irep = 1:nsamples, .combine = rbind) %dopar% {
      ms_irep <- c_chains_to_measure_as_list(c_chains[[irep]], k, m)
      data.frame(rep = rep(irep, length(ms_irep$weights)), MCMC = ms_irep$MCMC, weight = ms_irep$weights, atom = ms_irep$atoms)
    }
  } else {
    approximation <- data.frame()
    for (irep in 1:nsamples){
      ms_irep <- c_chains_to_measure_as_list(c_chains[[irep]], k, m)
      approximation <- rbind(approximation,
         data.frame(rep = rep(irep, length(ms_irep$weights)), MCMC = ms_irep$MCMC, weight = ms_irep$weights, atom = ms_irep$atoms))
    }
  }
  # normalize the weights
  approximation$weight <- approximation$weight / nsamples
  # prune identical atoms
  if (prune){
    approximation <- data.frame(prune_measure_(as.matrix(approximation[do.call(order, as.list(approximation[,c(1,2,4:ncol(approximation)),drop=F])),])))
  }
  names(approximation) <- c("rep", "MCMC", "weight", paste0("atom.", 1:(ncol(approximation)-3)))
  approximation <- approximation[do.call(order, as.list(approximation[,1:2])),]
  return(approximation)
}