Use a relative tolerance in is_constant

NEWS.md

@@ -11,7 +11,8 @@ selected.
 
 * Fix truncation operation for `uniformity_test()` (#462)
 * mid-p value adjusment for discrete `PRIT` test method  
-* Add additional unit tests for `uniformity_test()` 
+* Add additional unit tests for `uniformity_test()`
+* Use a tolerance relative to the magnitude of `x` in `is_constant()` (#469)
 
 # posterior 1.7.1
 

R/misc.R

@@ -184,16 +184,22 @@ isNA <- function(x) {
 
 #' Check if vector is constant
 #'
-#' check if a vector is constant, up to a defined tolerance.
+#' check if a vector is constant, up to a defined tolerance relative to the
+#' magnitude of its values.
 #'
 #' @param x (vector) vector to check if is constant
-#' @param tol (numeric) tolerance to consider two values equal.
+#' @param tol (numeric) relative tolerance to consider two values equal. `x`
+#'   is constant if `max(x) - min(x)` is at most `tol` times the largest
+#'   absolute value in `x`.
 #'   Default is `.Machine$double.eps`.
 #'
 #' @keywords internal
 #' @export
 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)))
 }
 
 # collapse variables via commas for pretty printing

tests/testthat/test-convergence.R

@@ -115,6 +115,42 @@ test_that("convergence diagnostics handle special cases correctly", {
   # expect_true(is.na(ess_basic(x)))
 })
 
+test_that("is_constant uses a tolerance relative to the magnitude of x", {
+  eps <- .Machine$double.eps
+  expect_true(is_constant(rep(0, 5)))
+  expect_true(is_constant(rep(-3e-300, 5)))
+  expect_true(is_constant(c(-Inf, -Inf)))
+  # values that differ by at most eps relative to their magnitude
+  expect_true(is_constant(c(1, 1 + eps)))
+  expect_true(is_constant(1e10 * c(1, 1 + eps)))
+  # small values that differ relative to their magnitude
+  expect_false(is_constant(c(0, 1e-300)))
+  expect_false(is_constant(1e-20 * c(1, 2)))
+  expect_false(is_constant(c(1, 1 + 4 * eps)))
+  expect_false(is_constant(c(-Inf, 2)))
+  expect_true(is.na(is_constant(c(1, NA))))
+})
+
+test_that("convergence diagnostics do not depend on the scale of the draws", {
+  set.seed(1234)
+  x <- matrix(exp(rnorm(400)), ncol = 4)
+  for (s in c(1e-10, 1e-100)) {
+    expect_equal(rhat_basic(s * x), rhat_basic(x))
+    expect_equal(ess_basic(s * x), ess_basic(x))
+    expect_equal(ess_mean(s * x), ess_mean(x))
+    expect_equal(ess_quantile(s * x, 0.9), ess_quantile(x, 0.9))
+    expect_equal(
+      rhat_nested(s * x, superchain_ids = c(1, 1, 2, 2)),
+      rhat_nested(x, superchain_ids = c(1, 1, 2, 2))
+    )
+  }
+  # ess_sd() works on squared deviations, and their autocovariance is of
+  # the order s^4, which underflows for s = 1e-100
+  for (s in c(1e-10, 1e-50)) {
+    expect_equal(ess_sd(s * x), ess_sd(x))
+  }
+})
+
 test_that("convergence diagnostics throw correct errors", {
   mu <- extract_variable_matrix(example_draws(), "mu")
   expect_error(ess_quantile(mu, probs = 1.2), "'probs' must contain values between 0 and 1")

tests/testthat/test-pareto_pit.R

@@ -446,3 +446,12 @@ test_that("pareto_pit.default passes weights through", {
   expect_length(result, 2)
   expect_true(all(result >= 0 & result <= 1))
 })
+
+test_that("pareto_pit does not depend on the scale of the draws", {
+  set.seed(42)
+  x <- matrix(exp(rnorm(4000)), ncol = 4)
+  y <- c(0.05, 1, 10, 30)
+  for (s in c(1e-20, 1e-200)) {
+    expect_equal(pareto_pit(s * x, s * y), pareto_pit(x, y))
+  }
+})

tests/testthat/test-pareto_smooth.R

@@ -364,3 +364,17 @@ test_that("ps_convergence_rate propagates missing shape parameters", {
     c(ps_convergence_rate(0.3, n), NA, ps_convergence_rate(0.8, n), NA, 1, 0)
   )
 })
+
+test_that("pareto diagnostics and smoothing do not depend on the scale of the draws", {
+  set.seed(1234)
+  x <- exp(rnorm(1000))
+  for (s in c(1e-20, 1e-200)) {
+    expect_equal(pareto_khat(s * x), pareto_khat(x))
+    expect_equal(pareto_diags(s * x), pareto_diags(x))
+    expect_equal(
+      pareto_smooth(s * x, return_k = FALSE, verbose = FALSE) / s,
+      pareto_smooth(x, return_k = FALSE, verbose = FALSE)
+    )
+    expect_equal(ps_tail(s * x, ndraws_tail = 50)$k, ps_tail(x, ndraws_tail = 50)$k)
+  }
+})