From 2af77630e0d5ede38878a256e7f435d3ed06d7c9 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Wed, 21 Oct 2026 12:43:13 +0300
Subject: [PATCH 1/2] Sort only the tail in do_psis_i()

---
 R/psis.R                   | 16 +++++++++++-----
 tests/testthat/test_psis.R | 20 ++++++++++++++++++++
 2 files changed, 31 insertions(+), 5 deletions(-)

diff --git a/R/psis.R b/R/psis.R
index 303da755..46c31f8b 100644
--- a/R/psis.R
+++ b/R/psis.R
@@ -215,9 +215,15 @@ do_psis_i <- function(log_ratios_i, tail_len_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]
+    # only the tail and the cutoff need sorting: a partial sort finds the
+    # cutoff value, and the values at or above it are few
+    cut_pos <- S - tail_len_i
+    cut_val <- sort.int(lw_i, partial = cut_pos)[cut_pos]
+    top <- which(lw_i >= cut_val)
+    top <- top[order(lw_i[top])]
+    top <- top[seq(length(top) - tail_len_i, length(top))] # cutoff and tail
+    tail_ids <- top[-1]
+    lw_tail <- lw_i[tail_ids]
     if (abs(max(lw_tail) - min(lw_tail)) < .Machine$double.eps / 100) {
       warning(
         "Can't fit generalized Pareto distribution ",
@@ -225,10 +231,10 @@ do_psis_i <- function(log_ratios_i, tail_len_i, ...) {
         call. = FALSE
       )
     } else {
-      cutoff <- ord$x[min(tail_ids) - 1] # largest value smaller than tail values
+      cutoff <- lw_i[top[1]] # largest value not in the tail
       smoothed <- psis_smooth_tail(lw_tail, cutoff)
       khat <- smoothed$k
-      lw_i[ord$ix[tail_ids]] <- smoothed$tail
+      lw_i[tail_ids] <- smoothed$tail
     }
   }
 
diff --git a/tests/testthat/test_psis.R b/tests/testthat/test_psis.R
index 9e033715..cd63699d 100644
--- a/tests/testthat/test_psis.R
+++ b/tests/testthat/test_psis.R
@@ -178,3 +178,23 @@ test_that("psis_smooth_tail returns original tail values if k is infinite", {
   expect_equal(val$k, Inf)
 })
 
+test_that("do_psis_i selects the same tail as a full sort, also with ties", {
+  # reference: the tail from a full sort, as in loo 2.10.1
+  do_psis_i_full_sort <- function(log_ratios_i, tail_len_i) {
+    S <- length(log_ratios_i)
+    lw_i <- log_ratios_i - max(log_ratios_i)
+    ord <- sort.int(lw_i, index.return = TRUE)
+    tail_ids <- seq(S - tail_len_i + 1, S)
+    smoothed <- psis_smooth_tail(ord$x[tail_ids], ord$x[min(tail_ids) - 1])
+    lw_i[ord$ix[tail_ids]] <- smoothed$tail
+    lw_i[lw_i > 0] <- 0
+    list(log_weights = lw_i + max(log_ratios_i), pareto_k = smoothed$k)
+  }
+  set.seed(260)
+  # rounding gives many ties, also at the cutoff
+  for (x in list(round(rnorm(4000), 1), round(rt(1000, df = 2), 2))) {
+    tail_len <- n_pareto(1, length(x))
+    expect_identical(do_psis_i(x, tail_len), do_psis_i_full_sort(x, tail_len))
+  }
+})
+

From 328af193ea82bcf0e1eb9716edcfa54d30d2c37b Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Wed, 21 Oct 2026 20:45:06 +0300
Subject: [PATCH 2/2] Use the partial-sort cutoff value directly in do_psis_i()

---
 R/psis.R | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/R/psis.R b/R/psis.R
index 46c31f8b..85f8883c 100644
--- a/R/psis.R
+++ b/R/psis.R
@@ -216,13 +216,12 @@ do_psis_i <- function(log_ratios_i, tail_len_i, ...) {
 
   if (enough_tail_samples(tail_len_i)) {
     # only the tail and the cutoff need sorting: a partial sort finds the
-    # cutoff value, and the values at or above it are few
+    # cutoff value, and the values at or above it are usually few
     cut_pos <- S - tail_len_i
     cut_val <- sort.int(lw_i, partial = cut_pos)[cut_pos]
     top <- which(lw_i >= cut_val)
     top <- top[order(lw_i[top])]
-    top <- top[seq(length(top) - tail_len_i, length(top))] # cutoff and tail
-    tail_ids <- top[-1]
+    tail_ids <- top[seq(length(top) - tail_len_i + 1, length(top))]
     lw_tail <- lw_i[tail_ids]
     if (abs(max(lw_tail) - min(lw_tail)) < .Machine$double.eps / 100) {
       warning(
@@ -231,7 +230,7 @@ do_psis_i <- function(log_ratios_i, tail_len_i, ...) {
         call. = FALSE
       )
     } else {
-      cutoff <- lw_i[top[1]] # largest value not in the tail
+      cutoff <- cut_val # largest value not in the tail
       smoothed <- psis_smooth_tail(lw_tail, cutoff)
       khat <- smoothed$k
       lw_i[tail_ids] <- smoothed$tail