diff --git a/NEWS.md b/NEWS.md index d43f875f..bcca5839 100644 --- a/NEWS.md +++ b/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. diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 70261de0..1b283eb5 100644 --- a/R/pareto_smooth.R +++ b/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) diff --git a/R/pit.R b/R/pit.R index d3789c47..a8b579ce 100644 --- a/R/pit.R +++ b/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 { diff --git a/tests/testthat/test-pareto_pit.R b/tests/testthat/test-pareto_pit.R index a6b4716f..d552302c 100644 --- a/tests/testthat/test-pareto_pit.R +++ b/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", { diff --git a/tests/testthat/test-pareto_smooth.R b/tests/testthat/test-pareto_smooth.R index 928aae5d..a79e0bd8 100644 --- a/tests/testthat/test-pareto_smooth.R +++ b/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)