R/evmSim.R

Defines functions initRNG texmexJumpConst texmexCheckMap evmSim

Documented in evmSim

#' MCMC simulation around an evmOpt fit
#'
#' @param o a fit \code{evmOpt} object
#' @param priorParameters A list with two components. The first should
#'     be a vector of means, the second should be a covariance matrix
#'     if the penalty/prior is "gaussian" or "quadratic" and a
#'     diagonal precision matrix if the penalty/prior is "lasso", "L1"
#'     or "Laplace".  If \code{method = "simulate"} then these
#'     represent the parameters in the Gaussian prior distribution.
#'     If \code{method = 'optimize'} then these represent the
#'     parameters in the penalty function.  If not supplied: all
#'     default prior means are zero; all default prior variances are
#'     \eqn{10^4}; all covariances are zero.
#' @param prop.dist The proposal distribution to use, either
#'     multivariate gaussian or a multivariate Cauchy.
#' @param jump.const Control parameter for the Metropolis algorithm.
#' @param jump.cov Covariance matrix for proposal distribution of
#'     Metropolis algorithm.  This is scaled by \code{jump.const}.
#' @param verbose Whether or not to print progress to screen. Defaults
#'     to \code{verbose=TRUE}.
#' @param iter Number of simulations to generate
#' @param start Starting values for the chain; if missing, defaults to
#'     the MAP/ML estimates in \code{o}.
#' @param burn The number of initial steps to be discarded.
#' @param thin The degree of thinning of the resulting Markov chains.
#' @param chains The number of Markov chains to run. Defaults to 1. If you run
#'   more, the function will try to figure out how to do it in parallel using
#'   the same number of cores as chains.
#' @param export Character vector of names of variables to export. See the
#'   help file for \code{parallel::export}. Defaults to \code{export = NULL}
#'   and most users will never need to use it. Only matters on Windows.
#' @param trace How frequently to talk to the user
#' @param theCall (internal use only)
#' @param ... ignored
#' @return an object of class \code{evmSim}:
#'
#' \item{call}{The call to \code{evmSim} that produced the object.}
#'
#' \item{threshold}{The threshold above which the model was fit.}
#'
#' \item{map}{The point estimates found by maximum penalized
#' likelihood and which were used as the starting point for the Markov
#' chain.  This is of class \code{evmOpt} and methods for this class
#' (such as resid and plot) may be useful.}
#'
#' \item{burn}{The number of steps of the Markov chain that are to be
#' treated as the burn-in and not used in inferences.}
#'
#' \item{thin}{The degree of thinning used.}
#'
#' \item{chains}{The entire Markov chain generated by the Metropolis
#' algorithm.}
#'
#' \item{y}{The response data above the threshold for
#' fitting.}
#'
#' \item{seed}{The seed used by the random number generator.}
#'
#' \item{param}{The remainder of the chain after deleting
#' the burn-in and applying any thinning.}
#' @note it is not expected that the user should call this directly
#' @export
evmSim <- function(o, priorParameters, prop.dist,
                   jump.const, jump.cov, iter, start,
                   thin, burn, chains, export = NULL,
                   verbose, trace, theCall, ...){
    if (!inherits(o, "evmOpt")){
        stop("o must be of class 'evmOpt'")
    }
    # Run checks and initialize algorithm
    wh <- texmexCheckMap(o)
    jump.const <- texmexJumpConst(jump.const, o)
    seed <- initRNG()
    cov <- if (missing(jump.cov)) { o$cov } else { jump.cov }

    # Initialize matrix to hold chain
    res <- matrix(ncol=length(o$coefficients), nrow=iter)
    res[1,] <- if (is.null(start)) { o$coefficients } else { start }

    ############################# Get prior and log-likelihood
    prior <- .make.mvn.prior(priorParameters)

    evm.log.lik <- o$family$log.lik(o$data, th=o$threshold, ...)

    log.lik <- function(param) {
      evm.log.lik(param) + prior(param)
    }

    # create proposals en bloc
    proposal.fn <- switch(prop.dist,
                          gaussian=rmvnorm,
                          cauchy=.rmvcauchy,
                          function () {stop("Bad proposal distribution")})

    proposals <- proposal.fn(iter,
                             double(length(o$coefficients)),
                             cov*jump.const)

    if (chains > 1){
      # Set up parallel chains depending in OS
      os <- .Platform$OS.type
      if (os == "windows"){
        getCluster <- function(n){
          wh <- try(requireNamespace("parallel"))
          if (!inherits(wh, "try-error")){
            if (is.null(n)) n <- parallel::detectCores()
            if (n == 1) { NULL }
            else parallel::makeCluster(n)
          }
          else NULL
        }
        cluster <- getCluster(min(c(chains, parallel::detectCores() - 1)))
        on.exit(if (!is.null(cluster)){ parallel::stopCluster(cluster) })

        if (!is.null(cluster)){
          if (!is.null(export)){
            parallel::clusterExport(cl = cluster, varlist = export)
          }
        }
      }
    } else {
      os <- "1chain"
    }

    # Seeds required by parallel::parLapply
    seeds <- as.integer(runif(chains, -(2^31 - 1), 2^31))

    ######################## Run the Metropolis algorithm...
    if (os == "unix"){
      res <- parallel::mclapply(1:chains, texmexMetropolis, o = res, log.lik = log.lik,
                                proposals = proposals, verbose = verbose,  trace = trace,
                                seeds = seeds, mc.cores = chains)
    } else if (os == "windows") {
      res <- parallel::parLapply(cluster, X=1:chains,
                                 texmexMetropolis, o = res, log.lik = log.lik,
                                 proposals = proposals, verbose = verbose,  trace = trace,
                                 seeds = seeds)
    } else { # 1 chain or something not handled
      res <- list(texmexMetropolis(1, res, log.lik, proposals, verbose, trace, seeds))
    }

    if (inherits(res[[1]], "try-error")){
      stop("Something went wrong: simulations not performed. Find this stop message in evmSim")
    }

    # XXX Tidy up the object below. Doesn't need any of the info in o, or acceptance
    res <- list(call=theCall, map = o,
                burn = burn, thin = thin,
                chains=res, seeds=seeds)

    oldClass(res) <- "evmSim"
    res <- thinAndBurn(res)
    res
}


texmexMetropolis <-
    # Metropolis algorithm.
    # x is a matrix, initialized to hold the chain. It's first row should be the
    #    starting point of the chain.
    # proposals is a matrix of proposals
function(X, o, log.lik, proposals, verbose, trace, seeds){
    s <- set.seed(seeds[[X]])

    last.cost <- log.lik(o[1,])
    if (!is.finite(last.cost)) {
      stop("Start is infeasible.")
    }

    acc <- 0
    for(i in 2:nrow(o)){
      if( verbose){
        if(i %% trace == 0) message(i, " steps taken" )
      }
      prop <- proposals[i - 1,] + o[i - 1,]
      top <- log.lik(prop)
      delta <- top - last.cost
      if (is.finite(top) && ((delta >= 0) ||
                             (runif(1) <= exp(delta)))) {
        o[i, ] <- prop
        last.cost <- top
        acc <- 1 + acc
      } else {
        o[i, ] <- o[i - 1,]
      }
    } # Close for(i in 2:nrow

    acc <- acc / nrow(o)
    if (acc < .1) {
        warning("Acceptance rate in Metropolis algorithm is low.")
    }
    if ((trace < nrow(o)) & verbose) {
        message("Acceptance rate:", round(acc , 3) , "\n")
    }

    attr(o, "acceptance") <- acc
    o
}

####### Support functions

# Need to check for convergence failure here. Otherwise, end up simulating
# proposals from distribution with zero variance in 1 dimension.
texmexCheckMap <- function(map){
    checkNA <- any(is.na(sqrt(diag(map$cov))))
    if (checkNA) {
        stop("MAP estimates have not converged or have converged on values for which the variance cannot be computed. Cannot proceed. Try a different prior" )
    }
    NULL
}

texmexJumpConst <- function(j, map){
    if (is.null(j)){
        j <- (2.4/sqrt(length(map$coefficients)))^2
    }
    j
}

initRNG <- function(){
    if (!exists(".Random.seed")){ runif(1)  }
    .Random.seed # Retain and add to output
}
harrysouthworth/texmex documentation built on March 8, 2024, 7:50 p.m.