R/game.R

Defines functions game

Documented in game

#' Gaussian Mixture Importance (GAME) sampling
#'
#' Implementation of the GAME algorithm for sampling of the marginal likelihood,
#' also called model evidence, for use in Bayesian model comparison and/or averaging
#' frameworks.
#'
#' @param theta \code{numeric}. An m-by-d matrix of posterior samples (e.g. generated by
#' \code{\link[HydroBayes]{dream_parallel}}).
#' @param theta_post \code{numeric}. An m-dimensional vector of posterior densities for
#' \code{theta} (e.g. generated by \code{\link[HydroBayes]{dream_parallel}}).
#' @param m_sub \code{integer}. The number of sub-samples from m in \code{theta} and
#' \code{theta_post} to be drawn to speed-up estimation of the Gaussian mixture distribution.
#' This sub-sample of m will be excluded when drawing the \code{m1} samples. Must be
#' \code{<= nrow(theta)-m1}. Default: NULL (all m samples are used).
#' @param m0 \code{integer}. Number of samples to be drawn from the mixture distribution
#' on wich the target function will be evaluated (an algorithmic parameter; only needed if \code{omega > 0};
#' also defines the number of times, \code{fun} needs to be evaluated and, thus, influences
#' the computing time, but also the sampling accuracy).
#' @param m1 \code{integer}. Number of sub-samples to be drawn from \code{theta} to
#' sample the model evidence. As such, it needs to be \code{<= nrow(theta)}.
#' @param Jmax \code{integer}. Maximum number of components to be used for the Gaussian
#' mixture distribution which is used as catalyst to sample the model evidence. Volpi
#' et al. (2017) suggest a value of 5 (default).
#' @param ic \code{character} string defining the Information Criterion to be used to
#' evaluate the calibration of the Gaussian mixture distribution. One of: 'bic' -
#' Bayesian Information Criterion; 'var' - variance of the ratio between the density of
#' the target distribution and the density of its approximation (the Gaussian mixture
#' distribution).
#' @param omega \code{numeric}. Parameter defining the sampling procedure to be applied:
#' 0: reciprocal importance sampling (RIS); 0 < omega < 1: bridge sampling with geometric
#' bridge (GB); 1: importance sampling (IS). Nedds to be within [0,1].
#' @param ob \code{logical}. If \code{TRUE}, bridge sampling with optimal bridge using
#' an initial value from GB or IS (depending on \code{omega}) will be applied. Default:
#' \code{FALSE}. If \code{omega == 0}, it will be ignored.
#' @param R \code{integer}. Number of iterations within bridge sampling with optimal
#' bridge. If \code{ob == TRUE} and \code{is.null(R)} (default), it will be set to
#' a value of 10 (suggested by Volpi et al., 2017). If \code{ob == FALSE}, it will be
#' ignored.
#' @param par.info \code{list}. See \code{\link[HydroBayes]{dream_parallel}}.
#' @param lik \code{integer}. See \code{\link[HydroBayes]{dream_parallel}}.
#' @param fun \code{character}. See \code{\link[HydroBayes]{dream_parallel}}.
#' @param ... Additional arguments to \code{fun}, see \code{\link[HydroBayes]{dream_parallel}}.
#' @param obs \code{numeric}. See \code{\link[HydroBayes]{dream_parallel}}.
#' @param abc_rho \code{character}. See \code{\link[HydroBayes]{dream_parallel}}.
#' @param abc_e \code{numeric}. See \code{\link[HydroBayes]{dream_parallel}}.
#' @param glue_shape \code{numeric}. See \code{\link[HydroBayes]{dream_parallel}}.
#' @param lik_fun \code{character}. See \code{\link[HydroBayes]{dream_parallel}}.
#' @param ncores \code{integer}. See \code{\link[HydroBayes]{dream_parallel}}.
#'
#' @return A named \code{list} with the following elements:
#'
#' \emph{Z}: Single value of the calculated marginal likelihood / Bayesian model evidence.
#'
#' \emph{J}: Number of components of the best calibrated Gaussian mixture distribution
#' (\code{<= Jmax}).
#'
#' \emph{ic_val}: Value of the selected Information Criterion of the best Gaussian mixture distribution.
#'
#' @note The arguments specific to function \code{\link[HydroBayes]{dream_parallel}}
#' are only needed if \code{omega > 0}, i.e. if additional evaluations of the target
#' function (\code{fun}) are necessary. Therefore, RIS might be most CPU efficient
#' (especially if an evaluation of \code{fun} requires a lot of computation time),
#' but could be inappropriate in sampling the model evidence!
#'
#' @details
#' Notation:
#'
#' - m: Number of MC samples.\cr
#' - d: number of parameters.\cr
#'
#' @references
#' Volpi, E., G. Schoups, G. Firmani and J. A. Vrugt: "Sworn testimony of the model
#' evidence: Gaussian Mixture Importance (GAME) sampling". Water Resources Research,
#' 2017, 53, \url{http://dx.doi.org/10.1002/2016WR020167}.
#'
#' @author Tobias Pilz \email{tpilz@@uni-potsdam.de}
#'
#' @export
game <- function(theta, theta_post, m_sub = NULL, m0 = NULL, m1, Jmax=5, ic, omega, ob = FALSE, R = NULL,
                 par.info = NULL, lik = NULL, fun = NULL, ..., obs = NULL, abc_rho = NULL,
                 abc_e = NULL, glue_shape = NULL, lik_fun = NULL, ncores = 1) {

  # argument checks
  if(nrow(theta) != length(theta_post)) stop("Number of rows of 'theta' and length of 'theta_post' must be the same!")
  if(omega < 0 | omega > 1) stop("Argument 'omega' must in [0,1]!")
  if(!grepl("bic|var", ic, ignore.case = T)) stop("Argument 'ic' must be one of {'bic', 'var'}!")
  if(m1 > nrow(theta)) stop("Argument 'm1' must not be larger than the number of rows of 'theta'!")
  if(!is.null(m_sub))
    if(m_sub > nrow(theta)-m1) stop("Argument 'm_sub' must be <= nrow(theta)-m1!")
  if(omega > 0) {
    if(is.null(m0)) stop("Argument 'omega' is > 0 and thus 'm0' needs to be given!")
    if(is.null(par.info)) stop("Argument 'omega' is > 0 and thus 'par.info' needs to be given!")
    if(is.null(lik)) stop("Argument 'omega' is > 0 and thus 'lik' needs to be given!")
    if(is.null(fun)) stop("Argument 'omega' is > 0 and thus 'fun' needs to be given!")
    if(ob & is.null(R)) R <- 10
  }
  if(ncores > 1)
    registerDoMC(ncores)

  # calibrate the mixture model for different J
  if(is.null(m_sub)) {
    m_subsamp <- 1:nrow(theta)
  } else {
    m_subsamp <- sample.int(nrow(theta), m_sub, replace = FALSE)
  }
  if(ncores > 1) {
    pmix_j <- mclapply(1:Jmax, function(j) em_pmix(theta[m_subsamp,], theta_post[m_subsamp], j, ic), mc.cores = ncores)
  } else {
    pmix_j <- lapply(1:Jmax, function(j) em_pmix(theta[m_subsamp,], theta_post[m_subsamp], j, ic))
  }

  # select optimal mixture distribution
  ic_vals <- sapply(pmix_j, function(x) x$ic_val)
  # sort out NaN and infinite values
  r_inf <- which(is.infinite(ic_vals) | is.nan(ic_vals))
  if(length(r_inf) == length(ic_vals))
    stop("Could not identify a suitable Gaussian mixture distribution. Change parameters 'Jmax' and/or 'ic' and try again. However, the distribution given within 'theta' might be not describable by a Gaussian mixture distribution!")
  ic_vals <- ic_vals[-r_inf]
  pmix_j <- pmix_j[-r_inf]
  pmix_optim <- pmix_j[[which.min(ic_vals)]]

  # draw sub-samples of posterior parameter realisations
  if(is.null(m_sub)) {
    r_samp <- sample.int(nrow(theta), size=m1, replace = FALSE)
  } else {
    r_samp <- sample((1:nrow(theta))[-m_subsamp], size=m1, replace = FALSE)
  }
  theta_m1 <- theta[r_samp, ]
  # get posterior densities for the sub-sample
  theta_m1_post <- theta_post[r_samp]
  # evaluate mixture density for the sub-samples
  theta_m1_pmix <- dmixmvn(theta_m1, pmix_optim$pmix)
  # avoid infinite values
  theta_m1_pmix <- pmin( pmax(theta_m1_pmix, .Machine$double.xmin), .Machine$double.xmax)

  # sampling of the evidence
  if(omega == 0) {
    # Reciprocal Importance Sampling (eq. 7 as special case of eq. 6)
    Z <- 1 / mean(theta_m1_pmix/theta_m1_post)

  } else if(omega > 0 & omega <= 1) {

    # draw samples from pmix (note that pmix is q0 in Volpi et al., 2017)
    sigma_full <- LTSigma2variance(pmix_optim$pmix$LTSigma) # get full covar matrix
    dist <- sample.int(pmix_optim$J, m0, replace = T, prob = pmix_optim$pmix$pi) # sampling from which of the J distributions?
    theta_m0_t <- t(sapply(dist, function(d) mvrnorm(1, pmix_optim$pmix$Mu[d,], sigma_full[,,d])))

    # check / bound parameters
    if(!is.null(par.info$min) && !is.null(par.info$max)) {
      theta_m0 <- apply(theta_m0_t, 1, function(xp) {
        sapply(1:length(xp), function(k) bound_par(xp[k], min = ifelse(length(par.info$min) > 1, par.info$min[k], par.info$min),
                                                max = ifelse(length(par.info$max) > 1, par.info$max[k], par.info$max),
                                                handle = par.info$bound))
      })
      theta_m0 <- t(theta_m0)
    }
    # make sure parameter names are retained
    if(!is.null(par.info$names)) colnames(theta_m0) <- par.info$names

    # evaluate mixture density for the samples
    theta_m0_pmix <- dmixmvn(theta_m0, pmix_optim$pmix)
    # avoid infinite values
    theta_m0_pmix <- pmin( pmax(theta_m0_pmix, .Machine$double.xmin), .Machine$double.xmax)

    # evaluate target density for the samples
    # calculate prior log-density
    lp_xp <- apply(theta_m0, 1, prior_pdf, par.info, lik)
    # evaluate fun for proposals
    if(ncores > 1) {
      res_fun_t <- mclapply(1:m0, function(j) get(fun)(theta_m0[j,], ...), mc.cores = ncores)
      res_fun_t <- matrix(unlist(res_fun_t), nrow=m0, byrow = T) # m0-by-[length of fun output]
    } else {
      res_fun_t <- apply(theta_m0, 1, get(fun), ...)
      if(!is.matrix(res_fun_t)) res_fun_t <- t(res_fun_t)
      res_fun_t <- t(res_fun_t) # m0-by-[length of fun output]
    }
    # calculate log-likelihood
    ll_xp <- apply(res_fun_t, 1, calc_ll, lik, obs, abc_rho, abc_e, glue_shape, lik_fun)
    # calculate posterior log-density
    lpost_xp <- ll_xp + lp_xp
    # normal space
    theta_m0_post <- exp(lpost_xp)

    if(omega == 1) {
      # Importance Sampling (eq. 10 as special case of eq. 6)
      Z <- mean(theta_m0_post/theta_m0_pmix)
    } else {
      # Bridge Sampling with Geometric Bridge (eq. 6)
      q_half_m0 <- theta_m0_pmix^(1-omega) * theta_m0_post^omega
      q_half_m1 <- theta_m1_pmix^(1-omega) * theta_m1_post^omega
      Z <- mean(q_half_m0/theta_m0_pmix) / mean(q_half_m1/theta_m1_post)
    }
    if(ob) {
      # Bridge Sampling with Optimal Bridge (eq. 11)
      s0 <- m0 / (m0 + m1)
      s1 <- m1 / (m0 + m1)
      r=1
      while(r<=R){
        Z <- mean( (theta_m0_post/theta_m0_pmix) / (s0*Z + s1*theta_m0_post/theta_m0_pmix) ) / mean( 1 / (s0*Z + s1*theta_m1_post/theta_m1_pmix) )
        r <- r+1
      }
    } # ob
  } # omega > 0 & omega <= 1

  # return output
  return(list(Z = Z,
              J = pmix_optim$J,
              ic_val = pmix_optim$ic_val))
} # EOF
tpilz/HydroBayes documentation built on May 6, 2019, 3:44 p.m.