R Under development (unstable) (2023-12-09 r85665 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library("matrixStats") > > x <- 1:5 > y <- weightedMedian(x) > y <- weightedMedian(x, w = c(NA, Inf, NA, Inf, NA), na.rm = TRUE) > print(y) [1] 3 > > y <- weightedMedian(x, w = c(NA, Inf, NA, Inf, NA), na.rm = FALSE) > print(y) [1] NA > stopifnot(is.na(y)) > > x <- 1:10 > n <- length(x) > > y1 <- median(x) # 5.5 > y2 <- weightedMedian(x) # 5.5 > stopifnot(all.equal(y1, y2)) > > > w <- rep(1, times = n) > y1 <- weightedMedian(x, w) # 5.5 (default) > y2a <- weightedMedian(x, ties = "weighted") # 5.5 (default) > y2b <- weightedMedian(x, ties = "min") # 5 > y2c <- weightedMedian(x, ties = "max") # 6 > stopifnot(all.equal(y2a, y1)) > > y3 <- weightedMedian(x, w) # 5.5 (default) > > > # Pull the median towards zero > w[1] <- 5 > y1 <- weightedMedian(x, w) # 3.5 > y <- c(rep(0, times = w[1]), x[-1]) # Only possible for integer weights > y2 <- median(y) # 3.5 > stopifnot(all.equal(y1, y2)) > > # Put even more weight on the zero > w[1] <- 8.5 > y <- weightedMedian(x, w) # 2 > > # All weight on the first value > w[1] <- Inf > y <- weightedMedian(x, w) # 1 > > # All weight on the last value > w[1] <- 1 > w[n] <- Inf > y <- weightedMedian(x, w) # 10 > > # All weights set to zero > w <- rep(0, times = n) > y <- weightedMedian(x, w) # NA > > x <- 1:4 > w <- rep(1, times = 4) > for (mode in c("integer", "double")) { + storage.mode(x) <- mode + for (ties in c("weighted", "mean", "min", "max")) { + cat(sprintf("ties = %s\n", ties)) + y <- weightedMedian(x, w, ties = ties) + } + } ties = weighted ties = mean ties = min ties = max ties = weighted ties = mean ties = min ties = max > > set.seed(0x42) > > y <- weightedMedian(x = double(0L)) > print(y) [1] NA > stopifnot(length(y) == 1L) > stopifnot(is.na(y)) > > y <- weightedMedian(x = x[1]) > print(y) [1] 1 > stopifnot(length(y) == 1L) > stopifnot(all.equal(y, x[1])) > > > n <- 1e3 > x <- runif(n) > w <- runif(n, min = 0, max = 1) > for (mode in c("integer", "double")) { + storage.mode(x) <- mode + for (ties in c("weighted", "mean", "min", "max")) { + y <- weightedMedian(x, w, ties = ties) + cat(sprintf("mode = %s, ties = %s, result = %g\n", mode, ties, y)) + } + } mode = integer, ties = weighted, result = 0 mode = integer, ties = mean, result = 0 mode = integer, ties = min, result = 0 mode = integer, ties = max, result = 0 mode = double, ties = weighted, result = 0 mode = double, ties = mean, result = 0 mode = double, ties = min, result = 0 mode = double, ties = max, result = 0 > > > # A large vector > n <- 1e5 > x <- runif(n) > w <- runif(n, min = 0, max = 1) > y <- weightedMedian(x, w) > > y <- weightedMedian(x, w, ties = "min") > > > # Single Number > xs <- c(1, NA_integer_) > ws <- c(1, NA_integer_) > for (x in xs) { + for (w in ws) { + y <- weightedMedian(x = x, w = w) + if (is.na(w)) z <- NA_real_ + else z <- x[1] + stopifnot(all.equal(y, z)) + } + } > > ## Logical > x1 <- c(TRUE, FALSE, TRUE) > w0 <- c(0, 0, 0) > stopifnot(!is.finite(weightedMedian(x1, w0)), + !is.infinite(weightedMedian(x1, w0))) > > w1 <- c(1, 1, 1) > stopifnot(weightedMedian(x1, w1) == 1) > > w2 <- c(1, 2, 3) > stopifnot(weightedMedian(x1, w2) == 1) > > ### NA > stopifnot(is.na(weightedMedian(c(TRUE, FALSE, NA), + c(1, 2, 3))), + all.equal(weightedMedian(c(TRUE, FALSE, NA), + c(1, 2, 3), + na.rm = TRUE), + weightedMedian(c(TRUE, FALSE), + c(1, 2)))) > ### Identical to as.integer() > x <- rcauchy(100) > w <- abs(rcauchy(100)) > stopifnot(all.equal(weightedMedian(x > 0, w), + weightedMedian(as.integer(x > 0), w))) > > > > proc.time() user system elapsed 0.20 0.03 0.21