R/rBrownResnick.R

In r-fndv/mvPot: Multivariate Peaks-over-Threshold Modelling for Spatial Extreme Events

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