# R/est_quancov.R In TLMoments: Calculate TL-Moments and Convert Them to Distribution Parameters

#### Documented in est_quancovest_quancov.matrixest_quancov.numericest_quancov.quantiles

```#' @title
#' Estimate the covariance matrix of quantile estimations
#' @description
#' Internal function. Use \link{est_cov}. Description not done yet.
#'
#' @param x numeric vector or matrix containing data.
#' @param distr character of length 1 giving the distribution if parametric assumption should be used.
#' @param p quantile levels from which the covariance should be calculated.
#' @param leftrim,rightrim lower and upper trimming parameter used for parameter calculation, have to be non-negative integers.
#' @param np.cov boolean, if TRUE no parametric assumptions are used to calculate the covariance matrix (default FALSE).
#' @param reg.weights numeric vector of weights for regionalized TLMoments.
#' @param set.n hypothetical data length n if theoretical values are given.
#'
#' @return numeric matrix
#'
#' @examples
#' ### Numeric vectors
#' x <- rgev(500, shape = .2)
#'
#' quantiles(parameters(TLMoments(x), "gev"), c(.9, .95, .99))
#' est_quancov(x, "gev", c(.9, .95, .99), 0, 0)
#' #cov(t(replicate(5000,
#' #  quantiles(parameters(TLMoments(rgev(500, shape = .2)), "gev"), c(.9, .95, .99))
#' #)))
#'
#' quantiles(parameters(TLMoments(x, rightrim = 1), "gev"), c(.9, .95, .99))
#' est_quancov(x, "gev", c(.9, .95, .99), 0, 1)
#' #cov(t(replicate(5000,
#' #  quantiles(
#' #    parameters(TLMoments(rgev(500, shape = .2), rightrim = 1), "gev"),
#' #    c(.9, .95, .99)
#' #  )
#' #)))
#'
#' ### Numeric matrices
#' x <- matrix(rgev(600, shape = .2), nc = 3)
#'
#' quantiles(parameters(TLMoments(x), "gev"), c(.9, .95, .99))
#' est_quancov(x, "gev", c(.9, .95, .99), 0, 0)
#'
#' est_quancov(x, "gev", .9, 0, 0)
#' #cov(t(replicate(5000,
#' # quantiles(
#' #   parameters(TLMoments(matrix(rgev(600, shape = .2), nc = 3)),
#' #  "gev"), .9)
#' #  )
#' #))
#'
#' ### quantiles object
#' q <- quantiles(as.parameters(loc = 3, scale = 2, shape = .4, distr = "gev"), c(.9, .99))
#' est_quancov(q)
#' est_quancov(q, leftrim = 0, rightrim = 0)
#' est_quancov(q, leftrim = 0, rightrim = 0, set.n = 10)
#'
#' @rdname est_quancov
#' @export
est_quancov <- function(x,
distr = "",
p = NULL,
leftrim = 0L,
rightrim = 0L,
...) {
if (inherits(x, "quantiles")) {
distr <- attr(x, "distribution")
p <- attr(x, "p")
if (!any(is.null(attr(x, "source")\$trimming))) {
leftrim <- attr(x, "source")\$trimming[1]
rightrim <- attr(x, "source")\$trimming[2]
}
}

# TODO: General error catches
if (is.na(leftrim) | is.na(rightrim))
stop("leftrim and/or right must be given. ")

if (!are.integer.like(leftrim, rightrim))
stop("leftrim and rightrim have to be integer types!")

if (!(distr %in% c("gev", "gum", "gpd", "ln3")))
stop("distr has to be \"gev\", \"gum\", \"gpd\", or \"ln3\". ")

if (!(distr %in% c("gev")))
stop("only GEV for now")

if (!is.numeric(p))
stop("Argument p must be numeric vector. ")

UseMethod("est_quancov")
}

#' @rdname est_quancov
#' @method est_quancov numeric
#' @export
est_quancov.numeric <- function(x,
distr,
p,
leftrim = 0L,
rightrim = 0L,
np.cov = FALSE,
...) {

if (np.cov) {
paramcov <- est_paramcov(x, distr, leftrim = leftrim, rightrim = rightrim, np.cov = TRUE)
} else {
paramcov <- est_paramcov(x, distr, leftrim = leftrim, rightrim = rightrim)
}
param <- parameters(TLMoments(x, leftrim = leftrim, rightrim = rightrim, na.rm = TRUE), distr)
A <- CovParamtoQuan(distr, param, p)
r <- A %*% paramcov %*% t(A)
rownames(r) <- colnames(r) <- buildNames("q", p)

attr(r, "distribution") <- distr
attr(r, "trimmings") <- c(leftrim, rightrim)
attr(r, "n") <- length(x)
r
}

#' @rdname est_quancov
#' @method est_quancov matrix
#' @export
est_quancov.matrix <- function(x,
distr,
p,
leftrim = 0L,
rightrim = 0L,
np.cov = FALSE,
reg.weights = NULL,
...) {

if (np.cov) {
paramcov <- est_paramcov(x, distr, leftrim, rightrim, np.cov = TRUE, reg.weights = reg.weights)
} else {
paramcov <- est_paramcov(x, distr, leftrim, rightrim, reg.weights = reg.weights)
}

if (is.null(reg.weights)) {
param <- parameters(TLMoments(x, leftrim = leftrim, rightrim = rightrim, na.rm = TRUE), distr)
A <- lapply(1:ncol(param), function(i) CovParamtoQuan(distr, param[, i], p))
A <- blockdiag_list(A)
r <- A %*% paramcov %*% t(A)
rownames(r) <- colnames(r) <- buildNames("q", p, 1:ncol(x))
} else {
param <- parameters(regionalize(TLMoments(x, leftrim = leftrim, rightrim = rightrim, na.rm = TRUE), reg.weights), distr)
A <- CovParamtoQuan(distr, param, p)
r <- A %*% paramcov %*% t(A)
rownames(r) <- colnames(r) <- buildNames("q", p)
}

attr(r, "distribution") <- distr
attr(r, "trimmings") <- c(leftrim, rightrim)
attr(r, "n") <- apply(x, 2, function(y) sum(!is.na(y)))
r
}

#' @rdname est_quancov
#' @method est_quancov quantiles
#' @export
est_quancov.quantiles <- function(x,
distr = attr(x, "distribution"),
p = attr(x, "p"),
leftrim = attr(x, "source")\$trimmings[1],
rightrim = attr(x, "source")\$trimmings[2],
set.n = NA,
...) {
if (!inherits(x, "numeric"))
stop("x must be a numeric quantiles-object. ")
if (!any(attr(x, "source")\$func %in% c("as.PWMs", "as.TLMoments", "as.parameters")))
stop("est_quancov.quantiles only for theoretical values. ")

if (are.integer.like(leftrim) && is.null(rightrim)) {
rightrim <- 0
}
if (are.integer.like(rightrim) && is.null(leftrim)) {
leftrim <- 0
}
if (any(is.null(c(leftrim, rightrim)))) {
#warning("No trimmings found, set to leftrim = 0, rightrim = 0. ")
leftrim <- rightrim <- 0
}

if (is.na(set.n) | !is.numeric(set.n)) {
#warning("Missing or invalid set.n argument. Giving results for n = 100. ")
n <- 100
} else { n <- set.n }

param <- as.parameters(attr(x, "source")\$parameters, distr = distr)
paramcov <- est_paramcov(param, leftrim = leftrim, rightrim = rightrim, set.n = n)
A <- CovParamtoQuan(distr, param, p)
r <- A %*% paramcov %*% t(A)
rownames(r) <- colnames(r) <- buildNames("q", p)

attr(r, "distribution") <- distr
attr(r, "trimmings") <- c(leftrim, rightrim)
attr(r, "n") <- n
r
}
```

## Try the TLMoments package in your browser

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

TLMoments documentation built on March 27, 2022, 5:07 p.m.