Report the Pareto k diagnostic for sis() and tis()

NEWS.md

@@ -1,5 +1,8 @@
 # loo (development version)
 
+* `sis()` and `tis()` now compute the Pareto k diagnostic in the same way as
+`psis()`, instead of setting all values to 0, and give the same warnings. The
+weights do not change. by @avehtari in #411
 * Improve numerical stability in `loo()`, `psis()`, model weighting, subsampling, 
 and moment matching in #395 
 * Fix `loo_compare()` when used with subsampling: compute model comparison by 

R/importance_sampling.R

@@ -196,10 +196,11 @@ do_importance_sampling <- function(log_ratios, r_eff, cores, method) {
   S <- nrow(log_ratios)
   k_threshold <- ps_khat_threshold(S)
   tail_len <- n_pareto(r_eff, S)
+  # all methods report the Pareto k diagnostic, which needs a long enough tail
+  throw_tail_length_warnings(tail_len)
 
   if (method == "psis") {
     is_fun <- do_psis_i
-    throw_tail_length_warnings(tail_len)
   } else if (method == "tis") {
     is_fun <- do_tis_i
   } else if (method == "sis") {

R/print.R

@@ -164,7 +164,7 @@ print_dims.importance_sampling_loo <- function(x, ...) {
   dims <- dim(x)
   cat(
     "Computed from", dims[1], "posterior draws and",
-    dims[2], "log-likelihood terms using", class(x)[1], ".\n"
+    dims[2], "log-likelihood terms using", paste0(class(x)[1], ".\n")
   )
 }
 

R/psis.R

@@ -305,6 +305,39 @@ enough_tail_samples <- function(tail_len, min_len = 5) {
   tail_len >= min_len
 }
 
+#' Pareto k diagnostic for a single vector of log ratios
+#'
+#' Fits the generalized Pareto distribution to the same tail as `do_psis_i()`,
+#' without smoothing, so that `sis()` and `tis()` report the same diagnostic as
+#' `psis()` for the same log ratios.
+#'
+#' @noRd
+#' @param log_ratios_i A vector of log importance ratios.
+#' @param tail_len_i An integer tail length.
+#' @return Scalar Pareto k estimate. `Inf` if the tail is too short or constant,
+#'   or if the fit fails, as in `do_psis_i()`.
+#'
+pareto_k_i <- function(log_ratios_i, tail_len_i) {
+  if (!enough_tail_samples(tail_len_i)) {
+    return(Inf)
+  }
+  S <- length(log_ratios_i)
+  lw_sorted <- sort.int(log_ratios_i - max(log_ratios_i))
+  tail_ids <- seq(S - tail_len_i + 1, S)
+  lw_tail <- lw_sorted[tail_ids]
+  if (abs(max(lw_tail) - min(lw_tail)) < .Machine$double.eps / 100) {
+    warning(
+      "Can't fit generalized Pareto distribution ",
+      "because all tail values are the same.",
+      call. = FALSE
+    )
+    return(Inf)
+  }
+  cutoff <- lw_sorted[min(tail_ids) - 1]
+  k <- posterior::gpdfit(exp_x_minus_exp_y(lw_tail, cutoff), sort_x = FALSE)$k
+  if (is.na(k)) Inf else k
+}
+
 
 #' Throw warnings about Pareto k estimates
 #'

R/sis.R

@@ -30,7 +30,9 @@
 #'   }
 #'  \item{`diagnostics`}{
 #'    A named list containing one vector:
-#'    * `pareto_k`: Not used in `sis`, all set to 0.
+#'    * `pareto_k`: Estimates of the shape parameter \eqn{k} of the
+#'    generalized Pareto distribution, computed from the log ratios as in
+#'    [psis()]. `sis()` uses them only as a diagnostic.
 #'    * `n_eff`: effective sample size estimates.
 #'  }
 #' }
@@ -45,7 +47,7 @@
 #'     If specified, the user's `r_eff` argument.
 #'   }
 #'   \item{`tail_len`}{
-#'     Not used for `sis`.
+#'     Vector of tail lengths used for the Pareto k diagnostic.
 #'   }
 #'   \item{`dims`}{
 #'     Integer vector of length 2 containing `S` (posterior sample size)
@@ -135,13 +137,16 @@ is.sis <- function(x) {
 #' @noRd
 #' @param log_ratios_i A vector of log importance ratios (for `loo()`, negative
 #'   log likelihoods).
+#' @param tail_len_i An integer tail length for the Pareto k diagnostic.
 #' @param ... Not used. Included to conform to PSIS API.
 #'
 #' @details Implementation standard importance sampling.
 #' @return A named list containing:
 #' * `lw`: vector of unnormalized log weights
-#' * `pareto_k`: scalar Pareto k estimate. For IS, this defaults to 0.
-do_sis_i <- function(log_ratios_i, ...) {
-  S <- length(log_ratios_i)
-  list(log_weights = log_ratios_i, pareto_k = 0)
+#' * `pareto_k`: scalar Pareto k estimate, the same as for `psis()`.
+do_sis_i <- function(log_ratios_i, tail_len_i, ...) {
+  list(
+    log_weights = log_ratios_i,
+    pareto_k = pareto_k_i(log_ratios_i, tail_len_i)
+  )
 }

R/tis.R

@@ -33,7 +33,9 @@
 #'   }
 #'  \item{`diagnostics`}{
 #'    A named list containing one vector:
-#'    * `pareto_k`: Not used in `tis`, all set to 0.
+#'    * `pareto_k`: Estimates of the shape parameter \eqn{k} of the
+#'    generalized Pareto distribution, computed from the log ratios before
+#'    truncation as in [psis()]. `tis()` uses them only as a diagnostic.
 #'    * `n_eff`: Effective sample size estimates.
 #'  }
 #' }
@@ -48,7 +50,7 @@
 #'     If specified, the user's `r_eff` argument.
 #'   }
 #'   \item{`tail_len`}{
-#'     Not used for `tis`.
+#'     Vector of tail lengths used for the Pareto k diagnostic.
 #'   }
 #'   \item{`dims`}{
 #'     Integer vector of length 2 containing `S` (posterior sample size)
@@ -139,6 +141,7 @@ is.tis <- function(x) {
 #' @noRd
 #' @param log_ratios_i A vector of log importance ratios (for `loo()`, negative
 #'   log likelihoods).
+#' @param tail_len_i An integer tail length for the Pareto k diagnostic.
 #' @param ... Not used. Included to conform to PSIS API.
 #'
 #' @details Implementation of Truncated importance sampling (TIS), a method for
@@ -148,13 +151,14 @@ is.tis <- function(x) {
 #'
 #' @return A named list containing:
 #' * `lw`: vector of unnormalized log weights
-#' * `pareto_k`: scalar Pareto k estimate. For 'tis', this defaults  to 0.
+#' * `pareto_k`: scalar Pareto k estimate of the log ratios before truncation,
+#'   the same as for `psis()`.
 #'
-do_tis_i <- function(log_ratios_i, ...) {
+do_tis_i <- function(log_ratios_i, tail_len_i, ...) {
   S <- length(log_ratios_i)
   log_Z <- logSumExp(log_ratios_i) - log(S) # Normalization term, c-hat in Ionides (2008) appendix
   log_cutpoint <- log_Z + 0.5 * log(S)
   lw_i <- pmin(log_ratios_i, log_cutpoint)
-  list(log_weights = lw_i, pareto_k = 0)
+  list(log_weights = lw_i, pareto_k = pareto_k_i(log_ratios_i, tail_len_i))
 }
 

tests/testthat/_snaps/tisis.md

@@ -17,7 +17,7 @@
       print(loo_tis)
     Output
 
-      Computed from 1000 posterior draws and 32 log-likelihood terms using tis_loo .
+      Computed from 1000 posterior draws and 32 log-likelihood terms using tis_loo.
 
                Estimate  SE
       elpd_loo    -83.6 4.3
@@ -31,7 +31,7 @@
       print(loo_sis)
     Output
 
-      Computed from 1000 posterior draws and 32 log-likelihood terms using sis_loo .
+      Computed from 1000 posterior draws and 32 log-likelihood terms using sis_loo.
 
                Estimate  SE
       elpd_loo    -83.6 4.3

tests/testthat/test_tisis.R

@@ -172,7 +172,12 @@ for (method_name in c("psis", "tis", "sis")) {
 
 test_that("explict test of values for 'sis' and 'tis'", {
   lw <- 1:16
-  expect_silent(tis_true <- tis(log_ratios = lw, r_eff = NA))
+  # too few draws for the Pareto k diagnostic, as in psis(); k is then Inf,
+  # which gives a second warning
+  suppressWarnings(expect_warning(
+    tis_true <- tis(log_ratios = lw, r_eff = NA),
+    "Not enough tail samples"
+  ))
   expect_equal(
     as.vector(weights(tis_true, log = TRUE, normalize = FALSE)),
     c(
@@ -196,7 +201,10 @@ test_that("explict test of values for 'sis' and 'tis'", {
       15.07238,
     tolerance = 0.001
   )
-  expect_silent(is_true <- sis(log_ratios = lw, r_eff = NA))
+  suppressWarnings(expect_warning(
+    is_true <- sis(log_ratios = lw, r_eff = NA),
+    "Not enough tail samples"
+  ))
   expect_equal(
     as.vector(weights(is_true, log = TRUE, normalize = FALSE)),
     lw,
@@ -216,7 +224,12 @@ test_that("explict test of values for 'sis' and 'tis'", {
     1.5599302
   )
 
-  expect_silent(tis_true <- tis(log_ratios = lw, r_eff = NA))
+  # too few draws for the Pareto k diagnostic, as in psis(); k is then Inf,
+  # which gives a second warning
+  suppressWarnings(expect_warning(
+    tis_true <- tis(log_ratios = lw, r_eff = NA),
+    "Not enough tail samples"
+  ))
   expect_equal(
     as.vector(weights(tis_true, log = TRUE, normalize = FALSE)),
     c(
@@ -234,7 +247,10 @@ test_that("explict test of values for 'sis' and 'tis'", {
       3.692668,
     tolerance = 0.001
   )
-  expect_silent(is_true <- sis(log_ratios = lw, r_eff = NA))
+  suppressWarnings(expect_warning(
+    is_true <- sis(log_ratios = lw, r_eff = NA),
+    "Not enough tail samples"
+  ))
   expect_equal(
     as.vector(weights(is_true, log = TRUE, normalize = FALSE)),
     lw,
@@ -257,3 +273,43 @@ test_that("tis_loo and sis_loo are returned", {
   expect_snapshot(print(loo_tis))
   expect_snapshot(print(loo_sis))
 })
+
+test_that("psis, sis and tis report the same Pareto k diagnostic", {
+  log_ratios <- -example_loglik_array()
+  r_eff <- relative_eff(exp(-log_ratios))
+  k_psis <- pareto_k_values(psis(log_ratios, r_eff = r_eff))
+  expect_true(all(is.finite(k_psis)))
+  expect_equal(pareto_k_values(sis(log_ratios, r_eff = r_eff)), k_psis)
+  expect_equal(pareto_k_values(tis(log_ratios, r_eff = r_eff)), k_psis)
+
+  # an observation with a thick tail gives a high k and a warning for all three
+  set.seed(227)
+  log_ratios_thick <- log_ratios
+  log_ratios_thick[, , 1] <- rt(length(log_ratios[, , 1]), df = 1)
+  for (method in list(psis, sis, tis)) {
+    expect_warning(
+      k <- pareto_k_values(method(log_ratios_thick, r_eff = r_eff)),
+      "Some Pareto k diagnostic values are too high"
+    )
+    expect_gt(k[1], 0.7)
+  }
+})
+
+test_that("psis, sis and tis report Inf for too few draws", {
+  # 3 iterations of 2 chains are too few for the tail of the diagnostic
+  log_ratios <- -example_loglik_array()[1:3, , , drop = FALSE]
+  r_eff <- relative_eff(exp(-log_ratios))
+  for (method in list(psis, sis, tis)) {
+    warnings <- character(0)
+    k <- withCallingHandlers(
+      pareto_k_values(method(log_ratios, r_eff = r_eff)),
+      warning = function(w) {
+        warnings <<- c(warnings, conditionMessage(w))
+        invokeRestart("muffleWarning")
+      }
+    )
+    expect_equal(k, rep(Inf, 32))
+    expect_match(warnings, "Not enough tail samples", all = FALSE)
+    expect_match(warnings, "Pareto k diagnostic values are too high", all = FALSE)
+  }
+})