Merge pull request #442 from stan-dev/fix-pareto_smooth

fix pareto_smooth in case of constant tail + cleanup

Author: avehtari

R/pareto_smooth.R

@@ -361,6 +361,10 @@ pareto_smooth.default <- function(x,
     x <- smoothed$x
   }
 
+  if (is.na(k)) {
+    return(pareto_diags_na(x, return_k, extra_diags))
+  }
+
   diags_list <- list(khat = k)
 
   if (extra_diags) {
@@ -470,38 +474,42 @@ pareto_convergence_rate.rvar <- function(x, ...) {
 #' 
 #' @export
 ps_tail <- function(x,
-                                ndraws_tail,
-                                smooth_draws = TRUE,
-                                tail = c("right", "left"),
-                                are_log_weights = FALSE,
-                                ...
-                                ) {
+                    ndraws_tail,
+                    smooth_draws = TRUE,
+                    tail = c("right", "left"),
+                    are_log_weights = FALSE,
+                    ...
+                    ) {
+
+  if (ndraws_tail < 5) {
+    warning_no_call(
+      "Can't fit generalized Pareto distribution ",
+      "because ndraws_tail is less than 5."
+    )
+    return(list(x = x, k = NA))
+  }
 
   if (are_log_weights) {
     # shift log values for safe exponentiation
     x <- x - max(x)
   }
 
   tail <- match.arg(tail)
-
-  ndraws <- length(x)
-  tail_ids <- seq(ndraws - ndraws_tail + 1, ndraws)
-
   if (tail == "left") {
     x <- -x
   }
 
+  ndraws <- length(x)
+  tail_ids <- seq(ndraws - ndraws_tail + 1, ndraws)
+
   ord <- sort.int(x, index.return = TRUE)
   draws_tail <- ord$x[tail_ids]
 
   if (is_constant(draws_tail)) {
-
     if (tail == "left") {
       x <- -x
     }
-
-    out <- list(x = x, k = NA)
-    return(out)
+    return(list(x = x, k = NA))
   }
 
   cutoff <- ord$x[min(tail_ids) - 1] # largest value smaller than tail values
@@ -511,43 +519,22 @@ ps_tail <- function(x,
   }
 
   max_tail <- max(draws_tail)
-  min_tail <- min(draws_tail)
-
-  if (ndraws_tail >= 5) {
-    ord <- sort.int(x, index.return = TRUE)
-    if (abs(max_tail - min_tail) < .Machine$double.eps / 100) {
-      warning_no_call(
-        "Can't fit generalized Pareto distribution ",
-        "because all tail values are the same."
-      )
-      smoothed <- NULL
-      k <- NA
-    } else {
-      # save time not sorting since x already sorted
-      if (are_log_weights) {
-        draws_tail <- exp(draws_tail)
-        cutoff <- exp(cutoff)
-      }
-      fit <- gpdfit(draws_tail - cutoff, sort_x = FALSE, ...)
-      k <- fit$k
-      sigma <- fit$sigma
-      if (is.finite(k) && smooth_draws) {
-        p <- (seq_len(ndraws_tail) - 0.5) / ndraws_tail
-        smoothed <- qgeneralized_pareto(p = p, mu = cutoff, k = k, sigma = sigma)
-        if (are_log_weights) {
-          smoothed <- log(smoothed)
-        }
-      } else {
-        smoothed <- NULL
-      }
+
+  if (are_log_weights) {
+    draws_tail <- exp(draws_tail)
+    cutoff <- exp(cutoff)
+  }
+  fit <- gpdfit(draws_tail - cutoff, sort_x = FALSE, ...)
+  k <- fit$k
+  sigma <- fit$sigma
+  if (is.finite(k) && smooth_draws) {
+    p <- (seq_len(ndraws_tail) - 0.5) / ndraws_tail
+    smoothed <- qgeneralized_pareto(p = p, mu = cutoff, k = k, sigma = sigma)
+    if (are_log_weights) {
+      smoothed <- log(smoothed)
     }
   } else {
-    warning_no_call(
-      "Can't fit generalized Pareto distribution ",
-      "because ndraws_tail is less than 5."
-    )
     smoothed <- NULL
-    k <- NA
   }
 
   # truncate at max of raw draws
@@ -598,7 +585,10 @@ ps_tail <- function(x,
 #' @return minimum sample size
 #' @export
 ps_min_ss <- function(k, ...) {
-  if (k < 1) {
+  if (is.na(k)) {
+    return(NA)
+  }
+  if (isTRUE(k < 1)) {
     out <- 10^(1 / (1 - max(0, k)))
   } else {
     out <- Inf

tests/testthat/test-pareto_smooth.R

@@ -13,6 +13,19 @@ test_that("pareto_khat handles constant tail correctly", {
 })
 
 
+test_that("pareto_smooth handles non-constant x with constant tail", {
+
+  # x is non-constant but the right tail is constant
+  x <- c(seq(-3, 3, length.out = 90), rep(5, 10))
+
+  expect_no_error(
+    ps <- pareto_smooth(x, tail = "right", ndraws_tail = 10, return_k = TRUE, verbose = FALSE)
+  )
+  expect_true(is.na(ps$diagnostics$khat))
+  expect_equal(ps$x, x)
+
+})
+
 test_that("pareto_khat handles tail argument", {
 
   # as tau is bounded (0, Inf) the left pareto k should be lower than
@@ -235,3 +248,37 @@ test_that("pareto_smooth works for log_weights", {
   expect_true(ps$diagnostics$khat > 0.7)
 
 })
+
+test_that("check ps_tail behavior for ndraws_tail less than 5", {
+  w <- c(1:25, 1e3, 1e3, 1e3)
+  lw <- log(w)
+
+  # prints correct warning
+  expect_warning(
+    tail <- ps_tail(lw, ndraws_tail = 4, tail = "right"),
+    "Can't fit generalized Pareto distribution because ndraws_tail is less than 5."
+  )
+  # output has expected shape and k = NA
+  expect_equal(names(tail), c("x", "k"))
+  expect_true(is.na(tail$k))
+})
+
+test_that("check ps_tail behavior for constant draws_tail", {
+  x <- log(replicate(10, 0.3))
+
+  tail <- ps_tail(x, ndraws_tail = 10, tail = "left")
+
+  # output has expected return values
+  expect_equal(x, tail$x)
+  expect_true(is.na(tail$k))
+})
+
+
+test_that("check ps_min_ss behavior special cases", {
+  # k = NA
+  expect_true(is.na(ps_min_ss(NA)))
+  # k > 1
+  expect_true(is.infinite(ps_min_ss(2)))
+  # k < 1
+  expect_equal(ps_min_ss(0.5), 10^(1 / (1 - max(0, 0.5))))
+})