Use a relative tolerance in is_constant

NEWS.md

@@ -9,6 +9,8 @@ comparison-model minus reference-model by @florence-bockting in #391
 * Fix `crps()` and `scrps()` ignoring the `permutations` argument for numeric
 input by @thanasibakis in #398
 * Speed up the Bayesian bootstrap in `pseudobma_weights()` by @VisruthSK in #368
+* Use a tolerance relative to the magnitude of `x` in `is_constant()` by
+@avehtari in #407
 
 # loo 2.10.1
 

R/helpers.R

@@ -277,6 +277,11 @@ release_questions <- function() {
 }
 # nocov end
 
+# `x` is constant if `max(x) - min(x)` is at most `tol` times the largest
+# absolute value in `x`. Default is `.Machine$double.eps`.
 is_constant <- function(x, tol = .Machine$double.eps) {
-  abs(max(x) - min(x)) < tol
+  r <- range(x)
+  spread <- r[2] - r[1]
+  # the equality check covers all-zero and all-infinite x
+  r[1] == r[2] || (is.finite(spread) && spread <= tol * max(abs(r)))
 }

tests/testthat/test_0_helpers.R

@@ -91,3 +91,15 @@ test_that("loo_cores works", {
   expect_warning(expect_equal(loo_cores(10), 2), "deprecated")
   options(loo.cores = NULL)
 })
+
+test_that("is_constant uses a tolerance relative to the magnitude of x", {
+  expect_true(is_constant(rep(0, 5)))
+  expect_true(is_constant(rep(1e-300, 5)))
+  expect_true(is_constant(rep(-Inf, 5)))
+  expect_true(is_constant(c(1, 1 + .Machine$double.eps)))
+  expect_false(is_constant(c(1e-20, 2e-20)))
+  expect_false(is_constant(c(1e-300, 2e-300)))
+  expect_false(is_constant(c(-Inf, 2)))
+  expect_false(is_constant(c(1, 2)))
+  expect_true(is.na(is_constant(c(1, NA))))
+})

tests/testthat/test_E_loo.R

@@ -245,3 +245,15 @@ test_that("weighted variance works", {
   x <- 1e8 + c(-1, 0, 1)
   expect_equal(.wvar(x, rep(1 / 3, 3)), 1)
 })
+
+test_that("E_loo pareto_k does not depend on the scale of x", {
+  # posterior::pareto_khat() of h * r also needs the relative tolerance
+  # with an older posterior it returns NA here
+  skip_if(posterior:::is_constant(c(1e-20, 2e-20)),
+          "posterior::is_constant() has an absolute tolerance")
+  k <- E_loo(x, psis_mat, type = "mean", log_ratios = log_rats)$pareto_k
+  for (s in c(1e-20, 1e-200)) {
+    k_s <- E_loo(s * x, psis_mat, type = "mean", log_ratios = log_rats)$pareto_k
+    expect_equal(k_s, k, tolerance = 1e-10)
+  }
+})