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

#### Documented in PWMsPWMs.data.framePWMs.listPWMs.matrixPWMs.numericPWMs.TLMomentsreturnPWMs

#' @title
#' Probability weighted moments
#' @description
#' Calculates probability weighted moments up to a specific order. Note that PWMs start with
#' order 0. Acceptable input types are numeric vectors, matrices, lists, and data.frames.
#'
#' @param x numeric vector or matrix, list, or data.frame of data OR an object of TLMoments.
#' @param formula if x is of type data.frame a formula has to be submitted.
#' @param max.order integer, maximal order of PWMs.
#' @param na.rm logical, indicates if NAs should be removed.
#'
#' @return numeric vector, matrix, list, or data.frame consisting of the PWMs and
#' with class \code{PWMs}.
#' The object contains the following attributes: \itemize{
#'  \item \code{order}: a integer vector with corresponding PWM orders
#'  \item \code{source}: a list with background information (used function, data, n, formula;
#'  mainly for internal purposes)
#' }
#' The attributes are hidden in the print-function for a clearer presentation.
#' @references Greenwood, J. A., Landwehr, J. M., Matalas, N. C., & Wallis, J. R. (1979). Probability weighted moments: definition and relation to parameters of several distributions expressable in inverse form. Water Resources Research, 15(5), 1049-1054.
#'
#' @examples
#' # Generating data sets:
#' xmat <- matrix(rnorm(100), nc = 4)
#' xvec <- xmat[, 3]
#' xlist <- lapply(1L:ncol(xmat), function(i) xmat[, i])
#' xdat <- data.frame(
#'  station = rep(letters[1:2], each = 50),
#'  season = rep(c("S", "W"), 50),
#'  hq = as.vector(xmat)
#' )
#'
#' # Calculating PWMs from data:
#' PWMs(xvec)
#' PWMs(xmat)
#' PWMs(xlist)
#' PWMs(xdat, formula = hq ~ station)
#' PWMs(xdat, formula = hq ~ season)
#' PWMs(xdat, formula = hq ~ .)
#' PWMs(xdat, formula = . ~ station + season)
#'
#' # Calculating PWMs from L-moments:
#' PWMs(TLMoments(xvec))
#' PWMs(TLMoments(xmat))
#' PWMs(TLMoments(xlist))
#' PWMs(TLMoments(xdat, hq ~ station))
#' PWMs(TLMoments(xdat, hq ~ season))
#' PWMs(TLMoments(xdat, hq ~ .))
#' PWMs(TLMoments(xdat, . ~ station + season))
#'
#' # In data.frame-mode invalid names are preceded by "."
#' xdat <- data.frame(
#'  beta0 = rep(letters[1:2], each = 50),
#'  beta1 = as.vector(xmat)
#' )
#' PWMs(xdat, formula = beta1 ~ beta0)
#'
#' @rdname PWMs
#' @export
PWMs <- function(x, ...) UseMethod("PWMs")

#' @title returnPWMs
#' @description Sets attributes to PWMs objects and returns them. This function is for internal use.
#' @param out -
#' @param order -
#' @param ... -
#' @return An object of class PWMs.
returnPWMs <- function(out, order, ...) {

class <- class(out)
args <- list(...)

# If no func attribute is set, set to
if (!exists("func", args)) args$func <- "PWMs" # If more than one func attributes are given, concatenate them if (sum(names(args) == "func") >= 2) { newfunc <- as.vector(unlist(args[names(args) == "func"])) args$func <- NULL
args$func <- newfunc } # Calculate n if not available and data exists. if (!exists("n", args) && exists("data", args)) {# && args$func == "PWMs") {
if (inherits(out, "numeric")) {
args$n <- sum(!is.na(args$data))
} else if (inherits(out, "matrix")) {
args$n <- apply(args$data, 2, function(y) sum(!is.na(y)))
} else if (inherits(out, "list")) {
args$n <- vapply(args$data, length, numeric(1))
} else if (inherits(out, "data.frame")) {
args$n <- aggregate(args$formula, args$data, length)[[getFormulaSides(args$formula)$lhs]] } } # Attributes of PWMs # order # source: func # data (if calculated) # n (if calculated) # formula (if data is data.frame) # lambdas (if coming from TLMoments) # trimmings (if coming from TLMoments) # class: "PWMs" attr(out, "order") <- order attr(out, "source") <- args class(out) <- c("PWMs", class) out } #' @rdname PWMs #' @method PWMs numeric #' @export PWMs.numeric <- function(x, max.order = 4L, na.rm = FALSE, ...) { out <- setNames( PWM(x, order = 0L:max.order, na.rm = na.rm), paste0("beta", 0:max.order) ) returnPWMs(out, 0L:max.order, data = x) } #' @rdname PWMs #' @method PWMs matrix #' @export PWMs.matrix <- function(x, max.order = 4L, na.rm = FALSE, ...) { out <- apply(x, 2, PWM, order = 0L:max.order, na.rm = na.rm) returnPWMs(out, 0L:max.order, data = x) } #' @rdname PWMs #' @method PWMs list #' @export PWMs.list <- function(x, max.order = 4L, na.rm = FALSE, ...) { out <- lapply(x, PWM, order = 0L:max.order, na.rm = na.rm) returnPWMs(out, 0L:max.order, data = x) } #' @rdname PWMs #' @method PWMs data.frame #' @export PWMs.data.frame <- function(x, formula, max.order = 4L, na.rm = FALSE, ...) { # Check for and repair invalid variables named beta[0-9] x <- correctNames(x, "beta[0-9]*", ".") formula <- correctNames(formula, "beta[0-9]*", ".") nam <- getFormulaSides(formula, names(x)) r <- aggregate(nam$new.formula, data = x, FUN = PWM, order = 0L:max.order, na.rm = na.rm)
out <- cbind(r[-length(r)], as.data.frame(r[[length(r)]]))

returnPWMs(out, 0L:max.order, data = x, formula = nam$new.formula) } #' @rdname PWMs #' @method PWMs TLMoments #' @export PWMs.TLMoments <- function(x, ...) { if (attr(x, "leftrim") != 0 | attr(x, "rightrim") != 0) stop("Transformation to PWMs only works for L-moments. ") if (any(diff(attr(x, "order"))) != 1) stop("Transformation to PWMs only runs for sequent L-moments. ") UseMethod("PWMs.TLMoments", x$lambdas)
}

#' @method PWMs.TLMoments numeric
#' @export
PWMs.TLMoments.numeric <- function(x, ...) {
max.order <- max(attr(x, "order"))
out <- as.numeric(solve(Z_C(max.order, 0, 0)) %*% x$lambdas) names(out) <- paste0("beta", 0:(max.order-1)) do.call( returnPWMs, c( list(out = out, order = 0L:(max.order-1)), func = "PWMs", lambdas = list(x$lambdas),
trimmings = list(c(attr(x, "leftrim"), attr(x, "rightrim"))),
tl.order = list(attr(x, "order")),
attr(x, "source")
)
)
}

#' @method PWMs.TLMoments matrix
#' @export
PWMs.TLMoments.matrix <- function(x, ...) {
max.order <- max(attr(x, "order"))
out <- solve(Z_C(max.order, 0, 0)) %*% x$lambdas rownames(out) <- paste0("beta", 0:(max.order-1)) do.call( returnPWMs, c( list(out = out, order = 0L:(max.order-1)), func = "PWMs", lambdas = list(x$lambdas),
trimmings = list(c(attr(x, "leftrim"), attr(x, "rightrim"))),
tl.order = list(attr(x, "order")),
attr(x, "source")
)
)
}

#' @method PWMs.TLMoments list
#' @export
PWMs.TLMoments.list <- function(x, ...) {
max.order <- max(attr(x, "order"))
A <- solve(Z_C(max.order, 0, 0))
out <- lapply(x$lambdas, function(x) setNames(as.numeric(A %*% x), paste0("beta", 0:(max.order-1)))) do.call( returnPWMs, c( list(out = out, order = 0L:(max.order-1)), func = "PWMs", lambdas = list(x$lambdas),
trimmings = list(c(attr(x, "leftrim"), attr(x, "rightrim"))),
tl.order = list(attr(x, "order")),
attr(x, "source")
)
)
}

#' @method PWMs.TLMoments data.frame
#' @export
PWMs.TLMoments.data.frame <- function(x, ...) {
max.order <- max(attr(x, "order"))

ls <- x$lambdas[, grep("L[0-9]*", names(x$lambdas))]
fac <- x$lambdas[, !grepl("L[0-9]*", names(x$lambdas)), drop = FALSE]

out <- as.data.frame(
t(solve(Z_C(max.order, 0, 0)) %*% t(as.matrix(ls)))
)
names(out) <- paste0("beta", 0:(max.order-1))
out <- cbind(fac, out)

do.call(
returnPWMs, c(
list(out = out, order = 0L:(max.order-1)),
func = "PWMs",
lambdas = list(x\$lambdas),
trimmings = list(c(attr(x, "leftrim"), attr(x, "rightrim"))),
tl.order = list(attr(x, "order")),
attr(x, "source")
)
)
}

#' @export
print.PWMs <- function(x, ...) {
if (inherits(x, "data.frame")) {
print.data.frame(x)
return(invisible(x))
}

tmp <- x
attributes(tmp) <- NULL

dim(tmp) <- dim(x)
names(tmp) <- names(x)
dimnames(tmp) <- dimnames(x)

print(tmp)
invisible(x)
}


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