Nothing
#' Simulation of Brown--Resnick random vectors
#'
#' \code{simulBrownResnick} provides \code{n} replicates of a Brown--Resnick max-stable process with semi-variogram \code{vario}
#' at locations \code{loc}.
#'
#' The algorithm used here is based on the spectral representation of the Brown--Resnick
#' model as described in Dombry et al. (2015). It provides \code{n} exact simulations
#' on the unit Frechet scale and requires, in average, for each max-stable vector, the simulation of d Pareto processes,
#' where d is the number of locations.
#'
#' @param n Number of replicates desired.
#' @param loc Matrix of coordinates as given by \code{expand.grid()}.
#' @param vario Semi-variogram function.
#' @param nCores Number of cores needed for the computation
#' @param cl Cluster instance as created by \code{makeCluster} of the \code{parallel} package. Make sure
#' the random number generator has been properly initialized with
#' \code{clusterSetRNGStream()}.
#' @return List of \code{n} random vectors drawn from a max-stable Brown--Resnick process
#' with semi-variogram \code{vario} at location \code{loc}.
#' @references Dombry, C., Engelke, S. and M. Oesting. Exact simulation of max-stable processes. Biometrika, 103(2), 303-317.
#' @examples
#' #Define semi-variogram function
#' vario <- function(h){
#' 1 / 2 * norm(h,type = "2")^1.5
#' }
#'
#' #Define locations
#' loc <- expand.grid(1:4, 1:4)
#'
#' #Simulate data
#' obs <- simulBrownResnick(10, loc, vario)
#' @export
simulBrownResnick <- function(n, loc, vario, nCores = 1, cl = NULL){
if(!inherits(loc, "data.frame")) {
stop('`loc` must be the data frame of coordinates as generated by `expand.grid()`')
}
dim <- nrow(loc)
if(!is.numeric(nCores) || nCores < 1) {
stop('`nCores` must a positive number of cores to use for parallel computing.')
}
if(nCores > 1 && !inherits(cl, "cluster")) {
stop('For parallel computation, `cl` must an cluster created by `makeCluster` of the package `parallel`.')
}
gamma <- tryCatch({
dists <- lapply(1:ncol(loc), function(i) {
outer(loc[,i],loc[,i], "-")
})
computeVarMat <- sapply(1:length(dists[[1]]), function(i){
h <- rep(0,ncol(loc))
for(j in 1:ncol(loc)){
h[j] = dists[[j]][i]
}
vario(h)
})
matrix(computeVarMat, dim, dim)
}, warning = function(war) {
war
}, error = function(err) {
stop('The semi-variogram provided is not valid for the provided locations.')
})
simFun <- function(i) {
print(i)
poisson <- stats::rexp(1)
brownResnick <- 1/poisson * paretoBR(1, dim, gamma)
for(i in 2:dim){
poisson <- stats::rexp(1)
while(1 / poisson > brownResnick[i]){
Y <- 1/poisson * paretoBR(i, dim, gamma)
if(all(brownResnick[1:(i - 1)] > Y[1: (i - 1)])) {
brownResnick <- pmax(brownResnick, Y)
}
poisson <- poisson + stats::rexp(1)
}
}
brownResnick
}
if(nCores > 1){
sims <- parallel::parLapply(cl, 1:n, simFun)
} else {
sims <- lapply(1:n, simFun)
}
return(sims)
}
paretoBR <- function(k, dim, gamma){
#Generate covariance and mean for the Gaussian field
mean <- - gamma[-k,k]
cov <- (outer(gamma[-k,k],gamma[-k,k], "+") - (gamma[-k,-k]))
dim <- nrow(gamma)
paretoProcess <- rep(1,dim)
paretoProcess[-k] <- exp(MASS::mvrnorm(n = 1, mu = mean, Sigma = cov))
return(paretoProcess)
}
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.