Add argument to force continuous or discrete Mode

R/point_interval.R

@@ -102,6 +102,11 @@ globalVariables(c("y", "ymin", "ymax"))
 #' and unbounded data.
 #' @param n For [hdi()] and [Mode()], the number of points to use to estimate highest-density
 #' intervals or modes.
+#' @param type For [Mode()], the type of estimator to use. Can be `"discrete"` to compute
+#' the discrete mode (the single most frequently occurring value) or `"continuous"`
+#' to compute the maximum a posteriori (MAP) estimate (the maximum value of the
+#' the estimated density). The default `"auto"` uses `"discrete"` if `x` is
+#' an integer or non-numeric (see [`rlang::is_integerish()`]) and `"continuous"` otherwise.
 #' @param weights For [Mode()], an optional vector, which (if not `NULL`) is of the same length
 #' as `x` and provides weights for each element of `x`.
 #' @return A data frame containing point summaries and intervals, with at least one column corresponding
@@ -522,20 +527,29 @@ hdi_.distribution = function(x, .width = .95, na.rm = FALSE, ..., density = dens
 #' @rdname point_interval
 #' @importFrom rlang is_integerish
 #' @importFrom stats density
-Mode = function(x, na.rm = FALSE, ...) {
+Mode = function(x, na.rm = FALSE, type = c("auto", "discrete", "continuous"), ...) {
   UseMethod("Mode")
 }
 #' @export
 #' @rdname point_interval
-Mode.default = function(x, na.rm = FALSE, ..., density = density_bounded(trim = TRUE), n = 2001, weights = NULL) {
+Mode.default = function(x, na.rm = FALSE, type = c("auto", "discrete", "continuous"), ..., density = density_bounded(trim = TRUE), n = 2001, weights = NULL) {
+  type = match.arg(type)
   if (na.rm) {
     x = x[!is.na(x)]
   } else if (anyNA(x)) {
     return(NA_real_)
   }
   density = match_function(density, "density_")
 
-  if (is_integerish(x)) {
+  if (type == "auto") {
+    if (is_integerish(x) || !is.numeric(x)) {
+      type = "discrete"
+    } else {
+      type = "continuous"
+    }
+  }
+
+  if (type == "discrete") {
     if (is.null(weights)) {
       # for the discrete case, based on https://stackoverflow.com/a/8189441
       ux = unique(x)
@@ -553,20 +567,20 @@ Mode.default = function(x, na.rm = FALSE, ..., density = density_bounded(trim =
 }
 #' @export
 #' @rdname point_interval
-Mode.rvar = function(x, na.rm = FALSE, ...) {
+Mode.rvar = function(x, na.rm = FALSE, type = c("auto", "discrete", "continuous"), ...) {
   draws = posterior::draws_of(x)
   .dim = dim(draws)
-  apply(draws, seq_along(.dim)[-1], Mode, na.rm = na.rm, weights = weights(x))
+  apply(draws, seq_along(.dim)[-1], Mode, na.rm = na.rm, type = type, weights = weights(x))
 }
 #' @importFrom stats optim
 #' @export
 #' @rdname point_interval
-Mode.distribution = function(x, na.rm = FALSE, ...) {
+Mode.distribution = function(x, na.rm = FALSE, type = c("auto", "discrete", "continuous"), ...) {
   find_mode = function(x) {
     if (anyNA(x)) {
       NA_real_
     } else if (distr_is_sample(x)) {
-      Mode(distr_get_sample(x), na.rm = na.rm, weights = distr_get_sample_weights(x))
+      Mode(distr_get_sample(x), na.rm = na.rm, type = type, weights = distr_get_sample_weights(x))
     } else if (distr_is_constant(x)) {
       quantile(x, 0.5)[[1]]
     } else if (distr_is_discrete(x)) {