From d88bbaf56ef6eff45d75aea879798e1a9c8a59a3 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Sat, 12 Sep 2026 23:56:20 +0300
Subject: [PATCH 01/22] fix: preserve extreme GPD tail probabilities

---
 R/gpd.R                                 | 27 +++++++++++--------------
 tests/testthat/test-gpd-distributions.R |  5 +++++
 2 files changed, 17 insertions(+), 15 deletions(-)

diff --git a/R/gpd.R b/R/gpd.R
index a9f2e1c8..ac7e99c9 100644
--- a/R/gpd.R
+++ b/R/gpd.R
@@ -20,15 +20,14 @@ qgeneralized_pareto <- function(p, mu = 0, sigma = 1, k = 0, lower.tail = TRUE,
     return(rep(NaN, length(p)))
   }
   if (log.p) {
-    p <- exp(p)
-  }
-  if (!lower.tail) {
-    p <- 1 - p
+    log_survival <- if (lower.tail) log(-expm1(p)) else p
+  } else {
+    log_survival <- if (lower.tail) log1p(-p) else log(p)
   }
   if (k == 0) {
-    q <-  mu - sigma * log1p(-p)
+    q <- mu - sigma * log_survival
   } else {
-    q <- mu + sigma * expm1(-k * log1p(-p)) / k
+    q <- mu + sigma * expm1(-k * log_survival) / k
   }
   q
 }
@@ -57,18 +56,16 @@ pgeneralized_pareto <- function(q, mu = 0, sigma = 1, k = 0, lower.tail = TRUE,
   z <- (q - mu) / sigma
   if (abs(k) < 1e-15) {
     # for very small values of indistinguishable in floating point accuracy from the case k=0
-    p <- -expm1(-z)
+    log_survival <- -z
   } else {
     # pmax handles values outside the support
-    p <- -expm1(log1p(pmax(k * z, -1)) / -k)
-  }
-  # force to [0, 1] for values outside the support
-  p <- pmin(pmax(p, 0), 1)
-  if (!lower.tail) {
-    p <- 1 - p
+    log_survival <- log1p(pmax(k * z, -1)) / -k
   }
-  if (log.p) {
-    p <- log(p)
+  log_survival <- pmin(log_survival, 0)
+  if (lower.tail) {
+    p <- if (log.p) log(-expm1(log_survival)) else -expm1(log_survival)
+  } else {
+    p <- if (log.p) log_survival else exp(log_survival)
   }
   p
 }
diff --git a/tests/testthat/test-gpd-distributions.R b/tests/testthat/test-gpd-distributions.R
index d1cbe01c..76c3deea 100644
--- a/tests/testthat/test-gpd-distributions.R
+++ b/tests/testthat/test-gpd-distributions.R
@@ -42,6 +42,9 @@ test_that("qgeneralized_pareto handles log.p = TRUE", {
   result <- qgeneralized_pareto(p, mu = 0, sigma = 1, k = 0.2)
   result_log <- qgeneralized_pareto(log(p), mu = 0, sigma = 1, k = 0.2, log.p = TRUE)
   expect_equal(result, result_log)
+
+  result_upper <- qgeneralized_pareto(-100, k = 0.2, lower.tail = FALSE, log.p = TRUE)
+  expect_equal(result_upper, expm1(20) / 0.2)
 })
 
 test_that("qgeneralized_pareto returns NaN for invalid sigma", {
@@ -117,6 +120,8 @@ test_that("pgeneralized_pareto handles lower.tail = FALSE", {
   result_lower <- pgeneralized_pareto(q, mu = 0, sigma = 1, k = 0.2)
   result_upper <- pgeneralized_pareto(q, mu = 0, sigma = 1, k = 0.2, lower.tail = FALSE)
   expect_equal(result_lower + result_upper, rep(1, 3))
+  expect_equal(pgeneralized_pareto(100, k = 0, lower.tail = FALSE), exp(-100))
+  expect_equal(pgeneralized_pareto(100, k = 0, lower.tail = FALSE, log.p = TRUE), -100)
 })
 
 test_that("pgeneralized_pareto handles log.p = TRUE", {

From 01f464a49fed600f76b445e6acfb96837de10d65 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Sat, 12 Sep 2026 23:58:00 +0300
Subject: [PATCH 02/22] fix: stabilize Pareto convergence rates

---
 R/pareto_smooth.R                   | 31 ++++++++++++++++++-----------
 tests/testthat/test-pareto_smooth.R | 15 ++++++++++++++
 2 files changed, 34 insertions(+), 12 deletions(-)

diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R
index b82accec..a989e6c8 100644
--- a/R/pareto_smooth.R
+++ b/R/pareto_smooth.R
@@ -633,20 +633,27 @@ ps_khat_threshold <- function(ndraws, ...) {
 ps_convergence_rate <- function(k, ndraws, ...) {
   # allow array of k's
   rate <- numeric(length(k))
-  # k<0 bounded distribution
-  rate[k < 0] <- 1
-  # k>0 non-finite mean
-  rate[k > 1] <- 0
-  # limit value at k=1/2
-  rate[k == 0.5] <- 1 - 1 / log(ndraws)
+  # k<=0 bounded distribution
+  rate[k <= 0] <- 1
+  # k>=1 non-finite mean
+  rate[k >= 1] <- 0
   # smooth approximation for the rest (see Appendix B of PSIS paper)
-  ki <- (k > 0 & k < 1 & k != 0.5)
+  ki <- k > 0 & k < 1
   kk <- k[ki]
-  rate[ki] <- pmax(
-    0,
-    (2 * (kk - 1) * ndraws^(2 * kk + 1) + (1 - 2 * kk) * ndraws^(2 * kk) + ndraws^2) /
-      ((ndraws - 1) * (ndraws - ndraws^(2 * kk)))
-  )
+  aa <- 2 * kk - 1
+  near_half <- abs(aa) <= abs(1 - aa)
+  rate_ki <- numeric(length(kk))
+  rate_ki[aa == 0] <- ndraws / (ndraws - 1) - 1 / log(ndraws)
+  use_half <- near_half & aa != 0
+  rate_ki[use_half] <-
+    ndraws / (ndraws - 1) -
+    aa[use_half] / -expm1(-aa[use_half] * log(ndraws))
+  use_one <- !near_half
+  bb <- 1 - aa[use_one]
+  rate_ki[use_one] <-
+    (bb * (ndraws - 1) - expm1(bb * log(ndraws))) /
+    ((ndraws - 1) * -expm1((bb - 1) * log(ndraws)))
+  rate[ki] <- pmax(0, rate_ki)
   rate
 }
 
diff --git a/tests/testthat/test-pareto_smooth.R b/tests/testthat/test-pareto_smooth.R
index d9dc10d1..9e6dd342 100644
--- a/tests/testthat/test-pareto_smooth.R
+++ b/tests/testthat/test-pareto_smooth.R
@@ -294,3 +294,18 @@ test_that("check ps_min_ss behavior special cases", {
   # k < 1
   expect_equal(ps_min_ss(0.5), 10^(1 / (1 - max(0, 0.5))))
 })
+
+
+test_that("ps_convergence_rate is stable at transition points", {
+  n <- 10
+  half_limit <- n / (n - 1) - 1 / log(n)
+
+  expect_equal(ps_convergence_rate(0, n), 1)
+  expect_equal(ps_convergence_rate(0.5, n), half_limit)
+  expect_equal(ps_convergence_rate(0.5 + .Machine$double.eps / 2, n), half_limit)
+  expect_equal(
+    ps_convergence_rate(1 - .Machine$double.eps / 2, n),
+    1.8359566012903116e-16
+  )
+  expect_equal(ps_convergence_rate(1, n), 0)
+})

From 2816f4c91708d740187a4b0482c7c2f3a156b562 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Sat, 12 Sep 2026 23:59:20 +0300
Subject: [PATCH 03/22] fix: stabilize Pareto tail differences

---
 R/pareto_smooth.R                   | 11 +++++++++--
 tests/testthat/test-pareto_smooth.R |  9 +++++++++
 2 files changed, 18 insertions(+), 2 deletions(-)

diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R
index a989e6c8..434d0ef9 100644
--- a/R/pareto_smooth.R
+++ b/R/pareto_smooth.R
@@ -454,6 +454,11 @@ pareto_convergence_rate.rvar <- function(x, ...) {
 }
 
 
+exp_x_minus_exp_y <- function(x, y) {
+  -exp(x) * expm1(y - x)
+}
+
+
 #' Pareto smooth tail
 #' function to Pareto smooth the tail of a vector. Exported
 #' for usage in other packages, not by users.
@@ -525,10 +530,12 @@ ps_tail <- function(x,
   max_tail <- max(draws_tail)
 
   if (are_log_weights) {
-    draws_tail <- exp(draws_tail)
+    tail_excesses <- exp_x_minus_exp_y(draws_tail, cutoff)
     cutoff <- exp(cutoff)
+  } else {
+    tail_excesses <- draws_tail - cutoff
   }
-  fit <- gpdfit(draws_tail - cutoff, sort_x = FALSE, ...)
+  fit <- gpdfit(tail_excesses, sort_x = FALSE, ...)
   k <- fit$k
   sigma <- fit$sigma
   if (is.finite(k) && smooth_draws) {
diff --git a/tests/testthat/test-pareto_smooth.R b/tests/testthat/test-pareto_smooth.R
index 9e6dd342..7caa2b81 100644
--- a/tests/testthat/test-pareto_smooth.R
+++ b/tests/testthat/test-pareto_smooth.R
@@ -261,6 +261,15 @@ test_that("pareto_smooth works for log_weights", {
 
 })
 
+test_that("exp_x_minus_exp_y is stable for nearby values", {
+  cutoff <- -30
+  x <- cutoff + 1e-14
+
+  expect_equal(exp_x_minus_exp_y(x, cutoff), 9.973486536736925e-28)
+  expect_equal(exp_x_minus_exp_y(0, -1000), 1)
+})
+
+
 test_that("check ps_tail behavior for ndraws_tail less than 5", {
   w <- c(1:25, 1e3, 1e3, 1e3)
   lw <- log(w)

From 5cdf5ffab83aeb36cfe1910f19043f39ea4d5179 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Sun, 13 Sep 2026 00:00:43 +0300
Subject: [PATCH 04/22] fix: stabilize MCSE standard deviation

---
 R/convergence.R                   | 7 ++++---
 tests/testthat/test-convergence.R | 5 +++++
 2 files changed, 9 insertions(+), 3 deletions(-)

diff --git a/R/convergence.R b/R/convergence.R
index ecd74b40..894c6bef 100644
--- a/R/convergence.R
+++ b/R/convergence.R
@@ -495,11 +495,12 @@ mcse_sd.default <- function(x, ...) {
   # has (X-E[X])^2. The following ESS is based on a relevant quantity
   # in the computation and is empirically a good choice.
   sims_c <- x - mean(x)
-  ess <- ess_mean((sims_c)^2)
+  sims_sq <- sims_c^2
+  ess <- ess_mean(sims_sq)
   # Variance of variance estimate by Kenney and Keeping (1951, p. 141),
   # which doesn't assume normality of sims.
-  Evar <- mean(sims_c^2)
-  varvar <- (mean(sims_c^4) - Evar^2) / ess
+  Evar <- mean(sims_sq)
+  varvar <- mean((sims_sq - Evar)^2) / ess
   # The first order Taylor series approximation of variance of sd.
   # Kenney and Keeping (1951, p. 141) write "...since fluctuations of
   # any moment are of order N^{-1/2}, squares and higher powers of
diff --git a/tests/testthat/test-convergence.R b/tests/testthat/test-convergence.R
index 0da94f5b..787f18c5 100644
--- a/tests/testthat/test-convergence.R
+++ b/tests/testthat/test-convergence.R
@@ -55,6 +55,11 @@ test_that("mcse diagnostics return reasonable values", {
   mcse <- mcse_sd(tau)
   expect_true(mcse > 0.15 & mcse < 0.25)
 
+  x <- rep(c(-1e8 - 1, -1e8 + 1, 1e8 - 1, 1e8 + 1), 100)
+  x_sq <- (x - mean(x))^2
+  expected <- sqrt(mean((x_sq - mean(x_sq))^2) / ess_mean(x_sq) / mean(x_sq) / 4)
+  expect_equal(mcse_sd(x), expected)
+
   mcse <- mcse_quantile(tau, probs = c(0.2, 0.8))
   expect_equal(names(mcse), c("mcse_q20", "mcse_q80"))
   expect_true(mcse[1] > 0.16 & mcse[1] < 0.21)

From a7aea6f036c724e706430a51876a1cc8fa524362 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Sun, 13 Sep 2026 00:02:00 +0300
Subject: [PATCH 05/22] fix: stabilize dissent logarithms

---
 R/discrete-summaries.R                   | 7 ++++++-
 tests/testthat/test-discrete-summaries.R | 4 ++++
 2 files changed, 10 insertions(+), 1 deletion(-)

diff --git a/R/discrete-summaries.R b/R/discrete-summaries.R
index 3c141c9d..26e56ef0 100644
--- a/R/discrete-summaries.R
+++ b/R/discrete-summaries.R
@@ -76,6 +76,11 @@ entropy.rvar <- function(x) {
 }
 
 
+log2_one_minus <- function(x) {
+  log1p(-x) / log(2)
+}
+
+
 #' Dissention
 #'
 #' Dissention, for measuring dispersion in draws from ordinal distributions.
@@ -148,7 +153,7 @@ dissent.default <- function(x) {
     out <- 0
   } else {
     x_i <- tab$x
-    out <- -sum(p * log2(1 - abs(x_i - mean(x)) / d))
+    out <- -sum(p * log2_one_minus(abs(x_i - mean(x)) / d))
   }
   out
 }
diff --git a/tests/testthat/test-discrete-summaries.R b/tests/testthat/test-discrete-summaries.R
index 78ce902b..b8496c67 100644
--- a/tests/testthat/test-discrete-summaries.R
+++ b/tests/testthat/test-discrete-summaries.R
@@ -57,6 +57,10 @@ test_that("entropy works on rvars", {
 
 # dissent -----------------------------------------------------------------
 
+test_that("log2_one_minus is stable near zero", {
+  expect_equal(log2_one_minus(1e-20), -1.4426950408889633e-20)
+})
+
 test_that("dissent works on vectors", {
   expect_equal(dissent(NULL), 0)
   expect_equal(dissent(1), 0)

From 96d8a950e59aa29a84412bffa6ffe86ea2c03ec6 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Sun, 13 Sep 2026 10:15:51 +0300
Subject: [PATCH 06/22] refactor: move new helper functions to misc.R

---
 R/discrete-summaries.R |  5 -----
 R/misc.R               | 10 ++++++++++
 R/pareto_smooth.R      |  5 -----
 3 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/R/discrete-summaries.R b/R/discrete-summaries.R
index 26e56ef0..937e2a55 100644
--- a/R/discrete-summaries.R
+++ b/R/discrete-summaries.R
@@ -76,11 +76,6 @@ entropy.rvar <- function(x) {
 }
 
 
-log2_one_minus <- function(x) {
-  log1p(-x) / log(2)
-}
-
-
 #' Dissention
 #'
 #' Dissention, for measuring dispersion in draws from ordinal distributions.
diff --git a/R/misc.R b/R/misc.R
index 1027ff9d..c36ca494 100644
--- a/R/misc.R
+++ b/R/misc.R
@@ -233,6 +233,16 @@ log_sum_exp <- function(x) {
   res
 }
 
+# numerically stable version of exp(x) - exp(y)
+exp_x_minus_exp_y <- function(x, y) {
+  -exp(x) * expm1(y - x)
+}
+
+# numerically stable version of log2(1 - x)
+log2_one_minus <- function(x) {
+  log1p(-x) / log(2)
+}
+
 # simple version of destructuring assignment
 `%<-%` <- function(vars, values, envir = parent.frame()) {
   vars <- as.character(substitute(vars)[-1])
diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R
index 434d0ef9..80070aff 100644
--- a/R/pareto_smooth.R
+++ b/R/pareto_smooth.R
@@ -454,11 +454,6 @@ pareto_convergence_rate.rvar <- function(x, ...) {
 }
 
 
-exp_x_minus_exp_y <- function(x, y) {
-  -exp(x) * expm1(y - x)
-}
-
-
 #' Pareto smooth tail
 #' function to Pareto smooth the tail of a vector. Exported
 #' for usage in other packages, not by users.

From a7799422ee5fd47c43aa106f79c217b2d770d7b3 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <aki.vehtari@aalto.fi>
Date: Mon, 14 Sep 2026 11:17:11 +0300
Subject: [PATCH 07/22] Add clarifying comment

Co-authored-by: Noa Kallioinen <33577035+n-kall@users.noreply.github.com>
---
 R/pareto_smooth.R | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R
index 80070aff..7d79d690 100644
--- a/R/pareto_smooth.R
+++ b/R/pareto_smooth.R
@@ -639,7 +639,8 @@ ps_convergence_rate <- function(k, ndraws, ...) {
   rate[k <= 0] <- 1
   # k>=1 non-finite mean
   rate[k >= 1] <- 0
-  # smooth approximation for the rest (see Appendix B of PSIS paper)
+  # smooth approximation for the rest 
+  # this is a numerically stable calculation of the equation on page 31, in Appendix B of PSIS paper
   ki <- k > 0 & k < 1
   kk <- k[ki]
   aa <- 2 * kk - 1

From d22e9a04850eaf1dd4e64690947b613f0f65263a Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Mon, 14 Sep 2026 11:26:45 +0300
Subject: [PATCH 08/22] fix: distinguish small nonzero GPD shapes

---
 R/gpd.R                                 |  3 +--
 tests/testthat/test-gpd-distributions.R | 15 +++++++++++++++
 2 files changed, 16 insertions(+), 2 deletions(-)

diff --git a/R/gpd.R b/R/gpd.R
index ac7e99c9..cfc193a5 100644
--- a/R/gpd.R
+++ b/R/gpd.R
@@ -54,8 +54,7 @@ pgeneralized_pareto <- function(q, mu = 0, sigma = 1, k = 0, lower.tail = TRUE,
     return(rep(NaN, length(q)))
   }
   z <- (q - mu) / sigma
-  if (abs(k) < 1e-15) {
-    # for very small values of indistinguishable in floating point accuracy from the case k=0
+  if (k == 0) {
     log_survival <- -z
   } else {
     # pmax handles values outside the support
diff --git a/tests/testthat/test-gpd-distributions.R b/tests/testthat/test-gpd-distributions.R
index 76c3deea..72f8c583 100644
--- a/tests/testthat/test-gpd-distributions.R
+++ b/tests/testthat/test-gpd-distributions.R
@@ -124,6 +124,21 @@ test_that("pgeneralized_pareto handles lower.tail = FALSE", {
   expect_equal(pgeneralized_pareto(100, k = 0, lower.tail = FALSE, log.p = TRUE), -100)
 })
 
+test_that("pgeneralized_pareto distinguishes small nonzero shape", {
+  k <- 9e-16
+  q <- 1e7
+  expected <- -log1p(k * q) / k
+  expect_equal(
+    pgeneralized_pareto(q, k = k, lower.tail = FALSE, log.p = TRUE),
+    expected
+  )
+
+  expect_equal(
+    pgeneralized_pareto(2e15, k = -k, lower.tail = FALSE, log.p = TRUE),
+    -Inf
+  )
+})
+
 test_that("pgeneralized_pareto handles log.p = TRUE", {
   q <- c(0.5, 1, 2)
   result <- pgeneralized_pareto(q, mu = 0, sigma = 1, k = 0.2)

From eecfde877bc774887fadb05d17b611733fd77fd9 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Tue, 15 Sep 2026 00:14:47 +0300
Subject: [PATCH 09/22] fix: return zero from exp_x_minus_exp_y for equal
 infinite inputs

---
 R/misc.R                            |  6 +++++-
 tests/testthat/test-pareto_smooth.R | 16 ++++++++++++++++
 2 files changed, 21 insertions(+), 1 deletion(-)

diff --git a/R/misc.R b/R/misc.R
index c36ca494..202fc9b8 100644
--- a/R/misc.R
+++ b/R/misc.R
@@ -235,7 +235,11 @@ log_sum_exp <- function(x) {
 
 # numerically stable version of exp(x) - exp(y)
 exp_x_minus_exp_y <- function(x, y) {
-  -exp(x) * expm1(y - x)
+  out <- -exp(x) * expm1(y - x)
+  # equal infinite inputs give 0 * NaN above; which() drops the NA comparisons
+  # that NA or NaN inputs would produce
+  out[which(x == y)] <- 0
+  out
 }
 
 # numerically stable version of log2(1 - x)
diff --git a/tests/testthat/test-pareto_smooth.R b/tests/testthat/test-pareto_smooth.R
index 7caa2b81..48aeb8dc 100644
--- a/tests/testthat/test-pareto_smooth.R
+++ b/tests/testthat/test-pareto_smooth.R
@@ -269,6 +269,22 @@ test_that("exp_x_minus_exp_y is stable for nearby values", {
   expect_equal(exp_x_minus_exp_y(0, -1000), 1)
 })
 
+test_that("exp_x_minus_exp_y returns zero for equal infinite inputs", {
+  expect_equal(exp_x_minus_exp_y(-Inf, -Inf), 0)
+  expect_equal(exp_x_minus_exp_y(c(-Inf, 0, 1), c(-Inf, 0, 0)), c(0, 0, exp(1) - 1))
+  expect_equal(exp_x_minus_exp_y(0, -Inf), 1)
+})
+
+test_that("ps_tail handles -Inf log weights below the cutoff", {
+  # more tail draws than finite log weights, so the cutoff is -Inf and the
+  # lowest tail draws equal it
+  lw <- c(rep(-Inf, 90), log(1:10))
+  tail <- ps_tail(lw, ndraws_tail = 20, tail = "right", are_log_weights = TRUE)
+
+  expect_false(anyNA(tail$x))
+  expect_true(all(is.infinite(tail$x[1:80]) & tail$x[1:80] < 0))
+})
+
 
 test_that("check ps_tail behavior for ndraws_tail less than 5", {
   w <- c(1:25, 1e3, 1e3, 1e3)

From 0542ad8ab62751d44c73981be77d24bd2f84a2c0 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Tue, 15 Sep 2026 00:15:30 +0300
Subject: [PATCH 10/22] fix: propagate missing shape parameters through
 ps_convergence_rate

---
 R/pareto_smooth.R                   | 15 ++++++++-------
 tests/testthat/test-pareto_smooth.R | 11 +++++++++++
 2 files changed, 19 insertions(+), 7 deletions(-)

diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R
index 7d79d690..70261de0 100644
--- a/R/pareto_smooth.R
+++ b/R/pareto_smooth.R
@@ -633,15 +633,16 @@ ps_khat_threshold <- function(ndraws, ...) {
 #' @return convergence rate
 #' @export
 ps_convergence_rate <- function(k, ndraws, ...) {
-  # allow array of k's
-  rate <- numeric(length(k))
-  # k<=0 bounded distribution
-  rate[k <= 0] <- 1
+  # allow array of k's; which() keeps NA and NaN shapes out of every branch,
+  # so they keep the NA initialization instead of silently reading as 0
+  rate <- rep(NA_real_, length(k))
+  # k<0 bounded distribution, k=0 exponential tail; both converge at rate 1
+  rate[which(k <= 0)] <- 1
   # k>=1 non-finite mean
-  rate[k >= 1] <- 0
-  # smooth approximation for the rest 
+  rate[which(k >= 1)] <- 0
+  # smooth approximation for the rest
   # this is a numerically stable calculation of the equation on page 31, in Appendix B of PSIS paper
-  ki <- k > 0 & k < 1
+  ki <- which(k > 0 & k < 1)
   kk <- k[ki]
   aa <- 2 * kk - 1
   near_half <- abs(aa) <= abs(1 - aa)
diff --git a/tests/testthat/test-pareto_smooth.R b/tests/testthat/test-pareto_smooth.R
index 48aeb8dc..f95bdb7e 100644
--- a/tests/testthat/test-pareto_smooth.R
+++ b/tests/testthat/test-pareto_smooth.R
@@ -334,3 +334,14 @@ test_that("ps_convergence_rate is stable at transition points", {
   )
   expect_equal(ps_convergence_rate(1, n), 0)
 })
+
+test_that("ps_convergence_rate propagates missing shape parameters", {
+  n <- 1000
+
+  expect_identical(ps_convergence_rate(NA_real_, n), NA_real_)
+  expect_identical(ps_convergence_rate(NaN, n), NA_real_)
+  expect_equal(
+    ps_convergence_rate(c(0.3, NA, 0.8, NaN, -1, 2), n),
+    c(ps_convergence_rate(0.3, n), NA, ps_convergence_rate(0.8, n), NA, 1, 0)
+  )
+})

From 61cdaf2b1ed364c157b0332a6399062ee999b21d Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Tue, 15 Sep 2026 00:16:52 +0300
Subject: [PATCH 11/22] fix: reject out-of-range probabilities in
 qgeneralized_pareto

---
 R/gpd.R                                 | 12 ++++++-
 man/qgeneralized_pareto.Rd              |  4 ++-
 tests/testthat/test-gpd-distributions.R | 46 +++++++++++++++++++++++++
 3 files changed, 60 insertions(+), 2 deletions(-)

diff --git a/R/gpd.R b/R/gpd.R
index cfc193a5..57eeaa5b 100644
--- a/R/gpd.R
+++ b/R/gpd.R
@@ -3,7 +3,9 @@
 #' Computes the quantile function for a generalized Pareto distribution
 #' with location `mu`, scale `sigma`, and shape `k`.
 #' 
-#' @param p Numeric vector of probabilities.
+#' @param p Numeric vector of probabilities, in `[0, 1]`, or of log
+#'   probabilities, in `[-Inf, 0]`, if `log.p` is `TRUE`. Values outside the
+#'   valid range return `NaN` with a warning.
 #' @param mu Location parameter.
 #' @param sigma Scale parameter (must be positive).
 #' @param k Shape parameter.
@@ -19,6 +21,14 @@ qgeneralized_pareto <- function(p, mu = 0, sigma = 1, k = 0, lower.tail = TRUE,
   if (is.na(sigma) || sigma <= 0) {
     return(rep(NaN, length(p)))
   }
+  # probabilities outside the valid range give NaN with a warning, as in the
+  # base R quantile functions. which() leaves NA and NaN inputs to propagate,
+  # and replacing them here avoids a second warning from log() or log1p().
+  invalid <- which(if (log.p) p > 0 else p < 0 | p > 1)
+  if (length(invalid)) {
+    warning_no_call("NaNs produced")
+    p[invalid] <- NaN
+  }
   if (log.p) {
     log_survival <- if (lower.tail) log(-expm1(p)) else p
   } else {
diff --git a/man/qgeneralized_pareto.Rd b/man/qgeneralized_pareto.Rd
index d9cda2ce..17a8261a 100644
--- a/man/qgeneralized_pareto.Rd
+++ b/man/qgeneralized_pareto.Rd
@@ -14,7 +14,9 @@ qgeneralized_pareto(
 )
 }
 \arguments{
-\item{p}{Numeric vector of probabilities.}
+\item{p}{Numeric vector of probabilities, in \code{[0, 1]}, or of log
+probabilities, in \code{[-Inf, 0]}, if \code{log.p} is \code{TRUE}. Values outside the
+valid range return \code{NaN} with a warning.}
 
 \item{mu}{Location parameter.}
 
diff --git a/tests/testthat/test-gpd-distributions.R b/tests/testthat/test-gpd-distributions.R
index 72f8c583..3d636785 100644
--- a/tests/testthat/test-gpd-distributions.R
+++ b/tests/testthat/test-gpd-distributions.R
@@ -47,6 +47,52 @@ test_that("qgeneralized_pareto handles log.p = TRUE", {
   expect_equal(result_upper, expm1(20) / 0.2)
 })
 
+test_that("qgeneralized_pareto returns the endpoints for infinite log probabilities", {
+  # log.p = TRUE with p = -Inf is probability 0: the lower endpoint for the
+  # lower tail, the upper endpoint of the support for the upper tail
+  for (k in c(-0.4, 0, 0.3)) {
+    expect_equal(
+      qgeneralized_pareto(-Inf, mu = 2, sigma = 3, k = k, log.p = TRUE),
+      2
+    )
+    upper <- if (k < 0) 2 - 3 / k else Inf
+    expect_equal(
+      qgeneralized_pareto(-Inf, mu = 2, sigma = 3, k = k,
+                          lower.tail = FALSE, log.p = TRUE),
+      upper
+    )
+  }
+
+  # log.p = TRUE with p = 0 is probability 1, i.e. the two endpoints swapped
+  expect_equal(qgeneralized_pareto(0, mu = 2, sigma = 3, k = -0.4, log.p = TRUE), 9.5)
+  expect_equal(
+    qgeneralized_pareto(0, mu = 2, sigma = 3, k = -0.4,
+                        lower.tail = FALSE, log.p = TRUE),
+    2
+  )
+})
+
+test_that("qgeneralized_pareto rejects probabilities outside the valid range", {
+  for (lower in c(TRUE, FALSE)) {
+    expect_warning(
+      out <- qgeneralized_pareto(c(-0.5, 0.5, 1.5), k = 0.2, lower.tail = lower),
+      "NaNs produced"
+    )
+    expect_equal(is.nan(out), c(TRUE, FALSE, TRUE))
+
+    expect_warning(
+      out <- qgeneralized_pareto(c(1, -1), k = 0.2, lower.tail = lower, log.p = TRUE),
+      "NaNs produced"
+    )
+    expect_equal(is.nan(out), c(TRUE, FALSE))
+  }
+
+  # missing values stay missing rather than becoming NaN
+  expect_identical(qgeneralized_pareto(NA_real_, k = 0.2), NA_real_)
+  expect_silent(qgeneralized_pareto(c(0, 0.5, 1), k = 0.2))
+  expect_silent(qgeneralized_pareto(c(-Inf, -1, 0), k = 0.2, log.p = TRUE))
+})
+
 test_that("qgeneralized_pareto returns NaN for invalid sigma", {
   result <- qgeneralized_pareto(0.5, mu = 0, sigma = -1, k = 0.2)
   expect_true(is.nan(result))

From fab0f076849bf43bba4621f2ddaf92f671c42f24 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Tue, 15 Sep 2026 00:17:51 +0300
Subject: [PATCH 12/22] test: cover -Inf contracts for log_sum_exp, GPD
 endpoints and pareto draws

---
 tests/testthat/test-gpd-distributions.R | 34 +++++++++++++++++++++++++
 tests/testthat/test-log_sum_exp.R       |  9 +++++++
 tests/testthat/test-pareto_smooth.R     | 19 ++++++++++++++
 3 files changed, 62 insertions(+)

diff --git a/tests/testthat/test-gpd-distributions.R b/tests/testthat/test-gpd-distributions.R
index 3d636785..b019d3c3 100644
--- a/tests/testthat/test-gpd-distributions.R
+++ b/tests/testthat/test-gpd-distributions.R
@@ -170,6 +170,40 @@ test_that("pgeneralized_pareto handles lower.tail = FALSE", {
   expect_equal(pgeneralized_pareto(100, k = 0, lower.tail = FALSE, log.p = TRUE), -100)
 })
 
+test_that("pgeneralized_pareto is -Inf at mu and below the support", {
+  for (k in c(-0.4, 0, 0.3)) {
+    # at the lower endpoint the lower-tail probability is exactly 0
+    expect_equal(pgeneralized_pareto(2, mu = 2, sigma = 3, k = k), 0)
+    expect_equal(
+      pgeneralized_pareto(2, mu = 2, sigma = 3, k = k, log.p = TRUE),
+      -Inf
+    )
+    expect_equal(
+      pgeneralized_pareto(2, mu = 2, sigma = 3, k = k, lower.tail = FALSE, log.p = TRUE),
+      0
+    )
+
+    # below the support the same holds, without escaping [0, 1]
+    below <- c(-Inf, -1e300, 1.5)
+    expect_equal(pgeneralized_pareto(below, mu = 2, sigma = 3, k = k), rep(0, 3))
+    expect_equal(
+      pgeneralized_pareto(below, mu = 2, sigma = 3, k = k, log.p = TRUE),
+      rep(-Inf, 3)
+    )
+    expect_equal(
+      pgeneralized_pareto(below, mu = 2, sigma = 3, k = k, lower.tail = FALSE),
+      rep(1, 3)
+    )
+  }
+
+  # above the support of a negative shape the upper tail is exactly 0
+  expect_equal(
+    pgeneralized_pareto(c(9.5, 20, Inf), mu = 2, sigma = 3, k = -0.4,
+                        lower.tail = FALSE, log.p = TRUE),
+    rep(-Inf, 3)
+  )
+})
+
 test_that("pgeneralized_pareto distinguishes small nonzero shape", {
   k <- 9e-16
   q <- 1e7
diff --git a/tests/testthat/test-log_sum_exp.R b/tests/testthat/test-log_sum_exp.R
index 3e343938..55621198 100644
--- a/tests/testthat/test-log_sum_exp.R
+++ b/tests/testthat/test-log_sum_exp.R
@@ -9,6 +9,15 @@ test_that("log_sum_exp of -Inf is -Inf", {
   expect_equal(log_sum_exp(x), -Inf)
 })
 
+test_that("log_sum_exp ignores -Inf terms among finite ones", {
+  x <- c(-Inf, log(2), -Inf, log(3))
+  expect_equal(log_sum_exp(x), log(5))
+  expect_equal(log_sum_exp(x), log_sum_exp(c(log(2), log(3))))
+
+  # a single finite term survives any number of zero-probability terms
+  expect_equal(log_sum_exp(c(rep(-Inf, 1000), -700)), -700)
+})
+
 test_that("log_sum_exp of empty input is -Inf", {
   x <- numeric(0)
   expect_equal(log_sum_exp(x), -Inf)
diff --git a/tests/testthat/test-pareto_smooth.R b/tests/testthat/test-pareto_smooth.R
index f95bdb7e..14f64b62 100644
--- a/tests/testthat/test-pareto_smooth.R
+++ b/tests/testthat/test-pareto_smooth.R
@@ -269,6 +269,25 @@ test_that("exp_x_minus_exp_y is stable for nearby values", {
   expect_equal(exp_x_minus_exp_y(0, -1000), 1)
 })
 
+test_that("user-facing pareto functions reject -Inf draws", {
+  # -Inf is a valid zero weight for the internal helpers, but the user-facing
+  # functions take draws rather than weights and decline to fit
+  set.seed(1)
+  x <- c(-Inf, rnorm(999))
+  msg <- "Input contains infinite or NA values, is constant or has constant tail"
+
+  expect_warning(khat <- pareto_khat(x), msg)
+  expect_identical(khat, NA_real_)
+
+  expect_warning(diags <- pareto_diags(x), msg)
+  expect_named(diags, c("khat", "min_ss", "khat_threshold", "convergence_rate"))
+  expect_true(all(is.na(unlist(diags))))
+
+  # smoothing returns the draws unchanged rather than silently altering them
+  expect_warning(smoothed <- pareto_smooth(x), msg)
+  expect_identical(smoothed, x)
+})
+
 test_that("exp_x_minus_exp_y returns zero for equal infinite inputs", {
   expect_equal(exp_x_minus_exp_y(-Inf, -Inf), 0)
   expect_equal(exp_x_minus_exp_y(c(-Inf, 0, 1), c(-Inf, 0, 0)), c(0, 0, exp(1) - 1))

From 01205e8be618af8bd03f3147756204104f8a38db Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Tue, 15 Sep 2026 10:00:08 +0300
Subject: [PATCH 13/22] docs: note that subsetting weighted draws renormalizes
 the weights

---
 R/weight_draws.R    | 7 +++++++
 man/weight_draws.Rd | 8 ++++++++
 2 files changed, 15 insertions(+)

diff --git a/R/weight_draws.R b/R/weight_draws.R
index a7419d8a..3c2298df 100644
--- a/R/weight_draws.R
+++ b/R/weight_draws.R
@@ -6,6 +6,13 @@
 #' `.log_weight`). See [weights.draws()] for details how to extract weights from
 #' `draws` objects.
 #'
+#' Because the stored log-weights are unnormalized and are normalized only when
+#' they are extracted, subsetting a weighted `draws` object conditions on the
+#' retained draws: the weights of the draws that remain are renormalized to sum
+#' to one. Any operation that drops draws has this effect, including
+#' [subset_draws()], [thin_draws()] and `[` indexing. Weights are therefore
+#' comparable only within one subset, not across subsets of the same object.
+#'
 #' @template args-methods-x
 #' @param weights (numeric vector) A vector of weights of length `ndraws(x)`.
 #'   Weights will be internally stored on the log scale (in a variable called
diff --git a/man/weight_draws.Rd b/man/weight_draws.Rd
index d866d466..4a7bcd1e 100644
--- a/man/weight_draws.Rd
+++ b/man/weight_draws.Rd
@@ -49,6 +49,14 @@ are stored in the form of unnormalized log-weights (in a variable called
 \code{.log_weight}). See \code{\link[=weights.draws]{weights.draws()}} for details how to extract weights from
 \code{draws} objects.
 }
+\details{
+Because the stored log-weights are unnormalized and are normalized only when
+they are extracted, subsetting a weighted \code{draws} object conditions on the
+retained draws: the weights of the draws that remain are renormalized to sum
+to one. Any operation that drops draws has this effect, including
+\code{\link[=subset_draws]{subset_draws()}}, \code{\link[=thin_draws]{thin_draws()}} and \code{[} indexing. Weights are therefore
+comparable only within one subset, not across subsets of the same object.
+}
 \examples{
 x <- example_draws()
 

From 88b8bb6b3e3157272fe9ed7ddb1d74112d76c089 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Tue, 15 Sep 2026 10:07:29 +0300
Subject: [PATCH 14/22] fix: reject weights that carry no probability mass

---
 R/weight_draws.R                   | 12 ++++++++++++
 tests/testthat/test-weight_draws.R | 23 +++++++++++++++++++++++
 2 files changed, 35 insertions(+)

diff --git a/R/weight_draws.R b/R/weight_draws.R
index 3c2298df..2a43eee1 100644
--- a/R/weight_draws.R
+++ b/R/weight_draws.R
@@ -200,6 +200,18 @@ validate_weights <- function(weights, draws, log = FALSE) {
     }
     weights <- log(weights)
   }
+  if (!any(is.finite(weights))) {
+    # classed so that callers working variable by variable, such as pit(), can
+    # handle one degenerate column without aborting the whole call
+    stop(errorCondition(
+      paste0(
+        "All weights are zero, so no draw carries any probability mass. ",
+        "If the weights underflowed to zero, pass them on the log scale ",
+        "with `log = TRUE`."
+      ),
+      class = "posterior_degenerate_weights_error"
+    ))
+  }
   weights
 }
 
diff --git a/tests/testthat/test-weight_draws.R b/tests/testthat/test-weight_draws.R
index fb6e6cc3..66a25879 100644
--- a/tests/testthat/test-weight_draws.R
+++ b/tests/testthat/test-weight_draws.R
@@ -1,3 +1,26 @@
+test_that("weight_draws rejects weights with no probability mass", {
+  x <- example_draws()
+  n <- ndraws(x)
+  msg <- "All weights are zero"
+
+  expect_error(weight_draws(x, rep(0, n)), msg)
+  expect_error(weight_draws(x, rep(-Inf, n), log = TRUE), msg)
+  # underflow is the usual cause, so the message points at log = TRUE
+  expect_error(weight_draws(x, exp(rep(-800, n))), "log = TRUE", fixed = TRUE)
+
+  # the error is classed, so callers can handle it per variable
+  expect_error(weight_draws(x, rep(0, n)), class = "posterior_degenerate_weights_error")
+
+  # a single draw carrying all the mass is still valid
+  w <- c(1, rep(0, n - 1))
+  expect_equal(unname(weights(weight_draws(x, w))), w)
+
+  # all formats reject it
+  for (fmt in list(as_draws_array, as_draws_df, as_draws_list, as_draws_rvars)) {
+    expect_error(weight_draws(fmt(x), rep(0, n)), msg)
+  }
+})
+
 test_that("weight_draws works on draws_matrix", {
   x <- as_draws_matrix(example_draws())
   weights <- rexp(ndraws(x))

From 8cc5acf5afd5c2c77ecdc6dae4fff1e4b8d60bfb Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Tue, 15 Sep 2026 10:08:06 +0300
Subject: [PATCH 15/22] fix: error instead of returning NaN when normalizing
 zero weights

---
 R/weight_draws.R                   |  5 +++++
 tests/testthat/test-weight_draws.R | 24 ++++++++++++++++++++++++
 2 files changed, 29 insertions(+)

diff --git a/R/weight_draws.R b/R/weight_draws.R
index 2a43eee1..cc9bc4b7 100644
--- a/R/weight_draws.R
+++ b/R/weight_draws.R
@@ -176,6 +176,11 @@ weights.draws <- function(object, log = FALSE, normalize = TRUE, ...) {
   }
   out <- extract_variable(object, ".log_weight")
   if (normalize) {
+    # weight_draws() rejects this at construction, but dropping draws from an
+    # object that was valid can still remove every weighted draw
+    if (!any(is.finite(out))) {
+      stop_no_call("All draws have zero weight.")
+    }
     out <- out - log_sum_exp(out)
   }
   if (!log) {
diff --git a/tests/testthat/test-weight_draws.R b/tests/testthat/test-weight_draws.R
index 66a25879..52c4e2a9 100644
--- a/tests/testthat/test-weight_draws.R
+++ b/tests/testthat/test-weight_draws.R
@@ -88,6 +88,30 @@ test_that("weight_draws works on draws_rvars", {
 
 # conversion preserves weights --------------------------------------------
 
+test_that("weights() errors when subsetting has removed all the mass", {
+  x <- example_draws()
+  n <- ndraws(x)
+  # valid at construction: the first 10 draws carry all the mass
+  xw <- weight_draws(x, c(rep(1, 10), rep(0, n - 10)))
+  expect_silent(weights(xw))
+
+  # dropping those draws leaves nothing behind
+  # (subsetting by draw merges chains, hence the message)
+  sub <- suppressMessages(subset_draws(xw, draw = 50:100))
+  expect_error(weights(sub), "All draws have zero weight")
+  expect_error(weights(sub, log = TRUE), "All draws have zero weight")
+
+  # the unnormalized weights are still well defined, so they are still returned
+  expect_equal(unname(weights(sub, normalize = FALSE)), rep(0, ndraws(sub)))
+  expect_equal(unname(weights(sub, normalize = FALSE, log = TRUE)),
+               rep(-Inf, ndraws(sub)))
+
+  # a subset that keeps some mass is unaffected
+  kept <- suppressMessages(subset_draws(xw, draw = 1:5))
+  expect_silent(w <- weights(kept))
+  expect_equal(sum(w), 1)
+})
+
 test_that("conversion between formats preserves weights", {
   draws <- list(
     matrix = weight_draws(draws_matrix(x = 1:10), 1:10),

From 8ee7f1e280702839299743c9ced85a27a1c58eba Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Tue, 15 Sep 2026 10:10:03 +0300
Subject: [PATCH 16/22] fix: report degenerate PIT weight columns per variable
 instead of aborting

---
 R/pit.R                   | 57 ++++++++++++++++++++++++++++++++-------
 tests/testthat/test-pit.R | 37 +++++++++++++++++++++++++
 2 files changed, 85 insertions(+), 9 deletions(-)

diff --git a/R/pit.R b/R/pit.R
index 2c421664..0fe4f49a 100644
--- a/R/pit.R
+++ b/R/pit.R
@@ -79,12 +79,12 @@ pit.default <- function(x, y, weights = NULL, log = FALSE, ...) {
 pit.draws_matrix <- function(x, y, weights = NULL, log = FALSE, ...) {
   y <- validate_y(y, x)
   if (!is.null(weights)) {
-    weights <- sapply(seq_len(nvariables(x)), function(var_idx) {
-      validate_weights(weights[, var_idx], x[, var_idx], log)
-    })
-    weights <- normalize_log_weights(weights)
+    weights <- validate_pit_weights(weights, x, log)
   }
   pit <- vapply(seq_len(ncol(x)), function(j) {
+    if (!is.null(weights) && anyNA(weights[, j])) {
+      return(NA_real_)
+    }
     sel_min <- x[, j] < y[j]
     if (!any(sel_min)) {
       pit <- 0
@@ -112,7 +112,7 @@ pit.draws_matrix <- function(x, y, weights = NULL, log = FALSE, ...) {
     pit
   }, FUN.VALUE = 1.0)
 
-  if (any(pit > 1 + 1e-10)) {
+  if (any(pit > 1 + 1e-10, na.rm = TRUE)) {
     warning_no_call(
       paste(
         "Some PIT values larger than 1. ",
@@ -241,10 +241,7 @@ pareto_pit.draws_matrix <- function(x, y, weights = NULL, log = FALSE,
 
   # validate and normalize weights to log scale (same as pit.draws_matrix)
   if (!is.null(weights)) {
-    weights <- sapply(seq_len(nvariables(x)), function(var_idx) {
-      validate_weights(weights[, var_idx], x[, var_idx], log)
-    })
-    weights <- normalize_log_weights(weights)
+    weights <- validate_pit_weights(weights, x, log)
   }
 
   ndraws <- ndraws(x)
@@ -270,6 +267,9 @@ pareto_pit.draws_matrix <- function(x, y, weights = NULL, log = FALSE,
   }
 
   pit_values <- vapply(seq_len(ncol(x)), function(j) {
+    if (!is.null(weights) && anyNA(weights[, j])) {
+      return(NA_real_)
+    }
     draws <- x[, j]
 
     # --- raw PIT (same logic as pit.draws_matrix) ---
@@ -425,3 +425,42 @@ validate_y <- function(y, x = NULL) {
 normalize_log_weights <- function(log_weights) {
   apply(log_weights, 2, function(col) col - log_sum_exp(col))
 }
+
+#' Validate and normalize a per-variable weights matrix
+#'
+#' Columns whose weights carry no probability mass are reported and returned as
+#' `NA`, so that one degenerate variable does not take down the whole call. This
+#' matches [pareto_khat()], which warns and returns `NA` rather than erroring.
+#' Weights that are malformed for any other reason still raise an error.
+#'
+#' @noRd
+#' @param weights A matrix of weights, one column per variable in `x`.
+#' @param x A `draws_matrix`.
+#' @param log (logical) Are the weights already on the log scale?
+#' @return A matrix of normalized log weights, `NA` for degenerate columns.
+#'
+validate_pit_weights <- function(weights, x, log) {
+  vars <- variables(x)
+  cols <- lapply(seq_along(vars), function(var_idx) {
+    tryCatch(
+      validate_weights(weights[, var_idx], x[, var_idx], log),
+      posterior_degenerate_weights_error = function(cnd) NULL
+    )
+  })
+  degenerate <- vapply(cols, is.null, logical(1))
+  if (any(degenerate)) {
+    warning_no_call(
+      "All draws have zero weight for ",
+      paste0("'", vars[degenerate], "'", collapse = ", "),
+      ". Returning NA for ",
+      if (sum(degenerate) > 1) "those variables." else "that variable."
+    )
+    cols[degenerate] <- list(rep(NA_real_, ndraws(x)))
+  }
+  out <- do.call(cbind, cols)
+  # normalize only the columns that carry mass; log_sum_exp() has no NA handling
+  if (any(!degenerate)) {
+    out[, !degenerate] <- normalize_log_weights(out[, !degenerate, drop = FALSE])
+  }
+  out
+}
diff --git a/tests/testthat/test-pit.R b/tests/testthat/test-pit.R
index d7af97d9..403386a8 100644
--- a/tests/testthat/test-pit.R
+++ b/tests/testthat/test-pit.R
@@ -44,6 +44,43 @@ test_that("normalize_log_weights returns log-normalized columns", {
 })
 
 # tests for pit.default
+test_that("pit and pareto_pit report degenerate weight columns per variable", {
+  set.seed(1)
+  x <- draws_matrix(a = rnorm(400), b = rnorm(400))
+  n <- ndraws(x)
+  y <- c(0, 0)
+  good <- rep(log(1 / n), n)
+  msg <- "All draws have zero weight for 'a'"
+
+  # one bad column must not take down the other variable
+  w <- cbind(rep(-Inf, n), good)
+  expect_warning(out <- pit(x, y, weights = w, log = TRUE), msg)
+  expect_true(is.na(out[["a"]]))
+  expect_false(is.na(out[["b"]]))
+  expect_warning(out <- pareto_pit(x, y, weights = w, log = TRUE), msg)
+  expect_true(is.na(out[["a"]]))
+  expect_false(is.na(out[["b"]]))
+
+  # a single degenerate variable warns and returns NA rather than erroring,
+  # matching pareto_khat() rather than aborting
+  x1 <- subset_draws(x, variable = "a")
+  expect_warning(out <- pit(x1, 0, weights = matrix(-Inf, n, 1), log = TRUE), msg)
+  expect_true(is.na(out))
+
+  # every column degenerate: warn once naming both, all NA
+  expect_warning(
+    out <- pit(x, y, weights = matrix(-Inf, n, 2), log = TRUE),
+    "'a', 'b'"
+  )
+  expect_true(all(is.na(out)))
+
+  # ordinary zero weights take the same path
+  expect_warning(pit(x, y, weights = matrix(0, n, 2)), "All draws have zero weight")
+
+  # genuinely malformed weights still abort
+  expect_error(pit(x, y, weights = matrix(-1, n, 2)), "non-negative")
+})
+
 test_that("pit works without weights", {
   x <- matrix(c(1, 2, 3, 5), nrow = 2)
   y <- c(3, 4)

From 7472fb8067487f350bf0dac6cc32ef506fcce00c Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Tue, 15 Sep 2026 10:11:28 +0300
Subject: [PATCH 17/22] refactor: skip the Pareto tail fit when the tail
 carries no weight

---
 R/pit.R                          | 11 +++++++++--
 tests/testthat/test-pareto_pit.R | 33 ++++++++++++++++++++++++++++++++
 2 files changed, 42 insertions(+), 2 deletions(-)

diff --git a/R/pit.R b/R/pit.R
index 0fe4f49a..fe09b2dc 100644
--- a/R/pit.R
+++ b/R/pit.R
@@ -312,9 +312,13 @@ pareto_pit.draws_matrix <- function(x, y, weights = NULL, log = FALSE,
     }
 
     # --- right tail ---
+    # A tail carrying no probability mass cannot inform the fit: the raw PIT is
+    # already the whole answer, and gpdfit() would be handed all-zero weights.
+    right_tail_empty <- !is.null(log_wt_sorted) && tail_proportion == 0
+
     right_replaced <- FALSE
     right_tail <- sorted[tail_ids]
-    if (!is_constant(right_tail)) {
+    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
@@ -344,8 +348,11 @@ pareto_pit.draws_matrix <- function(x, y, weights = NULL, log = FALSE,
       left_sorted <- left_ord$x
       log_wt_left_sorted <- if (!is.null(weights)) weights[left_ord$ix, j] else NULL
 
+      left_tail_empty <- !is.null(log_wt_left_sorted) &&
+        log_sum_exp(log_wt_left_sorted[tail_ids]) == -Inf
+
       left_tail <- left_sorted[tail_ids]
-      if (!is_constant(left_tail)) {
+      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
diff --git a/tests/testthat/test-pareto_pit.R b/tests/testthat/test-pareto_pit.R
index 2293e8f7..9a7e9d0f 100644
--- a/tests/testthat/test-pareto_pit.R
+++ b/tests/testthat/test-pareto_pit.R
@@ -25,6 +25,39 @@ test_that("pareto_pit bulk values match pit", {
   expect_equal(unname(refined), unname(raw))
 })
 
+test_that("pareto_pit keeps the raw PIT when the fitted tail has no mass", {
+  set.seed(1)
+  x <- draws_matrix(a = rnorm(400), b = rnorm(400))
+  n <- ndraws(x)
+  # weights built from each variable's own ordering
+  mass_on <- function(idx) sapply(seq_len(nvariables(x)), function(j) {
+    out <- rep(-Inf, n)
+    out[order(as.numeric(x[, j]))[idx]] <- log(1 / length(idx))
+    out
+  })
+  # y strictly between two draws, so the randomized-tie path cannot fire
+  y_between <- function(rank) vapply(seq_len(nvariables(x)), function(j) {
+    sorted <- sort(as.numeric(x[, j]))
+    mean(sorted[c(rank, rank + 1)])
+  }, numeric(1))
+
+  # mass in the middle: both fitted tails are empty, so the GPD can add
+  # nothing and the refined PIT must equal the raw weighted PIT
+  w <- mass_on(190:209)
+  y <- y_between(200)
+  expect_equal(unname(pareto_pit(x, y, weights = w, log = TRUE)), c(0.55, 0.55))
+  expect_equal(pareto_pit(x, y, weights = w, log = TRUE),
+               pit(x, y, weights = w, log = TRUE))
+
+  # all mass below y: the raw PIT is 1, up to the min_tail_prob clamp that
+  # pareto_pit applies to every result
+  w <- mass_on(1:20)
+  y <- y_between(300)
+  expect_equal(unname(pit(x, y, weights = w, log = TRUE)), c(1, 1))
+  expect_equal(unname(pareto_pit(x, y, weights = w, log = TRUE)),
+               rep(1 - 1 / n / 1e4, 2))
+})
+
 test_that("pareto_pit differs from pit in tails", {
   set.seed(42)
   ndraws <- 1000

From f242a2c7e0def2d66cc0f480518a29ee46f53fde Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Tue, 15 Sep 2026 10:17:06 +0300
Subject: [PATCH 18/22] test: cover mixed finite and -Inf log weights in
 weight_draws

---
 tests/testthat/test-weight_draws.R | 36 ++++++++++++++++++++++++++++++
 1 file changed, 36 insertions(+)

diff --git a/tests/testthat/test-weight_draws.R b/tests/testthat/test-weight_draws.R
index 52c4e2a9..3e0793ef 100644
--- a/tests/testthat/test-weight_draws.R
+++ b/tests/testthat/test-weight_draws.R
@@ -88,6 +88,42 @@ test_that("weight_draws works on draws_rvars", {
 
 # conversion preserves weights --------------------------------------------
 
+test_that("weight_draws handles a mix of finite and -Inf log weights", {
+  x <- example_draws()
+  n <- ndraws(x)
+  zero <- c(rep(TRUE, 3), rep(FALSE, n - 3))
+  log_wts <- log(seq_len(n))
+  log_wts[zero] <- -Inf
+  wts <- exp(log_wts)
+
+  for (fmt in list(as_draws_matrix, as_draws_array, as_draws_df,
+                   as_draws_list, as_draws_rvars)) {
+    xf <- fmt(x)
+    from_log <- weight_draws(xf, log_wts, log = TRUE)
+    from_ordinary <- weight_draws(xf, wts)
+
+    # -Inf log weights and ordinary zeros describe the same distribution
+    expect_equal(weights(from_log), weights(from_ordinary))
+    expect_equal(weights(from_log, log = TRUE),
+                 weights(from_ordinary, log = TRUE))
+
+    w <- weights(from_log)
+    expect_equal(unname(w[zero]), rep(0, sum(zero)))
+    expect_true(all(w[!zero] > 0))
+    expect_equal(sum(w), 1)
+
+    # the zero-weight draws stay at -Inf on the log scale rather than
+    # underflowing to something finite
+    expect_equal(unname(weights(from_log, log = TRUE)[zero]),
+                 rep(-Inf, sum(zero)))
+
+    # unnormalized weights round-trip back to the input
+    expect_equal(unname(weights(from_log, normalize = FALSE)), wts)
+    expect_equal(unname(weights(from_log, normalize = FALSE, log = TRUE)),
+                 log_wts)
+  }
+})
+
 test_that("weights() errors when subsetting has removed all the mass", {
   x <- example_draws()
   n <- ndraws(x)

From b6ef3aaf53aa11ab33e66064d69ddd1d1ddeef15 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Tue, 15 Sep 2026 10:17:21 +0300
Subject: [PATCH 19/22] test: cover mixed finite and -Inf log weights in pit

---
 tests/testthat/test-pit.R | 35 +++++++++++++++++++++++++++++++++++
 1 file changed, 35 insertions(+)

diff --git a/tests/testthat/test-pit.R b/tests/testthat/test-pit.R
index 403386a8..6b620e05 100644
--- a/tests/testthat/test-pit.R
+++ b/tests/testthat/test-pit.R
@@ -44,6 +44,41 @@ test_that("normalize_log_weights returns log-normalized columns", {
 })
 
 # tests for pit.default
+test_that("pit treats -Inf log weights as ordinary zero weights", {
+  set.seed(1)
+  x <- draws_matrix(a = rnorm(400), b = rnorm(400))
+  n <- ndraws(x)
+  # y strictly between two draws, so the randomized-tie path cannot fire
+  y <- vapply(seq_len(nvariables(x)), function(j) {
+    sorted <- sort(as.numeric(x[, j]))
+    mean(sorted[c(200, 201)])
+  }, numeric(1))
+
+  # uniform mass on a subset, expressed both ways
+  keep <- seq(1, n, by = 3)
+  log_w <- matrix(-Inf, n, 2)
+  log_w[keep, ] <- log(1 / length(keep))
+  w <- exp(log_w)
+
+  expect_equal(pit(x, y, weights = log_w, log = TRUE), pit(x, y, weights = w))
+
+  # zero-weight draws must not influence the result at all: the answer is the
+  # unweighted PIT of the retained draws alone
+  expect_equal(
+    unname(pit(x, y, weights = log_w, log = TRUE)),
+    unname(pit(x[keep, ], y))
+  )
+
+  # unequal weights on the retained draws, still with zeros elsewhere
+  log_w2 <- matrix(-Inf, n, 2)
+  log_w2[keep, ] <- log(seq_along(keep))
+  expect_equal(pit(x, y, weights = log_w2, log = TRUE),
+               pit(x, y, weights = exp(log_w2)))
+  expect_equal(unname(pit(x, y, weights = log_w2, log = TRUE)),
+               unname(pit(x[keep, ], y, weights = cbind(seq_along(keep),
+                                                        seq_along(keep)))))
+})
+
 test_that("pit and pareto_pit report degenerate weight columns per variable", {
   set.seed(1)
   x <- draws_matrix(a = rnorm(400), b = rnorm(400))

From 0c76a4c9fa8b54cfc73fbc9a6200dcbd4289c1c3 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Tue, 15 Sep 2026 10:17:50 +0300
Subject: [PATCH 20/22] test: cover mixed finite and -Inf log weights in
 pareto_pit

---
 tests/testthat/test-pareto_pit.R | 35 ++++++++++++++++++++++++++++++++
 1 file changed, 35 insertions(+)

diff --git a/tests/testthat/test-pareto_pit.R b/tests/testthat/test-pareto_pit.R
index 9a7e9d0f..c1a3fdc2 100644
--- a/tests/testthat/test-pareto_pit.R
+++ b/tests/testthat/test-pareto_pit.R
@@ -25,6 +25,41 @@ test_that("pareto_pit bulk values match pit", {
   expect_equal(unname(refined), unname(raw))
 })
 
+test_that("pareto_pit treats -Inf log weights as ordinary zero weights", {
+  set.seed(1)
+  x <- draws_matrix(a = rnorm(400), b = rnorm(400))
+  n <- ndraws(x)
+  # y strictly between two draws, so the randomized-tie path cannot fire
+  y <- vapply(seq_len(nvariables(x)), function(j) {
+    sorted <- sort(as.numeric(x[, j]))
+    mean(sorted[c(380, 381)])
+  }, numeric(1))
+
+  # zeros spread through the sample, leaving the tails populated
+  keep <- seq(1, n, by = 3)
+  log_w <- matrix(-Inf, n, 2)
+  log_w[keep, ] <- log(1 / length(keep))
+  expect_equal(pareto_pit(x, y, weights = log_w, log = TRUE),
+               pareto_pit(x, y, weights = exp(log_w)))
+
+  # unequal weights, still with zeros interleaved
+  log_w2 <- matrix(-Inf, n, 2)
+  log_w2[keep, ] <- log(seq_along(keep))
+  expect_equal(pareto_pit(x, y, weights = log_w2, log = TRUE),
+               pareto_pit(x, y, weights = exp(log_w2)))
+
+  # a tail refinement actually happened, so the equivalence is not vacuous
+  expect_false(isTRUE(all.equal(
+    unname(pareto_pit(x, y, weights = log_w, log = TRUE)),
+    unname(pit(x, y, weights = log_w, log = TRUE))
+  )))
+
+  # unlike pit(), dropping the zero-weight draws is NOT equivalent here,
+  # because the number of tail draws is derived from ndraws
+  expect_equal(ps_tail_length(n, 1), 60)
+  expect_equal(ps_tail_length(length(keep), 1), 26)
+})
+
 test_that("pareto_pit keeps the raw PIT when the fitted tail has no mass", {
   set.seed(1)
   x <- draws_matrix(a = rnorm(400), b = rnorm(400))

From 60d0d9b3b92f0fc4d0b31253144ec960225ce94e Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Tue, 15 Sep 2026 10:33:55 +0300
Subject: [PATCH 21/22] fix: shift log_sum_exp by the maximum instead of
 clamping it at zero

---
 R/misc.R                           |  8 ++++++--
 tests/testthat/test-log_sum_exp.R  | 14 ++++++++++++++
 tests/testthat/test-weight_draws.R | 11 +++++++++++
 3 files changed, 31 insertions(+), 2 deletions(-)

diff --git a/R/misc.R b/R/misc.R
index 202fc9b8..db2ccf37 100644
--- a/R/misc.R
+++ b/R/misc.R
@@ -221,9 +221,13 @@ escape_all <- function(x) {
 
 # numerically stable version of log(sum(exp(x)))
 log_sum_exp <- function(x) {
-  max <- max(as.numeric(x), warnings = FALSE)
+  x <- as.numeric(x)
+  if (length(x) == 0) {
+    return(-Inf)
+  }
+  max <- max(x)
   if (max == -Inf) {
-    res <- 0
+    res <- -Inf
   } else if (max == Inf) {
     res <- Inf
   } else {
diff --git a/tests/testthat/test-log_sum_exp.R b/tests/testthat/test-log_sum_exp.R
index 55621198..c10e1127 100644
--- a/tests/testthat/test-log_sum_exp.R
+++ b/tests/testthat/test-log_sum_exp.R
@@ -23,6 +23,20 @@ test_that("log_sum_exp of empty input is -Inf", {
   expect_equal(log_sum_exp(x), -Inf)
 })
 
+test_that("log_sum_exp shifts by the maximum", {
+  # exp(x - 0) underflows for very negative x, so the shift has to be max(x)
+  expect_equal(log_sum_exp(c(-800, -800)), -800 + log(2))
+  expect_equal(log_sum_exp(rep(-1000, 5)), -1000 + log(5))
+  expect_equal(log_sum_exp(c(-800, -Inf)), -800)
+  expect_equal(log_sum_exp(c(-745, -745)), -745 + log(2))
+
+  # the result must shift with the input rather than depend on its scale
+  x <- c(-1, -2, -3)
+  for (shift in c(0, 500, -500, -1000)) {
+    expect_equal(log_sum_exp(x + shift), log_sum_exp(x) + shift)
+  }
+})
+
 test_that("log_sum_exp works", {
   set.seed(1)
   x <- log(runif(10))
diff --git a/tests/testthat/test-weight_draws.R b/tests/testthat/test-weight_draws.R
index 3e0793ef..8640db29 100644
--- a/tests/testthat/test-weight_draws.R
+++ b/tests/testthat/test-weight_draws.R
@@ -124,6 +124,17 @@ test_that("weight_draws handles a mix of finite and -Inf log weights", {
   }
 })
 
+test_that("weights() normalizes very negative unnormalized log weights", {
+  x <- example_draws()
+  n <- ndraws(x)
+  # a log likelihood summed over a few hundred observations lands here, and
+  # every value is finite, so nothing upstream rejects it
+  xw <- weight_draws(x, rep(-1000, n), log = TRUE)
+
+  expect_equal(unname(weights(xw)), rep(1 / n, n))
+  expect_equal(unname(weights(xw, log = TRUE)), rep(-log(n), n))
+})
+
 test_that("weights() errors when subsetting has removed all the mass", {
   x <- example_draws()
   n <- ndraws(x)

From 700d18de419accc0f921a4bc0c584658c836f39a Mon Sep 17 00:00:00 2001
From: Florence Bockting <florence.bockting@aalto.fi>
Date: Fri, 18 Sep 2026 17:41:59 +0300
Subject: [PATCH 22/22] refactor: add misc-math to outsource helpers for
 numerical stability

---
 R/gpd.R                                       |  4 +-
 R/misc-math.R                                 | 41 +++++++++++++++++++
 R/misc.R                                      | 32 ---------------
 R/pit.R                                       |  2 +-
 R/weight_draws.R                              |  2 +-
 .../{test-log_sum_exp.R => test-misc-math.R}  | 32 +++++++++++++++
 6 files changed, 77 insertions(+), 36 deletions(-)
 create mode 100644 R/misc-math.R
 rename tests/testthat/{test-log_sum_exp.R => test-misc-math.R} (59%)

diff --git a/R/gpd.R b/R/gpd.R
index 57eeaa5b..ee0ea9be 100644
--- a/R/gpd.R
+++ b/R/gpd.R
@@ -30,7 +30,7 @@ qgeneralized_pareto <- function(p, mu = 0, sigma = 1, k = 0, lower.tail = TRUE,
     p[invalid] <- NaN
   }
   if (log.p) {
-    log_survival <- if (lower.tail) log(-expm1(p)) else p
+    log_survival <- if (lower.tail) log1m_exp(p) else p
   } else {
     log_survival <- if (lower.tail) log1p(-p) else log(p)
   }
@@ -72,7 +72,7 @@ pgeneralized_pareto <- function(q, mu = 0, sigma = 1, k = 0, lower.tail = TRUE,
   }
   log_survival <- pmin(log_survival, 0)
   if (lower.tail) {
-    p <- if (log.p) log(-expm1(log_survival)) else -expm1(log_survival)
+    p <- if (log.p) log1m_exp(log_survival) else -expm1(log_survival)
   } else {
     p <- if (log.p) log_survival else exp(log_survival)
   }
diff --git a/R/misc-math.R b/R/misc-math.R
new file mode 100644
index 00000000..554c5214
--- /dev/null
+++ b/R/misc-math.R
@@ -0,0 +1,41 @@
+# numerically stable version of log(sum(exp(x)))
+log_sum_exp <- function(x) {
+  x <- as.numeric(x)
+  if (length(x) == 0) {
+    return(-Inf)
+  }
+  max <- max(x)
+  if (max == -Inf) {
+    res <- -Inf
+  } else if (max == Inf) {
+    res <- Inf
+  } else {
+    sum <- sum(exp(x - max))
+    res <- max + log(sum)
+  }
+  res
+}
+
+# numerically stable version of exp(x) - exp(y)
+exp_x_minus_exp_y <- function(x, y) {
+  out <- -exp(x) * expm1(y - x)
+  # equal infinite inputs give 0 * NaN above; which() drops the NA comparisons
+  # that NA or NaN inputs would produce
+  out[which(x == y)] <- 0
+  out
+}
+
+# numerically stable version of log2(1 - x)
+log2_one_minus <- function(x) {
+  log1p(-x) / log(2)
+}
+
+# numerically stable version of log(1 - exp(x)) for x <= 0
+log1m_exp <- function(x) {
+  log(-expm1(x))
+}
+
+# numerically stable version of log(exp(x) / sum(exp(x)))
+log_normalize <- function(x) {
+  x - log_sum_exp(x)
+}
diff --git a/R/misc.R b/R/misc.R
index db2ccf37..e96ae29f 100644
--- a/R/misc.R
+++ b/R/misc.R
@@ -219,38 +219,6 @@ escape_all <- function(x) {
   gsub(specials, "\\\\\\1", x)
 }
 
-# numerically stable version of log(sum(exp(x)))
-log_sum_exp <- function(x) {
-  x <- as.numeric(x)
-  if (length(x) == 0) {
-    return(-Inf)
-  }
-  max <- max(x)
-  if (max == -Inf) {
-    res <- -Inf
-  } else if (max == Inf) {
-    res <- Inf
-  } else {
-    sum <- sum(exp(x - max))
-    res <- max + log(sum)
-  }
-  res
-}
-
-# numerically stable version of exp(x) - exp(y)
-exp_x_minus_exp_y <- function(x, y) {
-  out <- -exp(x) * expm1(y - x)
-  # equal infinite inputs give 0 * NaN above; which() drops the NA comparisons
-  # that NA or NaN inputs would produce
-  out[which(x == y)] <- 0
-  out
-}
-
-# numerically stable version of log2(1 - x)
-log2_one_minus <- function(x) {
-  log1p(-x) / log(2)
-}
-
 # simple version of destructuring assignment
 `%<-%` <- function(vars, values, envir = parent.frame()) {
   vars <- as.character(substitute(vars)[-1])
diff --git a/R/pit.R b/R/pit.R
index fe09b2dc..d3789c47 100644
--- a/R/pit.R
+++ b/R/pit.R
@@ -430,7 +430,7 @@ validate_y <- function(y, x = NULL) {
 }
 
 normalize_log_weights <- function(log_weights) {
-  apply(log_weights, 2, function(col) col - log_sum_exp(col))
+  apply(log_weights, 2, log_normalize)
 }
 
 #' Validate and normalize a per-variable weights matrix
diff --git a/R/weight_draws.R b/R/weight_draws.R
index cc9bc4b7..fee48807 100644
--- a/R/weight_draws.R
+++ b/R/weight_draws.R
@@ -181,7 +181,7 @@ weights.draws <- function(object, log = FALSE, normalize = TRUE, ...) {
     if (!any(is.finite(out))) {
       stop_no_call("All draws have zero weight.")
     }
-    out <- out - log_sum_exp(out)
+    out <- log_normalize(out)
   }
   if (!log) {
     out <- exp(out)
diff --git a/tests/testthat/test-log_sum_exp.R b/tests/testthat/test-misc-math.R
similarity index 59%
rename from tests/testthat/test-log_sum_exp.R
rename to tests/testthat/test-misc-math.R
index c10e1127..cf9e00c8 100644
--- a/tests/testthat/test-log_sum_exp.R
+++ b/tests/testthat/test-misc-math.R
@@ -42,3 +42,35 @@ test_that("log_sum_exp works", {
   x <- log(runif(10))
   expect_equal(log_sum_exp(x), log(sum(exp(x))))
 })
+
+test_that("log1m_exp agrees with log(1 - exp(x))", {
+  x <- c(-5, -1, -0.1)
+  expect_equal(log1m_exp(x), log(1 - exp(x)))
+})
+
+test_that("log1m_exp is stable near zero", {
+  # 1 - exp(-1e-20) rounds to 0, so the naive version gives -Inf
+  expect_equal(log1m_exp(-1e-20), log(1e-20))
+})
+
+test_that("log1m_exp handles the boundaries", {
+  expect_equal(log1m_exp(0), -Inf)
+  expect_equal(log1m_exp(-Inf), 0)
+})
+
+test_that("log_normalize gives weights that sum to 1", {
+  set.seed(1)
+  x <- log(runif(10))
+  expect_equal(sum(exp(log_normalize(x))), 1)
+  expect_equal(log_normalize(x), log(exp(x) / sum(exp(x))))
+})
+
+test_that("log_normalize is stable for large inputs", {
+  # exp(1000) overflows, so the naive version gives NaN
+  expect_equal(log_normalize(c(1000, 1000)), rep(log(0.5), 2))
+  expect_equal(log_normalize(c(-1000, -1000)), rep(log(0.5), 2))
+})
+
+test_that("log_normalize keeps zero weights at -Inf", {
+  expect_equal(log_normalize(c(-Inf, 0, 0)), c(-Inf, log(0.5), log(0.5)))
+})