Nothing
#' 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
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.