R/rMvtPar.R

Defines functions simulPareto

Documented in simulPareto

#' Simulate Pareto random vectors
#'
#' \code{simulPareto} provides \code{n} replicates of the multivariate Pareto distribution
#' associated to log-Gaussian random function with semi-variogram \code{vario}.
#'
#' 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} replicates conditioned
#' that \code{mean(x) > 1} on the unit Frechet scale.
#'
#' @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 used 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()}.
#' @references Dombry, C., Engelke, S. and M. Oesting. Exact simulation of max-stable processes. Biometrika, 103(2), 303-317.
#' @return List of \code{n} random vectors drawn from a multivariate Pareto distribution with semi-variogram \code{vario}.
#' @examples
#' #Define 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 <- simulPareto(100, loc, vario)
#' @export


simulPareto <- 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 valide for the provided locations.')
  })

  k <- sample(1:dim, n, replace = TRUE)
  d <- nrow(gamma)

  cov.mat <- (outer(gamma[-1,1],gamma[-1,1], "+") - (gamma[-1,-1]))
  chol.mat <- chol(cov.mat)

  proc <- matrix(0,d,n)
  proc[-1,] <- t(chol.mat)%*%matrix(rnorm((d- 1)*n), ncol=n)

  # # for(i in 1:n){
  # #   buffer <- exp(proc[,i] - proc[k[i],i] - gamma[,k[i]])
  # #   proc[,i] <-evd::rgpd(1, loc=1, scale=1, shape=1) * buffer / mean(buffer)
  # # }
  #
  # proc[proc == 0] = .Machine$double.xmin

  sims <- lapply(1:n, function(i){
    buffer <- exp(proc[,i] - proc[k[i],i] - gamma[,k[i]])
    proc[,i] <-evd::rgpd(1, loc=1, scale=1, shape=1) * buffer / mean(buffer)
    proc[proc == 0] = .Machine$double.xmin
    proc[,i]})
  return(sims)
}

Try the mvPot package in your browser

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

mvPot documentation built on Oct. 14, 2023, 1:06 a.m.