Use a relative tolerance in is_constant

[avehtari]

This PR was assisted by Claude

Summary

is_constant() uses an absolute tolerance: abs(max(x) - min(x)) < .Machine$double.eps. Draws whose values all have a small magnitude are then treated as constant, and the diagnostics that call should_return_NA() or is_constant() return NA. The result depends on the units of the variable.

This PR makes the tolerance relative to the largest absolute value in x. The diagnostics then give the same result for x and s * x.

Problem

library(posterior)
set.seed(1234)
x <- matrix(exp(rnorm(400)), ncol = 4)
ess_sd(x)                 # 437.6
ess_sd(1e-10 * x)         # master: NA, this PR: 437.6
rhat_basic(1e-20 * x)     # master: NA, this PR: 1.0031 (same as rhat_basic(x))
ess_mean(1e-20 * x)       # master: NA, this PR: 415.6 (same as ess_mean(x))
pareto_khat(1e-20 * x)    # master: NA with a warning, this PR: 0.2066 (same as pareto_khat(x))

ess_sd() works on squared deviations, so it already fails for draws with a standard deviation of about 1e-8 or less. Examples are parameters in SI units and natural-scale importance weights.

The affected functions call should_return_NA() or is_constant() on the raw draws: rhat_basic(), ess_basic(), ess_mean(), ess_sd(), ess_quantile() and the functions that use it (ess_tail(), ess_median()), the mcse_*() functions, rhat_nested(), pareto_khat(), pareto_diags(), pareto_smooth(), ps_tail() and pareto_pit(). Thus summarise_draws() with the default measures shows ess_tail = NA for such a variable. rhat() and ess_bulk() are not affected, because they call the check on rank-normalized values. With the same x and s = 1e-100, master gives NA for ess_tail(), ess_median() and mcse_mean(), and this PR gives 385.2, 330.3 and 1e-100 times 0.109, the same as for x.

Change

• x is constant if max(x) - min(x) <= tol * max(abs(x)). The result does not depend on the scale of x.
• x with all values equal is constant, also when all values are 0 or all are -Inf. On master, all -Inf gave NA.
• x with an infinite spread, for example c(-Inf, 2), is not constant, the same as on master.
• x with NA gives NA, the same as on master.
• The documentation of tol now says that it is relative.

Other behavior change: for values of magnitude above 1, values that differ by about 1 ulp are now constant. For example, is_constant(c(1, 1 + .Machine$double.eps)) is TRUE (on master FALSE). For magnitudes of 2 or more, master treated values that differ by 1 ulp as different, because 1 ulp is then larger than .Machine$double.eps.

Tests

• tests/testthat/test-convergence.R:
  • is_constant uses a tolerance relative to the magnitude of x: zeros, small constant values, -Inf, 1-ulp differences, small values that differ, infinite spread, NA.
  • convergence diagnostics do not depend on the scale of the draws: rhat_basic(), ess_basic(), ess_mean(), ess_quantile() and rhat_nested() for s = 1e-10 and 1e-100, and ess_sd() for s = 1e-10 and 1e-50. For s = 1e-100, the autocovariance of the squared deviations is of order s^4 and underflows.
• tests/testthat/test-pareto_smooth.R: pareto diagnostics and smoothing do not depend on the scale of the draws: pareto_khat(), pareto_diags(), pareto_smooth() and ps_tail() for s = 1e-20 and 1e-200.
• tests/testthat/test-pareto_pit.R: pareto_pit does not depend on the scale of the draws, for s = 1e-20 and 1e-200.

On unpatched master, the new tests give 22 failed expectations: 12 in test-convergence.R, 8 in test-pareto_smooth.R and 2 in test-pareto_pit.R. With this PR, all tests pass.

Effect on other input

On master and with this PR, these results are identical:

• summarise_draws(example_draws()), with the default measures and with default_convergence_measures().
• summarise_draws() with default_convergence_measures() for draws 1e10 + rnorm() and for Poisson draws (integer values with ties).
• pareto_khat(example_draws()) and pareto_pit(example_draws(), rep(1, 10)).

Related

#468 (removes the offset at a tied Pareto tail cutoff). That PR makes the Pareto smoothing independent of the scale of the tail excesses. This PR does the same for the checks for constant input.

Tested versions

• posterior: this branch, based on upstream master 6dd1e12383.
• R 4.6.1, testthat 3.3.2, Linux x86_64.
• devtools::test(): 2253 tests, 0 failures, 0 errors, 0 warnings, 4 skips (gbm is not installed).
• man/is_constant.Rd is generated with roxygen2 8.0.0. This PR does not include the changes that roxygen2 8.0.0 makes to other Rd files and to DESCRIPTION, because the package uses RoxygenNote 7.3.3.

Copyright and Licensing

By submitting this pull request, the copyright holder is agreeing to
license the submitted work under the following licenses:

• Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
• Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

[github-actions]

This is how benchmark results would change (along with a 95% confidence interval in relative change) if 641e907 is merged into master:

• ✔️as_draws_array: 150ms -> 150ms [-1.21%, +0.59%]
• ✔️as_draws_df: 73.1ms -> 73.9ms [-0.42%, +2.59%]
• ✔️as_draws_list: 173ms -> 174ms [-0.52%, +1.62%]
• ✔️as_draws_matrix: 25.4ms -> 25.3ms [-1.29%, +0.9%]
• ✔️as_draws_rvars: 117ms -> 118ms [-0.48%, +2.28%]
• ✔️summarise_draws_100_variables: 749ms -> 749ms [-0.24%, +0.26%]
• ✔️summarise_draws_10_variables: 82.9ms -> 82.6ms [-0.76%, +0.13%]
  Further explanation regarding interpretation and methodology can be found in the documentation.

[florence-bockting]

Thank you. Looks good to me. I just added a bullet point in the NEWS.md