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

#### Documented in est_pwmcovest_pwmcov.matrixest_pwmcov.numeric

```#' @title
#' Estimate the covariance matrix of PWM estimations
#' @description
#' Internal function. Use \link{est_cov}. Description not done yet.
#'
#' @param x numeric vector or matrix of data.
#' @param order numeric vector giving the orders that are returned.
#' @param distr character of length 1 which indicates a distribution if a
#' parametric assumption should be used.
#' @param distr.trim integer vector of length 2 indicating the trimming used
#' to calculate parameters if a parametric assumption is used (i.e. \code{distr} is set).
#'
#' @return numeric matrix
#'
#' @examples
#' ### Numeric vectors
#' x <- rgev(500, shape = .2)
#' est_pwmcov(x)
#' est_pwmcov(x, distr = "gev")
#'
#' ### Numeric matrices
#' x <- matrix(rgev(600, shape = .2), nc = 3)
#' est_pwmcov(x, order = 0:2)
#' est_pwmcov(x, order = 0:2, distr = "gev")
#'
#' @rdname est_pwmcov
#' @export
est_pwmcov <- function(x, order = 0:3, distr = NULL, distr.trim = c(0, 0)) {
if (inherits(x, c("PWMs", "TLMoments", "parameters", "quantiles")))
stop("est_pwmcov is only for data vectors or matrices. ")

if (is.double(order)) order <- as.integer(order)
if (!is.integer(order)) stop("order has to be integer type!")

if (!is.null(distr) && distr != "gev")
stop("distr has to be NULL (nonparametric) or \"gev\"")
if (!is.null(distr) && (!are.integer.like(distr.trim[1], distr.trim[2]) | length(distr.trim) != 2))
stop("If distr is not NULL, distr.trim has to be a integer-like vector of length 2. ")

UseMethod("est_pwmcov")
}

#' @rdname est_pwmcov
#' @method est_pwmcov numeric
#' @export
est_pwmcov.numeric <- function(x, order = 0:3, distr = NULL, distr.trim = c(0, 0)) {

if (is.null(distr)) {
# nonparametric estimation
r <- stats::cov(pseudo(x, order))
} else {
# parametric estimation
# scale and shape are estimated using L-moments.
p <- parameters(TLMoments(x, leftrim = distr.trim[1], rightrim = distr.trim[2], na.rm = TRUE), distr)

if (p["shape"] >= .5) {
warning("estimated shape parameter is too high. ")
r <- matrix(Inf, nrow = length(order), ncol = length(order))
} else {
r <- parametricPWMCov(distr, order, scale = p["scale"], shape = p["shape"])
}
}

rownames(r) <- colnames(r) <- buildNames("beta", order)
out <- r / sum(!is.na(x))
attr(out, "distribution") <- distr
out
}

#' @rdname est_pwmcov
#' @method est_pwmcov matrix
#' @export
est_pwmcov.matrix <- function(x, order = 0:3, distr = NULL, distr.trim = c(0, 0)) {

J <- ncol(x) # Anzahl Pegel
# Berechnung von Pseudobeobachtungen fuer jede Dimension
pseudos <- lapply(1L:J, function(i) pseudo(x[, i], order))
# Berechnung der Sigma-Matrix (ohne Korrekturvorfaktor)
sigma <- stats::cov(do.call("cbind", pseudos), use = "pairwise.complete.obs")
rownames(sigma) <- colnames(sigma) <- buildNames("beta", order, 1:J)

if (!is.null(distr)) {
# parametric estimation, overwrites non-diagonal block matrices
# scale and shape are estimated using L-moments.
xl <- lapply(1:ncol(x), function(i) x[, i])
p <- parameters(TLMoments(xl, leftrim = distr.trim[1], rightrim = distr.trim[2], na.rm = TRUE), distr)
r <- lapply(p, function(p) {
parametricPWMCov(distr, order, scale = p["scale"], shape = p["shape"])
})
sigma <- blockdiag_list(r, sigma)
}

# Calculation correction factor
rs <- apply(x, 2, function(y) sum(!is.na(y))/length(y))
vorf <- do.call(rbind, lapply(1:J, function(i) {
do.call(cbind, lapply(1:J, function(j) {
matrix(min(c(rs[i], rs[j]))/(rs[i] * rs[j]), nrow = length(order), ncol = length(order))
}))
}))

# Ausgeben
out <- vorf*sigma / max(apply(x, 2, function(y) sum(!is.na(y))))
attr(out, "distribution") <- distr
out
}
```

## 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.