Merge pull request #351 from n-kall/pareto_k_tail_fix

Pareto k tail fix

R/pareto_smooth.R

@@ -279,7 +279,7 @@ pareto_smooth.default <- function(x,
   if (are_log_weights) {
     tail <- "right"
   }
-  
+
   tail <- match.arg(tail)
   S <- length(x)
 
@@ -328,9 +328,9 @@ pareto_smooth.default <- function(x,
     )
     right_k <- smoothed$k
 
-    k <- max(left_k, right_k)
+    k <- max(left_k, right_k, na.rm = TRUE)
     x <- smoothed$x
-    
+
   } else {
 
     smoothed <- .pareto_smooth_tail(
@@ -444,7 +444,7 @@ pareto_convergence_rate.rvar <- function(x, ...) {
     # shift log values for safe exponentiation
     x <- x - max(x)
   }
-  
+
   tail <- match.arg(tail)
 
   S <- length(x)
@@ -457,11 +457,21 @@ pareto_convergence_rate.rvar <- function(x, ...) {
   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)
+  }
+
   cutoff <- ord$x[min(tail_ids) - 1] # largest value smaller than tail values
-  
+
   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) {
@@ -617,7 +627,7 @@ pareto_k_diagmsg <- function(diags, are_weights = FALSE, ...) {
   msg <- NULL
 
   if (!are_weights) {
-  
+
     if (khat > 1) {
       msg <- paste0(msg, " Mean does not exist, making empirical mean estimate of the draws not applicable.")
     } else {

tests/testthat/test-pareto_smooth.R

@@ -6,6 +6,22 @@ test_that("pareto_khat returns expected reasonable values", {
 
 })
 
+
+test_that("pareto_khat handles constant tail correctly", {
+
+  # left tail is constant, so khat should be NA, but for "both" it
+  # should be the same as the right tail
+  x <- c(rep(-100, 10), sort(rnorm(100)))
+
+  expect_true(is.na(pareto_khat(x, tail = "left", ndraws_tail = 10)))
+  expect_equal(
+    pareto_khat(x, tail = "right", ndraws_tail = 10),
+    pareto_khat(x, tail = "both", ndraws_tail = 10)
+  )
+
+})
+
+
 test_that("pareto_khat handles tail argument", {
 
   # as tau is bounded (0, Inf) the left pareto k should be lower than
@@ -180,14 +196,27 @@ test_that("pareto_khat functions work with rvars with and without chains", {
 
 })
 
-test_that("pareto_smooth returns x with smoothed tail", {
-  tau <- extract_variable_matrix(example_draws(), "tau")
+test_that("pareto_smooth returns x with smoothed tail(s)", {
+  mu <- extract_variable_matrix(example_draws(), "mu")
+
+  mu_smoothed_right <- pareto_smooth(mu, ndraws_tail = 10, tail = "right", return_k = TRUE)$x
+
+  mu_smoothed_left <- pareto_smooth(mu, ndraws_tail = 10, tail = "left", return_k = TRUE)$x
+
+  mu_smoothed_both <- pareto_smooth(mu, ndraws_tail = 10, tail = "both", return_k = TRUE)$x
+
+  expect_equal(sort(mu)[1:390], sort(mu_smoothed_right)[1:390])
+  expect_equal(sort(mu_smoothed_both)[11:400], sort(mu_smoothed_right)[11:400])
 
-  tau_smoothed <- pareto_smooth(tau, ndraws_tail = 10, tail = "right", return_k = TRUE)$x
+  expect_equal(sort(mu)[11:400], sort(mu_smoothed_left)[11:400])
+  expect_equal(sort(mu_smoothed_both)[1:390], sort(mu_smoothed_left)[1:390])
 
-  expect_equal(sort(tau)[1:390], sort(tau_smoothed)[1:390])
+  expect_false(isTRUE(all.equal(sort(mu), sort(mu_smoothed_left))))
+  expect_false(isTRUE(all.equal(sort(mu), sort(mu_smoothed_right))))
+  expect_false(isTRUE(all.equal(sort(mu), sort(mu_smoothed_both))))
 
-  expect_false(isTRUE(all.equal(sort(tau), sort(tau_smoothed))))
+  expect_false(isTRUE(all.equal(sort(mu_smoothed_both), sort(mu_smoothed_left))))
+  expect_false(isTRUE(all.equal(sort(mu_smoothed_both), sort(mu_smoothed_right))))
 
 })