Nothing
#' Simulation from Pareto processes using composition sampling
#'
#' The algorithm performs forward sampling by simulating first from a
#' mixture, then sample angles conditional on them being less than (max) or greater than (min) one.
#' The resulting sample from the angular distribution is then multiplied by
#' Pareto variates with tail index \code{shape}.
#'
#' For the moment, only exchangeable models and models based n elliptical processes are handled.
#'
#'
#' The parametrization of the Brown--Resnick is in terms of the matrix \code{Lambda}, which
#' is formed by evaluating the semivariogram \eqn{\gamma} at sites \eqn{s_i, s_j}, meaning that
#' \eqn{\Lambda_{i,j} = \gamma(s_i, s_j)/2}.
#'
#' @author Leo Belzile
#' @param n sample size.
#' @param Lambda parameter matrix for the Brown--Resnick model. See \bold{Details}.
#' @param Sigma correlation matrix if \code{model = 'xstud'}, otherwise
#' the covariance matrix formed from the stationary Brown-Resnick process.
#' @param df degrees of freedom for extremal Student process.
#' @param model string indicating the model family.
#' @param shape tail index of the Pareto variates (reciprocal shape parameter). Must be strictly positive.
#' @param param parameter value for the logistic or negative logistic model
#' @param d dimension of the multivariate model, only needed for logistic or negative logistic models
#' @param risk string indicating the risk functional. Only \code{max} and \code{min} are currently supported.
#' @param ... additional parameters, currently ignored
#' @details The argument \code{Sigma} is ignored for the Brown-Resnick model
#' if \code{Lambda} is provided by the user.
#' @export
#' @seealso \code{\link{rparp}} for general simulation of Pareto processes based on an accept-reject algorithm.
#' @return an \code{n} by \code{d} matrix of samples, where \code{d = ncol(Sigma)}, with \code{attributes} \code{mixt.weights}.
#' @examples
#' \dontrun{
#' #Brown-Resnick, Wadsworth and Tawn (2014) parametrization
#' D <- 20L
#' coord <- cbind(runif(D), runif(D))
#' semivario <- function(d, alpha = 1.5, lambda = 1){0.5 * (d/lambda)^alpha}
#' Lambda <- semivario(as.matrix(dist(coord))) / 2
#' rparpcs(n = 10, Lambda = Lambda, model = 'br', shape = 0.1)
#' #Extremal Student
#' Sigma <- stats::rWishart(n = 1, df = 20, Sigma = diag(10))[,,1]
#' rparpcs(n = 10, Sigma = cov2cor(Sigma), df = 3, model = 'xstud')
#' }
rparpcs <- function(n,
model = c("log", "neglog", "br", "xstud"),
risk = c("max", "min"),
param = NULL,
d,
Lambda = NULL,
Sigma = NULL,
df = NULL,
shape = 1,
...) {
args <- list(...)
if(!is.null(args$riskf)){
riskf <- args$riskf
} else{
riskf <- risk
}
riskf <- match.arg(arg = riskf,
choices = c("min", "max"),
several.ok = TRUE)[1]
model <- match.arg(model)
stopifnot(shape > 0)
if(model %in% c("br", "xstud")){
if (!requireNamespace("TruncatedNormal", quietly = TRUE)) {
stop(
"Package \"TruncatedNormal\" must be installed to use this function.",
call. = FALSE
)
}
if(!utils::packageVersion("TruncatedNormal") > "1.2"){
stop(
"Please update package \"TruncatedNormal\" to a more recent version.",
call. = FALSE
)
}
}
if(model %in% c("log", "neglog")){
if(missing(d)){
stop("Invalid dimension for the model.")
}
if (is.null(param) || param < 0 || d < 1) {
stop("Invalid parameter value")
}
if (length(param) != 1) {
warning("Only first entry of param vector considered")
param <- param[1]
}
if (model == "log") {
if (param < 1) {
param <- 1/param
}
}
weights <- rep(1/d, d)
out <- matrix(1, nrow = d, ncol = n)
# Exchangeable model, balanced mixture
if(model == "log"){
scale <- 1/gamma(1-1/param)
F0 <- scale*rgamma(n,
shape = 1-1/param)^(-1/param)
if(riskf == "min"){
Flow <- rep(
x = mev::pgev(
q = F0,
loc = scale,
scale = scale/param,
shape = 1/param),
each = d-1)
Fj <- mev::qgev(
p = Flow + (1-Flow) * runif(n * (d-1)),
loc = scale,
scale = scale/param,
shape = 1/param)
} else if(riskf == "max"){
Fj <- mev::qgev(
p = runif(n * (d-1)) *
rep(mev::pgev(F0,
loc = scale,
scale = scale/param,
shape = 1/param), d-1),
loc = scale,
scale = scale/param,
shape = 1/param)
}
} else if(model == "neglog"){
scale <- 1/gamma(1+1/param)
F0 <- scale*rgamma(n,
shape = 1+1/param)^(-1/param)
if(riskf == "min"){
Flow <- rep(
x = pweibull(
q = F0,
scale = scale,
shape = param),
each = d-1)
Fj <- qweibull(
p = Flow + (1-Flow) * runif(n * (d-1)),
scale = scale,
shape = param)
} else if(riskf == "max"){
Fj <- qweibull(
p = runif(n * (d-1)) *
rep(pweibull(
q = F0,
scale = scale,
shape = param), d-1),
scale = scale,
shape = param)
}
}
# pos <- (sample.int(n = d, size = n, replace = TRUE)-1) * n + 1:n
pos <- sample.int(n = d, size = n, replace = TRUE) +
d*(0:(n-1))
out[-pos] <- Fj/rep(F0, each = d - 1)
samp <- t(out) * runif(n)^(-shape)
} else if(model %in% c("br", "xstud")){
if (model == "xstud") {
if (is.null(df)) {
stop("Invalid degree of freedom argument")
}
stopifnot(!is.null(Sigma))
D <- nrow(Sigma)
if (!isTRUE(all.equal(as.vector(diag(Sigma)), rep(1, D)))) {
warning("Input \"Sigma\" must be a correlation matrix")
Sigma <- cov2cor(Sigma)
}
weights <- .weightsXstud(z = rep(1, D), Sigma = Sigma, df = df, riskf = riskf)
} else if (model == "br") {
if (!is.null(Lambda)) {
D <- nrow(Lambda)
weights <- .weightsBR(z = rep(1, D), Lambda = Lambda, riskf = riskf)
} else if (!is.null(Sigma)) {
D <- nrow(Sigma)
T1 <- cbind(rep(-1, D - 1), diag(D - 1))
weights <- .weightsBR_WT(z = rep(1, D), Sigma = Sigma, riskf = riskf)
}
}
weights <- weights/sum(weights)
id <- sample.int(n = D, size = n, replace = TRUE, prob = weights)
tabu <- table(id)
bins <- as.integer(names(tabu))
tabu <- as.vector(tabu)
# Matrix to store samples
ang <- matrix(1, nrow = D, ncol = n)
accu <- 0L
for (i in 1:length(tabu)) {
j <- bins[i]
if (riskf == "max") {
if (model == "xstud") {
ang[-j, (accu + 1L):(accu + tabu[i])] <- pmax(TruncatedNormal::mvrandt(n = tabu[i], l = rep(-Inf, D - 1), u = rep(1,
D - 1), df = df + 1, Sig = (Sigma[-j, -j] - Sigma[-j, j, drop = FALSE] %*% Sigma[j, -j, drop = FALSE])/(df + 1)),
0)^df
} else if (model == "br") {
if (!is.null(Lambda)) {
ang[-j, (accu + 1L):(accu + tabu[i])] <- exp(TruncatedNormal::mvrandn(n = tabu[i], l = rep(-Inf, D - 1), u = 2 *
Lambda[j, -j], Sig = 2 * (outer(Lambda[-j, j], Lambda[j, -j], FUN = "+") - Lambda[-j, -j])) - 2 * Lambda[j, -j])
} else if (!is.null(Sigma)) {
me <- diag(Sigma)[-j]/2 + Sigma[j, j]/2 - Sigma[j, -j]
Ti <- T1
Ti[, c(1L, j)] <- T1[, c(j, 1L)]
ang[-j, (accu + 1L):(accu + tabu[i])] <- exp(TruncatedNormal::mvrandn(n = tabu[i], l = rep(-Inf, D - 1), u = -me,
Sig = Ti %*% Sigma %*% t(Ti)) + me)
}
}
} else if (riskf == "min") {
if (model == "xstud") {
ang[-j, (accu + 1L):(accu + tabu[i])] <- pmax(TruncatedNormal::mvrandt(n = tabu[i], u = rep(Inf, D - 1), l = rep(1,
D - 1), df = df + 1, Sig = (Sigma[-j, -j] - Sigma[-j, j, drop = FALSE] %*% Sigma[j, -j, drop = FALSE])/(df + 1)),
0)^df
} else if (model == "br") {
if (!is.null(Lambda)) {
ang[-j, (accu + 1L):(accu + tabu[i])] <- exp(TruncatedNormal::mvrandn(n = tabu[i], u = rep(Inf, D - 1), Sig = 2 *
(outer(Lambda[-j, j], Lambda[j, -j], FUN = "+") - Lambda[-j, -j]), l = 2 * Lambda[j, -j]) - 2 * Lambda[j, -j])
} else if (!is.null(Sigma)) {
me <- diag(Sigma)[-j]/2 + Sigma[j, j]/2 - Sigma[j, -j]
Ti <- T1
Ti[, c(1L, j)] <- T1[, c(j, 1L)]
ang[-j, (accu + 1L):(accu + tabu[i])] <- exp(TruncatedNormal::mvrandn(n = tabu[i], u = rep(Inf, D - 1), l = -me,
Sig = Ti %*% Sigma %*% t(Ti)) + me)
}
}
}
accu <- accu + tabu[i]
}
samp <- runif(n)^(-shape) * t(ang[, sample.int(n, n, replace = FALSE)])
}
attr(samp, "mixt.weights") <- weights
return(samp)
}
#' Simulation of generalized Huesler-Reiss Pareto vectors via composition sampling
#'
#' Sample from the generalized Pareto process associated to Huesler-Reiss spectral profiles.
#' For the Huesler-Reiss Pareto vectors, the matrix \code{Sigma} is utilized to build \eqn{Q} viz.
#' \deqn{Q = \Sigma^{-1} - \frac{\Sigma^{-1}\mathbf{1}_d\mathbf{1}_d^\top\Sigma^{-1}}{\mathbf{1}_d^\top\Sigma^{-1}\mathbf{1}_d}.}
#' The location vector \code{m} and \code{Sigma} are the parameters of the underlying log-Gaussian process.
#'
#' @param n sample size
#' @param u vector of marginal location parameters (must be strictly positive)
#' @param alpha vector of shape parameters (must be strictly positive).
#' @param Sigma covariance matrix of process, used to define \eqn{Q}. See \bold{Details}.
#' @param m location vector of Gaussian distribution.
#' @return \code{n} by d matrix of observations
#' @references Ho, Z. W. O and C. Dombry (2019), Simple models for multivariate regular variations and the
#' Huesler-Reiss Pareto distribution, Journal of Multivariate Analysis (\bold{173}), p. 525-550, \doi{10.1016/j.jmva.2019.04.008}
#' @export
#' @examples
#' D <- 20L
#' coord <- cbind(runif(D), runif(D))
#' di <- as.matrix(dist(rbind(c(0, ncol(coord)), coord)))
#' semivario <- function(d, alpha = 1.5, lambda = 1){(d/lambda)^alpha}
#' Vmat <- semivario(di)
#' Sigma <- outer(Vmat[-1, 1], Vmat[1, -1], '+') - Vmat[-1, -1]
#' m <- Vmat[-1,1]
#' \dontrun{
#' samp <- rparpcshr(n = 100, u = c(rep(1, 10), rep(2, 10)),
#' alpha = seq(0.1, 1, length = 20), Sigma = Sigma, m = m)
#' }
rparpcshr <- function(n, u, alpha, Sigma, m) {
if (!requireNamespace("TruncatedNormal", quietly = TRUE)) {
stop(
"Package \"TruncatedNormal\" must be installed to use this function.",
call. = FALSE
)
}
D <- ncol(Sigma)
Sigmainv <- solve(Sigma)
q <- Sigmainv %*% rep(1, D)
Q <- (Sigmainv - q %*% t(q)/sum(q))
l <- c(Sigmainv %*% (((m %*% q - 1)/sum(q))[1] * rep(1, D) - m))
if (any(u <= 0)) {
stop("Threshold must be strictly positive")
}
if (any(alpha <= 0)) {
stop("Tail indices must all be strictly positive")
}
alpha <- rep(alpha, length = D)
u <- rep(u, length = D)
Lst <- l - Q %*% diag(alpha) %*% log(u) #prop 4.2 b/c u \neq 1
weights <- .weightsHR(z = rep(1, D), L = Lst, Q = Q)
weights <- weights/sum(weights)
id <- sample.int(n = D, size = n, replace = TRUE, prob = weights)
tabu <- table(id)
bins <- as.integer(names(tabu))
tabu <- as.vector(tabu)
# Matrix to store samples
ang <- matrix(1, nrow = D, ncol = n)
accu <- 0L
for (i in 1:length(tabu)) {
j <- bins[i]
Qjinv <- solve(Q[-j, -j])
ang[-j, (accu + 1L):(accu + tabu[i])] <- exp(TruncatedNormal::mvrandn(n = tabu[i], l = rep(-Inf, D - 1), u = c(-Qjinv %*%
Lst[-j, ]), Sig = Qjinv) + c(Qjinv %*% Lst[-j, ]))
accu <- accu + tabu[i]
}
R <- mev::rgp(n = n, loc = 1, scale = 1, shape = 1)
samp <- R * t(ang[, sample.int(n, n, replace = FALSE)])
for (j in 1:ncol(samp)) {
samp[, j] <- u[j] * (samp[, j]^(alpha[j])) #coro 4.3
}
return(samp)
}
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.