# R/weighted_quantile.R In helske/diagis: Diagnostic Plot and Multivariate Summary Statistics of Weighted Samples from Importance Sampling

#### Documented in weighted_quantile

```#' Weighted quantiles
#'
#' Computes a weighted quantiles of a vector or matrix. Based on the formula in
#' Wikipedia (see the vignette) which is one of many ways to compute weighted
#' quantiles.
#'
#' @note Compared to some other R functions, here the weights are regarded as
#' probability weights, not frequency weights.
#'
#' @export
#' @param x A numeric vector or matrix. For matrix, the quantiles are computed
#' for each column.
#' @param w A numeric vector of non-negative weights. Will be automatically
#' normalised to sum to one.
#' @param probs A numeric vector of probabilities with values between 0 and 1.
#' @param na.rm If \code{TRUE}, \code{NA} and \code{NaN} values in \code{x}
#' (and corresponding weights in \code{w}) are omitted from the computation.
#' Default is \code{FALSE}. Additional missing values in \code{w} are not
#' allowed.
#' @return A weighted variance.
weighted_quantile <- function(x, w, probs = probs, na.rm) {
UseMethod("weighted_quantile", x)
}
#' @export
#' @method weighted_quantile ts
weighted_quantile.ts <- function(x, w, probs, na.rm = FALSE) {
weighted_quantile(unclass(x), w, probs, na.rm)
}
#' @export
#' @method weighted_quantile mcmc
weighted_quantile.mcmc <- function(x, w, probs, na.rm = FALSE) {
dimx <- dim(x)
if (is.null(dimx)) {
weighted_quantile.numeric(x, w, probs, na.rm)
} else {
if (length(dimx) == 2) {
weighted_quantile.matrix(x, w, probs, na.rm)
} else {
weighted_quantile.numeric(x, w, probs, na.rm)
}
}
}

#' @export
#' @method weighted_quantile matrix
weighted_quantile.matrix <- function(x, w, probs, na.rm = FALSE) {
if (nrow(x) != length(w))
stop("Length of 'w' is not equal to the number of rows in 'x'. ")
apply(x, 2, weighted_quantile.numeric, w = w, probs = probs, na.rm = na.rm)
}

#' @export
#' @method weighted_quantile numeric
weighted_quantile.numeric <- function(x, w, probs, na.rm = FALSE) {

if (!(typeof(w) %in% c("integer", "double"))) {
stop("Argument 'w' must be of type 'integer' or 'double'. ")
}
if (!is.numeric(w) || any(w < 0)) stop ("Weight vector 'w' must contain only non-negative values. ")
if (!any(w > 0)) {
stop("No positive weights in 'w'.")
}
if (any(!is.finite(probs)))
stop("Nonfinite values in 'probs'. ")
if (any(probs < 0 | probs > 1))
stop("'probs' outside [0, 1].")
if (length(na.rm) > 1 || !is.logical(na.rm))
stop("Argument 'na.rm' should be a logical of length one.")

if (any(nas <- is.na(x))) {
if (na.rm) {
x <- x[!nas]
w <- w[!nas]
} else {
stop("Missing values or NaN in 'x' not allowed as 'na.rm' is FALSE.")
}
}
if (anyNA(w))
stop("Missing values or NaNs in 'w' not allowed.")

if (any(dups <- duplicated(x))) {
w <- cumsum(tapply(w, x, sum))
w <- w / w[length(w)]
x <- sort(x[!dups])

} else {
ord <- order(x)
x <- x[ord]
w <- cumsum(w[ord])
w <- w / w[length(w)]
}
out <- numeric(length(probs))
for (i in seq_along(probs)) {
k <- which(w >= probs[i])[1]
if (k == 1) {
out[i] <- x[1]
} else {
out[i] <- x[k - 1] +
(x[k] - x[k - 1]) * (probs[i] - w[k - 1]) / (w[k] - w[k - 1])
}
}
names(out) <- paste0(formatC(probs * 100, format = "fg"), "%")
out
}
```
helske/diagis documentation built on Sept. 9, 2023, 8:15 a.m.