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 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,
                   verbose, trace, theCall, ...){
    if (class(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)

    ######################## Run the Metropolis algorithm...
    res <- texmexMetropolis(res, log.lik, proposals, verbose, trace)

    # 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, seed=seed)

    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, log.lik, proposals, verbose, trace){
    last.cost <- log.lik(x[1,])
    if (!is.finite(last.cost)) {
      stop("Start is infeasible.")
    }

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

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

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

####### 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
}

Try the texmex package in your browser

Any scripts or data that you put into this service are public.

texmex documentation built on May 2, 2019, 5:41 a.m.