diff --git a/NEWS.md b/NEWS.md index 877c1afe..336166bb 100644 --- a/NEWS.md +++ b/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 diff --git a/R/importance_sampling.R b/R/importance_sampling.R index 46e61662..03d35a2e 100644 --- a/R/importance_sampling.R +++ b/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") { diff --git a/R/print.R b/R/print.R index 1f4030b8..6921d42a 100644 --- a/R/print.R +++ b/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") ) } diff --git a/R/psis.R b/R/psis.R index 303da755..156dc5c0 100644 --- a/R/psis.R +++ b/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 #' diff --git a/R/sis.R b/R/sis.R index fa64e358..84374650 100644 --- a/R/sis.R +++ b/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) + ) } diff --git a/R/tis.R b/R/tis.R index 39bf166d..8b226c1f 100644 --- a/R/tis.R +++ b/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)) } diff --git a/man/sis.Rd b/man/sis.Rd index d39af007..bb4f6202 100644 --- a/man/sis.Rd +++ b/man/sis.Rd @@ -61,7 +61,9 @@ class \code{sis}. \item{\code{diagnostics}}{ A named list containing one vector: \itemize{ -\item \code{pareto_k}: Not used in \code{sis}, all set to 0. +\item \code{pareto_k}: Estimates of the shape parameter \eqn{k} of the +generalized Pareto distribution, computed from the log ratios as in +\code{\link[=psis]{psis()}}. \code{sis()} uses them only as a diagnostic. \item \code{n_eff}: effective sample size estimates. } } @@ -77,7 +79,7 @@ used internally by the \code{weights} method to normalize the log weights. If specified, the user's \code{r_eff} argument. } \item{\code{tail_len}}{ -Not used for \code{sis}. +Vector of tail lengths used for the Pareto k diagnostic. } \item{\code{dims}}{ Integer vector of length 2 containing \code{S} (posterior sample size) diff --git a/man/tis.Rd b/man/tis.Rd index 1747a64f..c4a68fa2 100644 --- a/man/tis.Rd +++ b/man/tis.Rd @@ -63,7 +63,9 @@ class \code{tis}. \item{\code{diagnostics}}{ A named list containing one vector: \itemize{ -\item \code{pareto_k}: Not used in \code{tis}, all set to 0. +\item \code{pareto_k}: Estimates of the shape parameter \eqn{k} of the +generalized Pareto distribution, computed from the log ratios before +truncation as in \code{\link[=psis]{psis()}}. \code{tis()} uses them only as a diagnostic. \item \code{n_eff}: Effective sample size estimates. } } @@ -79,7 +81,7 @@ used internally by the \code{\link[=weights]{weights()}}method to normalize the If specified, the user's \code{r_eff} argument. } \item{\code{tail_len}}{ -Not used for \code{tis}. +Vector of tail lengths used for the Pareto k diagnostic. } \item{\code{dims}}{ Integer vector of length 2 containing \code{S} (posterior sample size) diff --git a/tests/testthat/_snaps/tisis.md b/tests/testthat/_snaps/tisis.md index c39c5e81..4b0112bf 100644 --- a/tests/testthat/_snaps/tisis.md +++ b/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 diff --git a/tests/testthat/test_tisis.R b/tests/testthat/test_tisis.R index 581d4898..580dd9b1 100644 --- a/tests/testthat/test_tisis.R +++ b/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) + } +})