From f834fde2306bca9b8856c15974f55ba862c3e818 Mon Sep 17 00:00:00 2001 From: VisruthSK <67435125+VisruthSK@users.noreply.github.com> Date: Thu, 12 Jun 2025 12:08:26 -0700 Subject: [PATCH 01/11] Rely on posterior for pareto smooth tails --- R/psis.R | 57 +++++++------------------------------- tests/testthat/test_psis.R | 11 -------- 2 files changed, 10 insertions(+), 58 deletions(-) diff --git a/R/psis.R b/R/psis.R index 1b321d70..4f32effd 100644 --- a/R/psis.R +++ b/R/psis.R @@ -206,62 +206,25 @@ do_psis_i <- function(log_ratios_i, tail_len_i, ...) { S <- length(log_ratios_i) # shift log ratios for safer exponentation lw_i <- log_ratios_i - max(log_ratios_i) - khat <- Inf - if (enough_tail_samples(tail_len_i)) { - ord <- sort.int(lw_i, index.return = TRUE) - tail_ids <- seq(S - tail_len_i + 1, S) - lw_tail <- ord$x[tail_ids] - if (abs(max(lw_tail) - min(lw_tail)) < .Machine$double.eps/100) { - warning( - "Can't fit generalized Pareto distribution ", - "because all tail values are the same.", - call. = FALSE - ) - } else { - cutoff <- ord$x[min(tail_ids) - 1] # largest value smaller than tail values - smoothed <- psis_smooth_tail(lw_tail, cutoff) - khat <- smoothed$k - lw_i[ord$ix[tail_ids]] <- smoothed$tail - } - } + smoothed <- posterior::pareto_smooth_tail( + x = lw_i, + ndraws_tail = tail_len_i, + tail = "right", + are_log_weights = TRUE + ) + + lw_i <- smoothed$x + khat <- smoothed$k # truncate at max of raw wts (i.e., 0 since max has been subtracted) lw_i[lw_i > 0] <- 0 # shift log weights back so that the smallest log weights remain unchanged lw_i <- lw_i + max(log_ratios_i) - list(log_weights = lw_i, pareto_k = khat) -} - -#' PSIS tail smoothing for a single vector -#' -#' @noRd -#' @param x Vector of tail elements already sorted in ascending order. -#' @return A named list containing: -#' * `tail`: vector same size as `x` containing the logs of the -#' order statistics of the generalized pareto distribution. -#' * `k`: scalar shape parameter estimate. -#' -psis_smooth_tail <- function(x, cutoff) { - len <- length(x) - exp_cutoff <- exp(cutoff) - - # save time not sorting since x already sorted - fit <- gpdfit(exp(x) - exp_cutoff, sort_x = FALSE) - k <- fit$k - sigma <- fit$sigma - if (is.finite(k)) { - p <- (seq_len(len) - 0.5) / len - qq <- qgpd(p, k, sigma) + exp_cutoff - tail <- log(qq) - } else { - tail <- x - } - list(tail = tail, k = k) + list(log_weights = lw_i, pareto_k = if (is.na(khat)) Inf else khat) } - #' Calculate tail lengths to use for fitting the GPD #' #' The number of weights (i.e., tail length) used to fit the generalized Pareto diff --git a/tests/testthat/test_psis.R b/tests/testthat/test_psis.R index ef4b14cd..cc6d11cc 100644 --- a/tests/testthat/test_psis.R +++ b/tests/testthat/test_psis.R @@ -150,14 +150,3 @@ test_that("do_psis_i throws warning if all tail values the same", { val <- expect_warning(do_psis_i(xx, tail_len_i = 6), "all tail values are the same") expect_equal(val$pareto_k, Inf) }) - -test_that("psis_smooth_tail returns original tail values if k is infinite", { - # skip on M1 Mac until we figure out why this test fails only on M1 Mac - skip_if(Sys.info()[["sysname"]] == "Darwin" && R.version$arch == "aarch64") - - xx <- c(1,2,3,4,4,4,4,4,4,4,4) - val <- suppressWarnings(psis_smooth_tail(xx, 3)) - expect_equal(val$tail, xx) - expect_equal(val$k, Inf) -}) - From 0a8ea9d2b8c7e93b85aebfcc3418c5b212e7fd87 Mon Sep 17 00:00:00 2001 From: VisruthSK <67435125+VisruthSK@users.noreply.github.com> Date: Thu, 12 Jun 2025 12:27:12 -0700 Subject: [PATCH 02/11] Added warning if tail values are constant to match original functionality --- R/psis.R | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/R/psis.R b/R/psis.R index 4f32effd..4d3b631e 100644 --- a/R/psis.R +++ b/R/psis.R @@ -207,6 +207,14 @@ do_psis_i <- function(log_ratios_i, tail_len_i, ...) { # shift log ratios for safer exponentation lw_i <- log_ratios_i - max(log_ratios_i) + if (length(unique(tail(log_ratios_i, -tail_len_i))) == 1) { + warning( + "Can't fit generalized Pareto distribution ", + "because all tail values are the same.", + call. = FALSE + ) + } + smoothed <- posterior::pareto_smooth_tail( x = lw_i, ndraws_tail = tail_len_i, From e7db5b56298c7319764307fd8dd4e50539aab018 Mon Sep 17 00:00:00 2001 From: VisruthSK <67435125+VisruthSK@users.noreply.github.com> Date: Fri, 13 Jun 2025 09:50:50 -0700 Subject: [PATCH 03/11] Bump version of posterior --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 5fdb524d..7021cb49 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -36,7 +36,7 @@ Imports: checkmate, matrixStats (>= 0.52), parallel, - posterior (>= 1.5.0), + posterior (>= 1.6.2), stats Suggests: bayesplot (>= 1.7.0), From 698e484c35ff55e6094e1a44d9716538e5d031eb Mon Sep 17 00:00:00 2001 From: VisruthSK <67435125+VisruthSK@users.noreply.github.com> Date: Fri, 13 Jun 2025 10:15:08 -0700 Subject: [PATCH 04/11] Reverted version --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 7021cb49..d4a128e4 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -36,7 +36,7 @@ Imports: checkmate, matrixStats (>= 0.52), parallel, - posterior (>= 1.6.2), + posterior (>= 1.6.1), stats Suggests: bayesplot (>= 1.7.0), From 2c1ed4c14d7abc2fc62015e2ea5811592821fcaf Mon Sep 17 00:00:00 2001 From: Visruth Srimath Kandali Date: Tue, 7 Apr 2026 10:11:32 -0700 Subject: [PATCH 05/11] Update DESCRIPTION Bumped posterior dependency --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 4881488a..6456d907 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -37,7 +37,7 @@ Imports: checkmate, matrixStats (>= 0.52), parallel, - posterior (>= 1.6.1), + posterior (>= 1.7.0), stats Suggests: bayesplot (>= 1.7.0), From a8f0a18b4a8d8ccc9049463630bb72d0f472670f Mon Sep 17 00:00:00 2001 From: VisruthSK Date: Tue, 7 Apr 2026 10:24:28 -0700 Subject: [PATCH 06/11] Fixed function name --- R/psis.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/psis.R b/R/psis.R index ee183593..e30601b4 100644 --- a/R/psis.R +++ b/R/psis.R @@ -221,7 +221,7 @@ do_psis_i <- function(log_ratios_i, tail_len_i, ...) { ) } - smoothed <- posterior::pareto_smooth_tail( + smoothed <- posterior::ps_tail( x = lw_i, ndraws_tail = tail_len_i, tail = "right", From 2c7613426505be4df1f0b6d37890bb259b2b1a15 Mon Sep 17 00:00:00 2001 From: VisruthSK Date: Tue, 7 Apr 2026 10:24:35 -0700 Subject: [PATCH 07/11] Tentatively updating snapshots --- tests/testthat/_snaps/psis.md | 64 +++++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) diff --git a/tests/testthat/_snaps/psis.md b/tests/testthat/_snaps/psis.md index 9f2deaf3..2ce0adba 100644 --- a/tests/testthat/_snaps/psis.md +++ b/tests/testthat/_snaps/psis.md @@ -4801,6 +4801,70 @@ Warning: Not enough tail samples to fit the generalized Pareto distribution in some or all columns of matrix of log importance ratios. Skipping the following columns: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, ... [22 more not printed]. Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: + Can't fit generalized Pareto distribution because ndraws_tail is less than 5. + Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details. Output Computed from 10 by 32 log-weights matrix. From 99e0f6f07273a48b1dddd1cc1377b731a3a12205 Mon Sep 17 00:00:00 2001 From: Visruth Srimath Kandali Date: Tue, 7 Apr 2026 13:06:50 -0700 Subject: [PATCH 08/11] Update R/psis.R Co-authored-by: Jonah Gabry --- R/psis.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/psis.R b/R/psis.R index e30601b4..6ee30ce8 100644 --- a/R/psis.R +++ b/R/psis.R @@ -213,7 +213,7 @@ do_psis_i <- function(log_ratios_i, tail_len_i, ...) { # shift log ratios for safer exponentation lw_i <- log_ratios_i - max(log_ratios_i) - if (length(unique(tail(log_ratios_i, -tail_len_i))) == 1) { + if (length(unique(utils::tail(log_ratios_i, -tail_len_i))) == 1) { warning( "Can't fit generalized Pareto distribution ", "because all tail values are the same.", From 02f44664e652f2c71d39927d6b30b78d511395f8 Mon Sep 17 00:00:00 2001 From: VisruthSK Date: Tue, 7 Apr 2026 15:09:15 -0700 Subject: [PATCH 09/11] Supressed warnings --- R/psis.R | 4 +-- tests/testthat/_snaps/psis.md | 64 ----------------------------------- 2 files changed, 2 insertions(+), 66 deletions(-) diff --git a/R/psis.R b/R/psis.R index 6ee30ce8..b084ef7c 100644 --- a/R/psis.R +++ b/R/psis.R @@ -221,12 +221,12 @@ do_psis_i <- function(log_ratios_i, tail_len_i, ...) { ) } - smoothed <- posterior::ps_tail( + smoothed <- suppressWarnings(posterior::ps_tail( x = lw_i, ndraws_tail = tail_len_i, tail = "right", are_log_weights = TRUE - ) + )) lw_i <- smoothed$x khat <- smoothed$k diff --git a/tests/testthat/_snaps/psis.md b/tests/testthat/_snaps/psis.md index 2ce0adba..9f2deaf3 100644 --- a/tests/testthat/_snaps/psis.md +++ b/tests/testthat/_snaps/psis.md @@ -4801,70 +4801,6 @@ Warning: Not enough tail samples to fit the generalized Pareto distribution in some or all columns of matrix of log importance ratios. Skipping the following columns: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, ... [22 more not printed]. Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: - Can't fit generalized Pareto distribution because ndraws_tail is less than 5. - Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details. Output Computed from 10 by 32 log-weights matrix. From 8850e1410f2745d89758319bc0bc91758521e02f Mon Sep 17 00:00:00 2001 From: VisruthSK Date: Tue, 7 Apr 2026 21:36:22 -0700 Subject: [PATCH 10/11] Sort tail values --- R/psis.R | 2 +- tests/testthat/test_psis.R | 7 +++++++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/R/psis.R b/R/psis.R index b084ef7c..639f462f 100644 --- a/R/psis.R +++ b/R/psis.R @@ -213,7 +213,7 @@ do_psis_i <- function(log_ratios_i, tail_len_i, ...) { # shift log ratios for safer exponentation lw_i <- log_ratios_i - max(log_ratios_i) - if (length(unique(utils::tail(log_ratios_i, -tail_len_i))) == 1) { + if (length(unique(utils::tail(sort(log_ratios_i), tail_len_i))) == 1) { warning( "Can't fit generalized Pareto distribution ", "because all tail values are the same.", diff --git a/tests/testthat/test_psis.R b/tests/testthat/test_psis.R index 9ae108a5..f67bd05d 100644 --- a/tests/testthat/test_psis.R +++ b/tests/testthat/test_psis.R @@ -150,4 +150,11 @@ test_that("do_psis_i throws warning if all tail values the same", { "all tail values are the same" ) expect_equal(val$pareto_k, Inf) + + xx <- c(4, 1, 4, 4, 4, 4, 4, 2, 3, 4, 4) + expect_warning( + val <- do_psis_i(xx, tail_len_i = 6), + "all tail values are the same" + ) + expect_equal(val$pareto_k, Inf) }) From b16435a32167b7a68a03a496995cc1744dd951f2 Mon Sep 17 00:00:00 2001 From: VisruthSK Date: Wed, 8 Apr 2026 16:06:14 -0700 Subject: [PATCH 11/11] Reverted tail logic/warning handling --- R/psis.R | 38 ++++++++++++++++++++++++-------------- 1 file changed, 24 insertions(+), 14 deletions(-) diff --git a/R/psis.R b/R/psis.R index 639f462f..eb6fb647 100644 --- a/R/psis.R +++ b/R/psis.R @@ -212,24 +212,34 @@ do_psis_i <- function(log_ratios_i, tail_len_i, ...) { S <- length(log_ratios_i) # shift log ratios for safer exponentation lw_i <- log_ratios_i - max(log_ratios_i) + khat <- Inf + smooth_tail <- TRUE - if (length(unique(utils::tail(sort(log_ratios_i), tail_len_i))) == 1) { - warning( - "Can't fit generalized Pareto distribution ", - "because all tail values are the same.", - call. = FALSE - ) + if ( + enough_tail_samples(tail_len_i) + ) { + lw_tail <- utils::tail(sort.int(lw_i), tail_len_i) + if (abs(max(lw_tail) - min(lw_tail)) < .Machine$double.eps / 100) { + warning( + "Can't fit generalized Pareto distribution ", + "because all tail values are the same.", + call. = FALSE + ) + smooth_tail <- FALSE + } } - smoothed <- suppressWarnings(posterior::ps_tail( - x = lw_i, - ndraws_tail = tail_len_i, - tail = "right", - are_log_weights = TRUE - )) + if (smooth_tail) { + smoothed <- suppressWarnings(posterior::ps_tail( + x = lw_i, + ndraws_tail = tail_len_i, + tail = "right", + are_log_weights = TRUE + )) - lw_i <- smoothed$x - khat <- smoothed$k + lw_i <- smoothed$x + khat <- smoothed$k + } # truncate at max of raw wts (i.e., 0 since max has been subtracted) lw_i[lw_i > 0] <- 0