add generalised mean + tests

NAMESPACE

@@ -57,6 +57,7 @@ export(Screen)
 export(Treat)
 export(a_amean)
 export(a_copeland)
+export(a_genmean)
 export(a_gmean)
 export(a_hmean)
 export(approx_df)

R/aggregate.R

@@ -547,7 +547,7 @@ Aggregate <- function(x, ...){
 #'
 #' The vector of weights `w` is relative since the formula is:
 #'
-#' \deqn{ y = 1(\sum w) \sum wx }
+#' \deqn{ y = \frac{1}{\sum w_i} \sum w_i x_i }
 #'
 #' If `x` contains `NA`s, these `x` values and the corresponding `w` values are removed before applying the
 #' formula above.
@@ -689,6 +689,69 @@ a_hmean <- function(x, w = NULL){
 }
 
 
+#' Weighted generalised mean
+#'
+#' Weighted generalised mean of a vector. `NA` are skipped by default.
+#'
+#' The generalised mean is as follows:
+#'
+#' \deqn{ y = \left( \frac{1}{\sum w_i} \sum w_i x_i^p  \right)^{1/p} }
+#'
+#' where `p` is a coefficient specified in the function argument here. Note that:
+#'
+#' - For negative `p`, all `x` values must be positive
+#' - Setting `p = 0` will result in an error due to the negative exponent. This case
+#' is equivalent to the geometric mean in the limit, so use [a_gmean()] instead.
+#'
+#' @param x A numeric vector of positive values.
+#' @param w A vector of weights, which should have length equal to `length(x)`. Weights are relative
+#' and will be re-scaled to sum to 1. If `w` is not specified, defaults to equal weights.
+#' @param p Coefficient - see details.
+#'
+#' @examples
+#' # a vector of values
+#' x <- 1:10
+#' # a vector of weights
+#' w <- runif(10)
+#' # weighted harmonic mean
+#' a_genmean(x,w)
+#'
+#' @return Weighted harmonic mean, as a numeric value.
+#'
+#' @export
+a_genmean <- function(x, w = NULL, p){
+
+  if(is.null(w)){
+    # default equal weights
+    w <- rep(1,length(x))
+    message("No weights specified, using equal weights.")
+  }
+
+  if(p==0){
+    stop("Setting p = 0 results in an infinite exponent. In the limit, this case is equal to the geometric mean: use a_gmean() instead.")
+  }
+
+  if(any(!is.na(x))){
+
+    if(any(x <= 0, na.rm = TRUE) && p < 0){
+      stop("Zero or negative values found when applying generalised mean with negative p: cannot be calculated.")
+    }
+
+    # have to set any weights to NA to correspond to NAs in x
+    w[is.na(x)] <- NA
+
+    gm <- (sum(w*x^p)/sum(w))^(1/p)
+
+  } else {
+    gm <- NA
+  }
+
+  gm
+
+}
+
+
+
 #' Outranking matrix
 #'
 #' Constructs an outranking matrix based on a data frame of indicator data and corresponding weights.

tests/testthat/test-aggregate.R

@@ -125,19 +125,25 @@ test_that("agg_functions", {
 
   # geometric
   x <- c(1, 2, 3)
-  y1 <- a_gmean(x)
+  yg <- a_gmean(x)
 
-  expect_equal(y1, (1*2*3)^(1/3))
+  expect_equal(yg, (1*2*3)^(1/3))
   expect_error(a_gmean(c(1, 2, -1)))
   expect_error(a_gmean(c(1, 2, 0)))
 
-  y1 <- a_gmean(x, c(1,1,2))
-  expect_equal(y1, (1*2*3^2)^(1/4))
+  yg2 <- a_gmean(x, c(1,1,2))
+  expect_equal(yg2, (1*2*3^2)^(1/4))
 
   # harmonic
   x <- c(1, 2, 3)
-  y1 <- a_hmean(x, c(2, 1, 1))
-  expect_equal(y1, 4/(2/1 + 1/2 + 1/3))
+  yh <- a_hmean(x, c(2, 1, 1))
+  expect_equal(yh, 4/(2/1 + 1/2 + 1/3))
+
+  # generalised
+  ygen <- a_genmean(x, c(2,1,1), p = -1)
+  expect_equal(ygen, yh)
+  ygen <- a_genmean(x, p = 1)
+  expect_equal(ygen, 2) # simple arithmetic av.
 
 })