Remove the offset at a tied Pareto tail cutoff

NEWS.md

@@ -11,9 +11,11 @@ selected.
 
 * Fix truncation operation for `uniformity_test()` (#462)
 * mid-p value adjusment for discrete `PRIT` test method  
-* Add additional unit tests for `uniformity_test()`
+* Add additional unit tests for `uniformity_test()` 
+* Remove the offset at a tied Pareto tail cutoff in `ps_tail()` and `pareto_pit()` (#468)
 * Use a tolerance relative to the magnitude of `x` in `is_constant()` (#469)
 
+
 # posterior 1.7.1
 
 * Minor release to accommodate the JOSS paper.

R/pareto_smooth.R

@@ -516,11 +516,9 @@ ps_tail <- function(x,
     return(list(x = x, k = NA))
   }
 
-  cutoff <- ord$x[min(tail_ids) - 1] # largest value smaller than tail values
-  if (cutoff == ord$x[min(tail_ids)]) {
-    # cutoff is not smaller than the tail values
-    cutoff <- cutoff - .Machine$double.eps
-  }
+  # largest value not in the tail; tail values equal to it give zero
+  # excesses, which are in the support of the generalized Pareto distribution
+  cutoff <- ord$x[min(tail_ids) - 1]
 
   max_tail <- max(draws_tail)
 

R/pit.R

@@ -320,9 +320,6 @@ pareto_pit.draws_matrix <- function(x, y, weights = NULL, log = FALSE,
     right_tail <- sorted[tail_ids]
     if (!right_tail_empty && !is_constant(right_tail)) {
       right_cutoff <- sorted[min(tail_ids) - 1]
-      if (right_cutoff == right_tail[1]) {
-        right_cutoff <- right_cutoff - .Machine$double.eps
-      }
       right_tail_wt <- if (!is.null(log_wt_sorted)) {
         exp(log_wt_sorted[tail_ids])
       } else {
@@ -354,9 +351,6 @@ pareto_pit.draws_matrix <- function(x, y, weights = NULL, log = FALSE,
       left_tail <- left_sorted[tail_ids]
       if (!left_tail_empty && !is_constant(left_tail)) {
         left_cutoff <- left_sorted[min(tail_ids) - 1]
-        if (left_cutoff == left_tail[1]) {
-          left_cutoff <- left_cutoff - .Machine$double.eps
-        }
         left_tail_wt <- if (!is.null(log_wt_left_sorted)) {
           exp(log_wt_left_sorted[tail_ids])
         } else {

tests/testthat/test-pareto_pit.R

@@ -188,6 +188,33 @@ test_that("pareto_pit uses randomization for discrete observations in bulk", {
   expect_true(all(results >= 0 & results <= 1))
 })
 
+test_that("pareto_pit uses the raw PIT at a tail cutoff that ties with tail draws", {
+  # the tail cutoff (0) ties with the lowest 10 tail draws
+  draws <- c(rep(0, 350), exp(-seq(0, 5, length.out = 50)))
+  min_pit <- 1 / length(draws) / 1e4
+
+  # right tail: the result does not depend on the location of the draws
+  x0 <- as_draws_matrix(matrix(draws, ncol = 1))
+  x5 <- as_draws_matrix(matrix(draws + 5, ncol = 1))
+  set.seed(1)
+  pit0 <- pareto_pit(x0, 0)
+  set.seed(1)
+  pit5 <- pareto_pit(x5, 5)
+  expect_equal(unname(pit0), unname(pit5))
+  # an observation below all draws is not in the right tail
+  expect_equal(unname(pareto_pit(x0, -1e-300)), min_pit)
+
+  # left tail: same checks on the mirrored draws
+  xl0 <- as_draws_matrix(matrix(-draws, ncol = 1))
+  xl5 <- as_draws_matrix(matrix(-draws - 5, ncol = 1))
+  set.seed(1)
+  pitl0 <- pareto_pit(xl0, 0)
+  set.seed(1)
+  pitl5 <- pareto_pit(xl5, -5)
+  expect_equal(unname(pitl0), unname(pitl5))
+  expect_equal(unname(pareto_pit(xl0, 1e-300)), 1 - min_pit)
+})
+
 # ---- ndraws_tail argument ----
 
 test_that("pareto_pit respects ndraws_tail argument", {

tests/testthat/test-pareto_smooth.R

@@ -304,6 +304,32 @@ test_that("ps_tail handles -Inf log weights below the cutoff", {
   expect_true(all(is.infinite(tail$x[1:80]) & tail$x[1:80] < 0))
 })
 
+test_that("ps_tail does not smooth tail values below the cutoff", {
+  # many weights underflow to 0, so the cutoff ties with the lowest tail
+  # values; the smoothed tail must not go below the cutoff (here below 0)
+  w <- c(rep(0, 350), exp(-seq(0, 60, length.out = 50)))
+  tail <- ps_tail(w, ndraws_tail = 60, tail = "right")
+  expect_true(is.finite(tail$k))
+  expect_true(all(tail$x >= 0))
+  # the order of the draws does not change
+  expect_true(all(diff(tail$x[order(w)]) >= 0))
+
+  # same for the left tail
+  tail <- ps_tail(-w, ndraws_tail = 60, tail = "left")
+  expect_true(is.finite(tail$k))
+  expect_true(all(tail$x <= 0))
+})
+
+test_that("pareto_smooth does not return negative weights for ties at 0", {
+  w <- c(rep(0, 350), exp(-seq(0, 60, length.out = 50)))
+  ws <- pareto_smooth(w, tail = "right", r_eff = 1, return_k = FALSE, verbose = FALSE)
+  expect_true(all(ws >= 0))
+  x <- as_draws_df(data.frame(theta = seq_along(ws)))
+  expect_no_error(
+    resample_draws(x, ndraws = 4, weights = ws, method = "simple_no_replace")
+  )
+})
+
 
 test_that("check ps_tail behavior for ndraws_tail less than 5", {
   w <- c(1:25, 1e3, 1e3, 1e3)