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