R/mcmc.r

Defines functions mcmcMh

Documented in mcmcMh

#' Metropolis-Hastings MCMC
#'
#' Run \code{nIterations} of a Metropolis-Hastings MCMC to sample from the
#' target distribution using a gaussian proposal kernel.
#' Two optional optimizations are also implemented: truncated gaussian proposal
#' (to match the support of the target distribution, i.e. boundary of the
#' parameters) and adaptive gaussian proposal (to match the size and the shape
#' of the target distribution).
#' @param target \R-function that takes a single argument: \code{theta} (named
#'   numeric vector of parameter values) and returns a list of 2 elements:
#' \itemize{
#' \item \code{logDensity} the logged value of the target density, evaluated at
#'   \code{theta}.
#' \item \code{trace} a named numeric vector of values to be printed in the
#'   \code{trace} data.frame returned by \code{mcmcMh}.
#' }
#' @param initTheta named vector of initial parameter values to start the
#'   chain.
#' @param proposalSd vector of standard deviations. If this is given and covmat
#'   is not, a diagonal matrix will be built from this to use as covariance
#'   matrix of the multivariate Gaussian proposal distribution. By default, this
#'   is set to \code{initTheta/10}.
#' @param nIterations number of iterations to run the MCMC chain.
#' @param covmat named numeric covariance matrix of the multivariate Gaussian
#'   proposal distribution. Must have named rows and columns with at least all
#'   estimated theta. If \code{proposalSd} is given, this is ignored.
#' @param limits limits for the - potentially truncated - multi-variate normal
#'   proposal distribution of the MCMC. Contains 2 elements:
#' \itemize{
#'      \item \code{lower} named numeric vector. Lower truncation points in each
#'   dimension of the Gaussian proposal distribution. By default they are set to
#'   \code{-Inf}.
#'      \item \code{upper} named numeric vector. Upper truncation points in each
#'   dimension of the Gaussian proposal distribution. By default they are set to
#'   \code{Inf}.
#' }
#' @param adaptSizeStart number of iterations to run before adapting the size
#'   of the proposal covariance matrix (see note below). Set to NULL (default)
#'   if size is not to be adapted.
#' @param adaptSizeCooling cooling factor for the scaling factor of the
#'   covariance matrix during size adaptation (see note below).
#' @param adaptShapeStart number of accepted jumps before adapting the shape of
#'   the proposal covariance matrix (see note below). Set to NULL (default) if
#'   shape is not to be adapted
#' @param adaptShapeStop number of iterations to run with adaptations
#'   of the shape of the proposal covariance matrix  before stopping. Se
#'   to NULL (default) if never to stop.
#' @param printInfoEvery frequency of information on the chain: acceptance
#'   rate and state of the chain. Default value to \code{nIterations/100}. Set
#'   to \code{NULL} to avoid any info.
#' @param verbose logical. If \code{TRUE}, information are printed.
#' @param maxScalingSd numeric. Maximum value for the scaling factor of the
#'   covariance matrix. Avoid too high values for the scaling factor, which
#'   might happen due to the exponential update scheme. In this case, the
#'   covariance matrix becomes too wide and the sampling from the truncated
#'   proposal kernel becomes highly inefficient
#' @note The size of the proposal covariance matrix is adapted using the
#'   following formulae: \deqn{\Sigma_{n+1}=\sigma_n * \Sigma_n} with
#'   \eqn{\sigma_n=\sigma_{n-1}*exp(\alpha^n*(acc - 0.234))}, where \eqn{\alpha}
#'   is equal to \code{adaptSizeCooling} and \eqn{acc} is the acceptance rate
#'   of the chain.
#'
#' The shape of the proposal covariance matrix is adapted using the following
#'   formulae: \deqn{\Sigma_{n+1}=2.38^2/d * \Sigma_n} with \eqn{\Sigma_n} the
#'   empirical covariance matrix and \eqn{d} is the number of estimated
#'   parameters in the model.
#' @references Roberts GO, Rosenthal JS. Examples of adaptive MCMC. Journal of
#'   Computational and Graphical Statistics. Taylor & Francis;
#'   2009;18(2):349-67.
#' @export
#' @importFrom lubridate as.period
#' @importFrom tmvtnorm rtmvnorm dtmvnorm
#' @importFrom stats runif
#' @return a list with 3 elements:
#' \itemize{
#'      \item \code{trace} a \code{data.frame}. Each row contains a state of the
#'   chain (as returned by \code{target}, and an extra column for the
#'   logDensity).
#'      \item \code{acceptanceRate} acceptance rate of the MCMC chain.
#'      \item \code{covmatProposal} last covariance matrix used for proposals.
#' }
# nolint start: cyclocomp_linter
mcmcMh <- function(
    target, initTheta, proposalSd = NULL,
    nIterations, covmat = NULL,
    limits = list(lower = NULL, upper = NULL),
    adaptSizeStart = NULL, adaptSizeCooling = 0.99,
    adaptShapeStart = NULL, adaptShapeStop = NULL,
    printInfoEvery = nIterations / 100,
    verbose = FALSE, maxScalingSd = 50) {
  # initialise theta
  thetaCurrent <- initTheta
  thetaPropose <- initTheta

  # extract theta of gaussian proposal
  covmatProposal <- covmat
  lowerProposal <- limits$lower
  upperProposal <- limits$upper

  # reorder vector and matrix by names, set to default if necessary
  thetaNames <- names(initTheta)
  if (!is.null(proposalSd) && is.null(names(proposalSd))) {
    names(proposalSd) <- thetaNames
  }

  if (is.null(covmatProposal)) {
    if (is.null(proposalSd)) {
      proposalSd <- initTheta / 10
    }
    covmatProposal <-
      matrix(diag(proposalSd[thetaNames]^2, nrow = length(thetaNames)),
        nrow = length(thetaNames),
        dimnames = list(thetaNames, thetaNames)
      )
  } else {
    covmatProposal <- covmatProposal[thetaNames, thetaNames]
  }

  if (is.null(lowerProposal)) {
    lowerProposal <- initTheta
    lowerProposal[] <- -Inf
  } else {
    lowerProposal <- lowerProposal[thetaNames]
  }

  if (is.null(upperProposal)) {
    upperProposal <- initTheta
    upperProposal[] <- Inf
  } else {
    upperProposal <- upperProposal[thetaNames]
  }

  # covmat init
  covmatProposalInit <- covmatProposal

  adaptingSize <- FALSE # will be set to TRUE once we start
  # adapting the size

  adaptingShape <- 0 # will be set to the iteration at which
  # adaptation starts

  # find estimated theta
  thetaEstimatedNames <- names(which(diag(covmatProposal) > 0))

  # evaluate target at theta init
  targetThetaCurrent <- target(thetaCurrent)

  # if return value is a vector, set log.density and trace
  if (is(targetThetaCurrent, "numeric")) {
    targetThetaCurrent <- list(
      logDensity = targetThetaCurrent,
      trace = thetaCurrent
    )
  }

  if (!is.null(printInfoEvery)) {
    message(
      Sys.time(), ", Init: ",
      printNamedVector(thetaCurrent[thetaEstimatedNames]),
      ", target: ", targetThetaCurrent[["logDensity"]]
    )
  }

  # trace
  trace <- matrix(
    ncol = length(targetThetaCurrent[["trace"]]) + 1,
    nrow = nIterations,
    0
  )
  colnames(trace) <- c(names(targetThetaCurrent[["trace"]]), "logDensity")

  # acceptance rate
  acceptanceRate <- 0

  # scaling factor for covmat size
  scalingSd <- 1

  # scaling multiplier
  scalingMultiplier <- 1

  # empirical covariance matrix (0 everywhere initially)
  covmatEmpirical <- covmatProposal
  covmatEmpirical[, ] <- 0

  # empirical mean vector
  thetaMean <- thetaCurrent

  # if printInfoEvery is null never print info
  if (is.null(printInfoEvery)) {
    printInfoEvery <- nIterations + 1
  }

  for (iIteration in seq_len(nIterations)) {
    # adaptive step
    if (!is.null(adaptSizeStart) && iIteration >= adaptSizeStart &&
      (is.null(adaptShapeStart) ||
       acceptanceRate * iIteration < adaptShapeStart)) {
      if (!adaptingSize) {
        message("\n---> Start adapting size of covariance matrix")
        adaptingSize <- TRUE
      }
      # adapt size of covmat until we get enough accepted jumps
      scalingMultiplier <-
        exp(adaptSizeCooling^(iIteration - adaptSizeStart) *
            (acceptanceRate - 0.234))
      scalingSd <- scalingSd * scalingMultiplier
      scalingSd <- min(c(scalingSd, maxScalingSd))
      # only scale if it doesn't reduce the covariance matrix to 0
      covmatProposalNew <- scalingSd^2 * covmatProposalInit
      if (!(any(diag(covmatProposalNew)[thetaEstimatedNames] <
        .Machine$double.eps))) {
        covmatProposal <- covmatProposalNew
      }
    } else if (!is.null(adaptShapeStart) &&
      acceptanceRate * iIteration >= adaptShapeStart &&
      (adaptingShape == 0 || is.null(adaptShapeStop) ||
        iIteration < adaptingShape + adaptShapeStop)) {
      if (!adaptingShape) {
        message("\n---> Start adapting shape of covariance matrix")
        adaptingShape <- iIteration
      }

      ## adapt shape of covmat using optimal scaling factor for multivariate
      ## target distributions
      scalingSd <- 2.38 / sqrt(length(thetaEstimatedNames))

      covmatProposal <- scalingSd^2 * covmatEmpirical
    } else if (adaptingShape > 0) {
      message("\n---> Stop adapting shape of covariance matrix")
      adaptingShape <- -1
    }

    # print info
    if (iIteration %% ceiling(printInfoEvery) == 0) {
      message(Sys.time(), ", Iteration: ", iIteration, "/", nIterations,
        ", acceptance rate: ",
        sprintf("%.3f", acceptanceRate),
        appendLF = FALSE
      )
      if (!is.null(adaptSizeStart) || !is.null(adaptShapeStart)) {
        message(", scalingSd: ", sprintf("%.3f", scalingSd),
          ", scalingMultiplier: ", sprintf("%.3f", scalingMultiplier),
          appendLF = FALSE
        )
      }
      message(", state: ", printNamedVector(targetThetaCurrent[["trace"]]))
      message(", logdensity: ", targetThetaCurrent[["logDensity"]])
    }

    # propose another parameter set
    if (any(diag(covmatProposal)[thetaEstimatedNames] <
      .Machine$double.eps)) {
      print(covmatProposal[thetaEstimatedNames, thetaEstimatedNames])
      stop("non-positive definite covmat", call. = FALSE)
    }
    if (length(thetaEstimatedNames) > 0) {
      thetaPropose[thetaEstimatedNames] <-
        as.vector(rtmvnorm(1,
          mean =
            thetaCurrent[thetaEstimatedNames],
          sigma =
            covmatProposal[thetaEstimatedNames, thetaEstimatedNames],
          lower =
            lowerProposal[thetaEstimatedNames],
          upper = upperProposal[thetaEstimatedNames]
        ))
    }

    # evaluate posterior of proposed parameter
    targetThetaPropose <- target(thetaPropose)
    # if return value is a vector, set logDensity and trace
    if (is(targetThetaPropose, "numeric")) {
      targetThetaPropose <- list(
        logDensity = targetThetaPropose,
        trace = thetaPropose
      )
    }


    if (!is.finite(targetThetaPropose[["logDensity"]])) {
      # if posterior is 0 then do not compute anything else and don't accept
      logAcceptance <- -Inf
    } else {
      # compute Metropolis-Hastings ratio (acceptance probability)
      logAcceptance <- targetThetaPropose[["logDensity"]] -
        targetThetaCurrent[["logDensity"]]
      logAcceptance <- logAcceptance +
        dtmvnorm(
          x = thetaCurrent[thetaEstimatedNames],
          mean =
            thetaPropose[thetaEstimatedNames],
          sigma =
            covmatProposal[
              thetaEstimatedNames,
              thetaEstimatedNames
            ],
          lower =
            lowerProposal[thetaEstimatedNames],
          upper =
            upperProposal[thetaEstimatedNames],
          log = TRUE
        )
      logAcceptance <- logAcceptance -
        dtmvnorm(
          x = thetaPropose[thetaEstimatedNames],
          mean = thetaCurrent[thetaEstimatedNames],
          sigma =
            covmatProposal[
              thetaEstimatedNames,
              thetaEstimatedNames
            ],
          lower =
            lowerProposal[thetaEstimatedNames],
          upper =
            upperProposal[thetaEstimatedNames],
          log = TRUE
        )
    }

    if (verbose) {
      message("Propose: ", thetaPropose[thetaEstimatedNames],
        ", target: ", targetThetaPropose,
        ", acc prob: ", exp(logAcceptance), ", ",
        appendLF = FALSE
      )
    }

    if (isAccepted <- (log(runif(1)) < logAcceptance)) {
      # accept proposed parameter set
      thetaCurrent <- thetaPropose
      targetThetaCurrent <- targetThetaPropose
      if (verbose) {
        message("accepted")
      }
    } else if (verbose) {
      message("rejected")
    }
    trace[iIteration, ] <- c(
      targetThetaCurrent[["trace"]], targetThetaCurrent[["logDensity"]]
    )

    # update acceptance rate
    if (iIteration == 1) {
      acceptanceRate <- isAccepted
    } else {
      acceptanceRate <- acceptanceRate +
        (isAccepted - acceptanceRate) / iIteration
    }

    # update empirical covariance matrix
    if (adaptingShape >= 0) {
      tmp <- updateCovmat(
        covmatEmpirical, thetaMean,
        thetaCurrent, iIteration
      )
      covmatEmpirical <- tmp$covmat
      thetaMean <- tmp$thetaMean
    }
  }

  return(list(
    trace = trace,
    acceptanceRate = acceptanceRate,
    covmatEmpirical = covmatEmpirical
  ))
}
# nolint end: cyclocomp_linter
sbfnk/fitR documentation built on July 18, 2023, 3:28 p.m.