diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 98cd6543..976296c6 100644 --- a/R/pareto_smooth.R +++ b/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 { diff --git a/tests/testthat/test-pareto_smooth.R b/tests/testthat/test-pareto_smooth.R index 89635600..7a6c9393 100644 --- a/tests/testthat/test-pareto_smooth.R +++ b/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)))) })