From 894267b480ce2e0edc33d92c8d8b4720e9c73afb Mon Sep 17 00:00:00 2001 From: n-kall Date: Tue, 14 Nov 2023 12:17:08 +0200 Subject: [PATCH 01/14] Cleanup pareto documentation --- R/convergence.R | 2 ++ R/pareto_smooth.R | 2 ++ man/diagnostics.Rd | 2 ++ man/pareto_khat.Rd | 3 +++ 4 files changed, 9 insertions(+) diff --git a/R/convergence.R b/R/convergence.R index 01c7cd16..b3c83b8b 100644 --- a/R/convergence.R +++ b/R/convergence.R @@ -21,6 +21,8 @@ #' | [mcse_mean()] | Monte Carlo standard error for the mean | #' | [mcse_quantile()] | Monte Carlo standard error for quantiles | #' | [mcse_sd()] | Monte Carlo standard error for standard deviations | +#' | [pareto_khat()] | Pareto khat diagnostic for tail(s) | +#' | [pareto_diags()] | Additional diagnostics related to Pareto khat | #' | [rhat_basic()] | Basic version of Rhat | #' | [rhat()] | Improved, rank-based version of Rhat | #' | [rstar()] | R* diagnostic | diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 228aede9..295fe567 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -9,6 +9,8 @@ #' @template args-methods-dots #' @template ref-vehtari-paretosmooth-2022 #' @return `khat` estimated Generalized Pareto Distribution shape parameter k +#' +#' @seealso [`pareto_diags`] for additional related diagnostics. #' @examples #' mu <- extract_variable_matrix(example_draws(), "mu") #' pareto_khat(mu) diff --git a/man/diagnostics.Rd b/man/diagnostics.Rd index 2d5d7e31..b4ef67c5 100644 --- a/man/diagnostics.Rd +++ b/man/diagnostics.Rd @@ -21,6 +21,8 @@ A list of available diagnostics and links to their individual help pages. \code{\link[=mcse_mean]{mcse_mean()}} \tab Monte Carlo standard error for the mean \cr \code{\link[=mcse_quantile]{mcse_quantile()}} \tab Monte Carlo standard error for quantiles \cr \code{\link[=mcse_sd]{mcse_sd()}} \tab Monte Carlo standard error for standard deviations \cr + \code{\link[=pareto_khat]{pareto_khat()}} \tab Pareto khat diagnostic for tail(s) \cr + \code{\link[=pareto_diags]{pareto_diags()}} \tab Additional diagnostics related to Pareto khat \cr \code{\link[=rhat_basic]{rhat_basic()}} \tab Basic version of Rhat \cr \code{\link[=rhat]{rhat()}} \tab Improved, rank-based version of Rhat \cr \code{\link[=rstar]{rstar()}} \tab R* diagnostic \cr diff --git a/man/pareto_khat.Rd b/man/pareto_khat.Rd index c3383eb4..1a594085 100644 --- a/man/pareto_khat.Rd +++ b/man/pareto_khat.Rd @@ -72,3 +72,6 @@ Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao and Jonah Gabry (2022). Pareto Smoothed Importance Sampling. arxiv:arXiv:1507.02646 } +\seealso{ +\code{\link{pareto_diags}} for calculating additional related diagnostics. +} From 68617a52a917fad44acaab1da79542363679ff86 Mon Sep 17 00:00:00 2001 From: n-kall Date: Tue, 14 Nov 2023 17:25:02 +0200 Subject: [PATCH 02/14] add log option for pareto-smooth and smooth option to weight_draws --- R/pareto_smooth.R | 30 +++++++++++++++++++++++++----- R/weight_draws.R | 27 ++++++++++++++++++++++----- man/ess_basic.Rd | 2 ++ man/ess_bulk.Rd | 2 ++ man/ess_quantile.Rd | 2 ++ man/ess_sd.Rd | 2 ++ man/ess_tail.Rd | 2 ++ man/mcse_mean.Rd | 2 ++ man/mcse_quantile.Rd | 2 ++ man/mcse_sd.Rd | 2 ++ man/pareto_diags.Rd | 20 ++++++++++++++++++++ man/pareto_khat.Rd | 19 ++++++++++++++++++- man/pareto_smooth.Rd | 4 ++++ man/rhat.Rd | 2 ++ man/rhat_basic.Rd | 2 ++ man/rhat_nested.Rd | 2 ++ man/rstar.Rd | 2 ++ man/weight_draws.Rd | 15 +++++++++------ 18 files changed, 122 insertions(+), 17 deletions(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 295fe567..baaa5833 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -5,12 +5,14 @@ #' the number of fractional moments that is useful for convergence #' diagnostics. For further details see Vehtari et al. (2022). #' +#' @family diagnostics #' @template args-pareto #' @template args-methods-dots #' @template ref-vehtari-paretosmooth-2022 #' @return `khat` estimated Generalized Pareto Distribution shape parameter k #' -#' @seealso [`pareto_diags`] for additional related diagnostics. +#' @seealso [`pareto_diags`] for additional related diagnostics, and +#' [`pareto_smooth`] for Pareto smoothed draws. #' @examples #' mu <- extract_variable_matrix(example_draws(), "mu") #' pareto_khat(mu) @@ -67,6 +69,7 @@ pareto_khat.rvar <- function(x, ...) { #' replacing tail draws by order statistics of a generalized Pareto #' distribution fit to the tail(s). #' +#' @family diagnostics #' @template args-pareto #' @template args-methods-dots #' @template ref-vehtari-paretosmooth-2022 @@ -102,6 +105,8 @@ pareto_khat.rvar <- function(x, ...) { #' when the sample size is increased, compared to the central limit #' theorem convergence rate. See Appendix B in Vehtari et al. (2022). #' +#' @seealso [`pareto_khat`] for only calculating khat, and +#' [`pareto_smooth`] for Pareto smoothed draws. #' @examples #' mu <- extract_variable_matrix(example_draws(), "mu") #' pareto_diags(mu) @@ -191,6 +196,8 @@ pareto_diags.rvar <- function(x, ...) { #' Pareto smoothed estimates #' * `convergence_rate`: Relative convergence rate for Pareto smoothed estimates #' +#' @seealso [`pareto_khat`] for only calculating khat, and +#' [`pareto_diags`] for additional diagnostics. #' @examples #' mu <- extract_variable_matrix(example_draws(), "mu") #' pareto_smooth(mu) @@ -351,26 +358,32 @@ pareto_smooth.default <- function(x, ndraws_tail, smooth_draws = TRUE, tail = c("right", "left"), + log = FALSE, ... ) { + if (log) { + # shift log values for safe exponentiation + x <- x - max(x) + } + tail <- match.arg(tail) S <- length(x) tail_ids <- seq(S - ndraws_tail + 1, S) - if (tail == "left") { x <- -x } ord <- sort.int(x, index.return = TRUE) draws_tail <- ord$x[tail_ids] - cutoff <- ord$x[min(tail_ids) - 1] # largest value smaller than tail values + cutoff <- ord$x[min(tail_ids) - 1] # largest value smaller than tail values + max_tail <- max(draws_tail) min_tail <- min(draws_tail) - + if (ndraws_tail >= 5) { ord <- sort.int(x, index.return = TRUE) if (abs(max_tail - min_tail) < .Machine$double.eps / 100) { @@ -382,12 +395,19 @@ pareto_smooth.default <- function(x, k <- NA } else { # save time not sorting since x already sorted - fit <- gpdfit(draws_tail - cutoff, sort_x = FALSE) + if (log) { + draws_tail <- exp(draws_tail) + cutoff <- exp(cutoff) + } + fit <- gpdfit(draws_tail - cutoff, sort_x = FALSE, ...) k <- fit$k sigma <- fit$sigma if (is.finite(k) && smooth_draws) { p <- (seq_len(ndraws_tail) - 0.5) / ndraws_tail smoothed <- qgeneralized_pareto(p = p, mu = cutoff, k = k, sigma = sigma) + if (log) { + smoothed <- log(smoothed) + } } else { smoothed <- NULL } diff --git a/R/weight_draws.R b/R/weight_draws.R index 52326392..69debc1c 100644 --- a/R/weight_draws.R +++ b/R/weight_draws.R @@ -14,6 +14,8 @@ #' @param log (logical) Are the weights passed already on the log scale? The #' default is `FALSE`, that is, expecting `weights` to be on the standard #' (non-log) scale. +#' @param pareto_smooth (logical) Should the weights be Pareto-smoothed? +#' The default is `FALSE`. #' @template args-methods-dots #' @template return-draws #' @@ -50,9 +52,12 @@ weight_draws <- function(x, weights, ...) { #' @rdname weight_draws #' @export -weight_draws.draws_matrix <- function(x, weights, log = FALSE, ...) { +weight_draws.draws_matrix <- function(x, weights, log = FALSE, pareto_smooth = FALSE, ...) { log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) + if (pareto_smooth) { + log_weights <- pareto_smooth(log_weights, tail = "right", return_k = FALSE, log = TRUE) + } if (".log_weight" %in% variables(x, reserved = TRUE)) { # overwrite existing weights x[, ".log_weight"] <- log_weights @@ -66,9 +71,12 @@ weight_draws.draws_matrix <- function(x, weights, log = FALSE, ...) { #' @rdname weight_draws #' @export -weight_draws.draws_array <- function(x, weights, log = FALSE, ...) { +weight_draws.draws_array <- function(x, weights, log = FALSE, pareto_smooth = FALSE,...) { log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) + if (pareto_smooth) { + log_weights <- pareto_smooth(log_weights, tail = "right", return_k = FALSE, log = TRUE) + } if (".log_weight" %in% variables(x, reserved = TRUE)) { # overwrite existing weights x[, , ".log_weight"] <- log_weights @@ -82,18 +90,24 @@ weight_draws.draws_array <- function(x, weights, log = FALSE, ...) { #' @rdname weight_draws #' @export -weight_draws.draws_df <- function(x, weights, log = FALSE, ...) { +weight_draws.draws_df <- function(x, weights, log = FALSE, pareto_smooth = FALSE, ...) { log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) + if (pareto_smooth) { + log_weights <- pareto_smooth(log_weights, tail = "right", return_k = FALSE, log = TRUE) + } x$.log_weight <- log_weights x } #' @rdname weight_draws #' @export -weight_draws.draws_list <- function(x, weights, log = FALSE, ...) { +weight_draws.draws_list <- function(x, weights, log = FALSE, pareto_smooth = FALSE, ...) { log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) + if (pareto_smooth) { + log_weights <- pareto_smooth(log_weights, tail = "right", return_k = FALSE, log = TRUE) + } niterations <- niterations(x) for (i in seq_len(nchains(x))) { sel <- (1 + (i - 1) * niterations):(i * niterations) @@ -104,9 +118,12 @@ weight_draws.draws_list <- function(x, weights, log = FALSE, ...) { #' @rdname weight_draws #' @export -weight_draws.draws_rvars <- function(x, weights, log = FALSE, ...) { +weight_draws.draws_rvars <- function(x, weights, log = FALSE, pareto_smooth = FALSE, ...) { log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) + if (pareto_smooth) { + log_weights <- pareto_smooth(log_weights, tail = "right", return_k = FALSE, log = TRUE) + } x$.log_weight <- rvar(log_weights) x } diff --git a/man/ess_basic.Rd b/man/ess_basic.Rd index e300ad5e..867076ca 100755 --- a/man/ess_basic.Rd +++ b/man/ess_basic.Rd @@ -79,6 +79,8 @@ Other diagnostics: \code{\link{mcse_mean}()}, \code{\link{mcse_quantile}()}, \code{\link{mcse_sd}()}, +\code{\link{pareto_diags}()}, +\code{\link{pareto_khat}()}, \code{\link{rhat_basic}()}, \code{\link{rhat_nested}()}, \code{\link{rhat}()}, diff --git a/man/ess_bulk.Rd b/man/ess_bulk.Rd index adf3faf8..c1456be3 100755 --- a/man/ess_bulk.Rd +++ b/man/ess_bulk.Rd @@ -72,6 +72,8 @@ Other diagnostics: \code{\link{mcse_mean}()}, \code{\link{mcse_quantile}()}, \code{\link{mcse_sd}()}, +\code{\link{pareto_diags}()}, +\code{\link{pareto_khat}()}, \code{\link{rhat_basic}()}, \code{\link{rhat_nested}()}, \code{\link{rhat}()}, diff --git a/man/ess_quantile.Rd b/man/ess_quantile.Rd index 6bfc3cdf..aa85c909 100755 --- a/man/ess_quantile.Rd +++ b/man/ess_quantile.Rd @@ -81,6 +81,8 @@ Other diagnostics: \code{\link{mcse_mean}()}, \code{\link{mcse_quantile}()}, \code{\link{mcse_sd}()}, +\code{\link{pareto_diags}()}, +\code{\link{pareto_khat}()}, \code{\link{rhat_basic}()}, \code{\link{rhat_nested}()}, \code{\link{rhat}()}, diff --git a/man/ess_sd.Rd b/man/ess_sd.Rd index 2344211a..38475d2a 100755 --- a/man/ess_sd.Rd +++ b/man/ess_sd.Rd @@ -66,6 +66,8 @@ Other diagnostics: \code{\link{mcse_mean}()}, \code{\link{mcse_quantile}()}, \code{\link{mcse_sd}()}, +\code{\link{pareto_diags}()}, +\code{\link{pareto_khat}()}, \code{\link{rhat_basic}()}, \code{\link{rhat_nested}()}, \code{\link{rhat}()}, diff --git a/man/ess_tail.Rd b/man/ess_tail.Rd index f211f7aa..8f959718 100755 --- a/man/ess_tail.Rd +++ b/man/ess_tail.Rd @@ -72,6 +72,8 @@ Other diagnostics: \code{\link{mcse_mean}()}, \code{\link{mcse_quantile}()}, \code{\link{mcse_sd}()}, +\code{\link{pareto_diags}()}, +\code{\link{pareto_khat}()}, \code{\link{rhat_basic}()}, \code{\link{rhat_nested}()}, \code{\link{rhat}()}, diff --git a/man/mcse_mean.Rd b/man/mcse_mean.Rd index 9afaa7b7..c75935b1 100755 --- a/man/mcse_mean.Rd +++ b/man/mcse_mean.Rd @@ -63,6 +63,8 @@ Other diagnostics: \code{\link{ess_tail}()}, \code{\link{mcse_quantile}()}, \code{\link{mcse_sd}()}, +\code{\link{pareto_diags}()}, +\code{\link{pareto_khat}()}, \code{\link{rhat_basic}()}, \code{\link{rhat_nested}()}, \code{\link{rhat}()}, diff --git a/man/mcse_quantile.Rd b/man/mcse_quantile.Rd index cc4f9685..2d05f626 100755 --- a/man/mcse_quantile.Rd +++ b/man/mcse_quantile.Rd @@ -78,6 +78,8 @@ Other diagnostics: \code{\link{ess_tail}()}, \code{\link{mcse_mean}()}, \code{\link{mcse_sd}()}, +\code{\link{pareto_diags}()}, +\code{\link{pareto_khat}()}, \code{\link{rhat_basic}()}, \code{\link{rhat_nested}()}, \code{\link{rhat}()}, diff --git a/man/mcse_sd.Rd b/man/mcse_sd.Rd index 7e322864..671ef249 100755 --- a/man/mcse_sd.Rd +++ b/man/mcse_sd.Rd @@ -68,6 +68,8 @@ Other diagnostics: \code{\link{ess_tail}()}, \code{\link{mcse_mean}()}, \code{\link{mcse_quantile}()}, +\code{\link{pareto_diags}()}, +\code{\link{pareto_khat}()}, \code{\link{rhat_basic}()}, \code{\link{rhat_nested}()}, \code{\link{rhat}()}, diff --git a/man/pareto_diags.Rd b/man/pareto_diags.Rd index c2558a9a..38e55b5a 100644 --- a/man/pareto_diags.Rd +++ b/man/pareto_diags.Rd @@ -103,3 +103,23 @@ Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao and Jonah Gabry (2022). Pareto Smoothed Importance Sampling. arxiv:arXiv:1507.02646 } +\seealso{ +\code{\link{pareto_khat}} for only calculating khat, and +\code{\link{pareto_smooth}} for Pareto smoothed draws. + +Other diagnostics: +\code{\link{ess_basic}()}, +\code{\link{ess_bulk}()}, +\code{\link{ess_quantile}()}, +\code{\link{ess_sd}()}, +\code{\link{ess_tail}()}, +\code{\link{mcse_mean}()}, +\code{\link{mcse_quantile}()}, +\code{\link{mcse_sd}()}, +\code{\link{pareto_khat}()}, +\code{\link{rhat_basic}()}, +\code{\link{rhat_nested}()}, +\code{\link{rhat}()}, +\code{\link{rstar}()} +} +\concept{diagnostics} diff --git a/man/pareto_khat.Rd b/man/pareto_khat.Rd index 1a594085..1c0d0e92 100644 --- a/man/pareto_khat.Rd +++ b/man/pareto_khat.Rd @@ -73,5 +73,22 @@ Jonah Gabry (2022). Pareto Smoothed Importance Sampling. arxiv:arXiv:1507.02646 } \seealso{ -\code{\link{pareto_diags}} for calculating additional related diagnostics. +\code{\link{pareto_diags}} for additional related diagnostics, and +\code{\link{pareto_smooth}} for Pareto smoothed draws. + +Other diagnostics: +\code{\link{ess_basic}()}, +\code{\link{ess_bulk}()}, +\code{\link{ess_quantile}()}, +\code{\link{ess_sd}()}, +\code{\link{ess_tail}()}, +\code{\link{mcse_mean}()}, +\code{\link{mcse_quantile}()}, +\code{\link{mcse_sd}()}, +\code{\link{pareto_diags}()}, +\code{\link{rhat_basic}()}, +\code{\link{rhat_nested}()}, +\code{\link{rhat}()}, +\code{\link{rstar}()} } +\concept{diagnostics} diff --git a/man/pareto_smooth.Rd b/man/pareto_smooth.Rd index 259421de..4a62a547 100644 --- a/man/pareto_smooth.Rd +++ b/man/pareto_smooth.Rd @@ -93,3 +93,7 @@ Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao and Jonah Gabry (2022). Pareto Smoothed Importance Sampling. arxiv:arXiv:1507.02646 } +\seealso{ +\code{\link{pareto_khat}} for only calculating khat, and +\code{\link{pareto_diags}} for additional diagnostics. +} diff --git a/man/rhat.Rd b/man/rhat.Rd index fed2c14a..263561cd 100755 --- a/man/rhat.Rd +++ b/man/rhat.Rd @@ -67,6 +67,8 @@ Other diagnostics: \code{\link{mcse_mean}()}, \code{\link{mcse_quantile}()}, \code{\link{mcse_sd}()}, +\code{\link{pareto_diags}()}, +\code{\link{pareto_khat}()}, \code{\link{rhat_basic}()}, \code{\link{rhat_nested}()}, \code{\link{rstar}()} diff --git a/man/rhat_basic.Rd b/man/rhat_basic.Rd index 8a94efb3..16ffd332 100755 --- a/man/rhat_basic.Rd +++ b/man/rhat_basic.Rd @@ -75,6 +75,8 @@ Other diagnostics: \code{\link{mcse_mean}()}, \code{\link{mcse_quantile}()}, \code{\link{mcse_sd}()}, +\code{\link{pareto_diags}()}, +\code{\link{pareto_khat}()}, \code{\link{rhat_nested}()}, \code{\link{rhat}()}, \code{\link{rstar}()} diff --git a/man/rhat_nested.Rd b/man/rhat_nested.Rd index 2e23242d..f2536efd 100644 --- a/man/rhat_nested.Rd +++ b/man/rhat_nested.Rd @@ -83,6 +83,8 @@ Other diagnostics: \code{\link{mcse_mean}()}, \code{\link{mcse_quantile}()}, \code{\link{mcse_sd}()}, +\code{\link{pareto_diags}()}, +\code{\link{pareto_khat}()}, \code{\link{rhat_basic}()}, \code{\link{rhat}()}, \code{\link{rstar}()} diff --git a/man/rstar.Rd b/man/rstar.Rd index 87e8e372..c9479902 100644 --- a/man/rstar.Rd +++ b/man/rstar.Rd @@ -115,6 +115,8 @@ Other diagnostics: \code{\link{mcse_mean}()}, \code{\link{mcse_quantile}()}, \code{\link{mcse_sd}()}, +\code{\link{pareto_diags}()}, +\code{\link{pareto_khat}()}, \code{\link{rhat_basic}()}, \code{\link{rhat_nested}()}, \code{\link{rhat}()} diff --git a/man/weight_draws.Rd b/man/weight_draws.Rd index ecd611c9..dfc5d2ba 100644 --- a/man/weight_draws.Rd +++ b/man/weight_draws.Rd @@ -11,15 +11,15 @@ \usage{ weight_draws(x, weights, ...) -\method{weight_draws}{draws_matrix}(x, weights, log = FALSE, ...) +\method{weight_draws}{draws_matrix}(x, weights, log = FALSE, pareto_smooth = FALSE, ...) -\method{weight_draws}{draws_array}(x, weights, log = FALSE, ...) +\method{weight_draws}{draws_array}(x, weights, log = FALSE, pareto_smooth = FALSE, ...) -\method{weight_draws}{draws_df}(x, weights, log = FALSE, ...) +\method{weight_draws}{draws_df}(x, weights, log = FALSE, pareto_smooth = FALSE, ...) -\method{weight_draws}{draws_list}(x, weights, log = FALSE, ...) +\method{weight_draws}{draws_list}(x, weights, log = FALSE, pareto_smooth = FALSE, ...) -\method{weight_draws}{draws_rvars}(x, weights, log = FALSE, ...) +\method{weight_draws}{draws_rvars}(x, weights, log = FALSE, pareto_smooth = FALSE, ...) } \arguments{ \item{x}{(draws) A \code{draws} object or another \R object for which the method @@ -32,9 +32,12 @@ can be returned via the \code{\link[=weights.draws]{weights.draws()}} method lat \item{...}{Arguments passed to individual methods (if applicable).} -\item{log}{(logicla) Are the weights passed already on the log scale? The +\item{log}{(logical) Are the weights passed already on the log scale? The default is \code{FALSE}, that is, expecting \code{weights} to be on the standard (non-log) scale.} + +\item{pareto_smooth}{(logical) Should the weights be Pareto-smoothed? +The default is \code{FALSE}.} } \value{ A \code{draws} object of the same class as \code{x}. From 16843c1f0f2d641297897cfa81aaced878bfb37e Mon Sep 17 00:00:00 2001 From: n-kall Date: Wed, 15 Nov 2023 13:27:08 +0200 Subject: [PATCH 03/14] add log_weights option to pareto functions --- R/pareto_smooth.R | 40 +++++++++++++++++++++++++-------------- man-roxygen/args-pareto.R | 6 +++++- man/pareto_diags.Rd | 8 +++++++- man/pareto_khat.Rd | 8 +++++++- man/pareto_smooth.Rd | 8 +++++++- 5 files changed, 52 insertions(+), 18 deletions(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index baaa5833..27b6895a 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -29,6 +29,7 @@ pareto_khat.default <- function(x, r_eff = NULL, ndraws_tail = NULL, verbose = FALSE, + log_weights = FALSE, ...) { smoothed <- pareto_smooth.default( x, @@ -38,6 +39,7 @@ pareto_khat.default <- function(x, verbose = verbose, return_k = TRUE, smooth_draws = FALSE, + log_weights = log_weights, ...) return(smoothed$diagnostics) } @@ -120,11 +122,12 @@ pareto_diags <- function(x, ...) UseMethod("pareto_diags") #' @rdname pareto_diags #' @export pareto_diags.default <- function(x, - tail = c("both", "right", "left"), - r_eff = NULL, - ndraws_tail = NULL, - verbose = FALSE, - ...) { + tail = c("both", "right", "left"), + r_eff = NULL, + ndraws_tail = NULL, + verbose = FALSE, + log_weights = FALSE, + ...) { smoothed <- pareto_smooth.default( x, @@ -135,6 +138,7 @@ pareto_diags.default <- function(x, extra_diags = TRUE, verbose = verbose, smooth_draws = FALSE, + log_weights = FALSE, ...) return(smoothed$diagnostics) @@ -234,8 +238,8 @@ pareto_smooth.rvar <- function(x, return_k = TRUE, extra_diags = FALSE, ...) { ) } out <- list( - x = rvar(apply(draws_diags, margins, function(x) x[[1]]$x), nchains = nchains(x)), - diagnostics = diags + x = rvar(apply(draws_diags, margins, function(x) x[[1]]$x), nchains = nchains(x)), + diagnostics = diags ) } else { out <- rvar(apply(draws_diags, margins, function(x) x[[1]]), nchains = nchains(x)) @@ -252,6 +256,7 @@ pareto_smooth.default <- function(x, return_k = TRUE, extra_diags = FALSE, verbose = FALSE, + log_weights = FALSE, ...) { checkmate::assert_number(ndraws_tail, null.ok = TRUE) @@ -259,6 +264,7 @@ pareto_smooth.default <- function(x, checkmate::assert_logical(extra_diags) checkmate::assert_logical(return_k) checkmate::assert_logical(verbose) + checkmate::assert_logical(log_weights) # check for infinite or na values if (should_return_NA(x)) { @@ -266,6 +272,10 @@ pareto_smooth.default <- function(x, return(list(x = x, diagnostics = NA_real_)) } + if (log_weights) { + tail = "right" + } + tail <- match.arg(tail) S <- length(x) @@ -299,6 +309,7 @@ pareto_smooth.default <- function(x, x, ndraws_tail = ndraws_tail, tail = "left", + log_weights = log_weights, ... ) left_k <- smoothed$k @@ -308,6 +319,7 @@ pareto_smooth.default <- function(x, x = smoothed$x, ndraws_tail = ndraws_tail, tail = "right", + log_weights = log_weights, ... ) right_k <- smoothed$k @@ -358,11 +370,11 @@ pareto_smooth.default <- function(x, ndraws_tail, smooth_draws = TRUE, tail = c("right", "left"), - log = FALSE, + log_weights = FALSE, ... ) { - if (log) { + if (log_weights) { # shift log values for safe exponentiation x <- x - max(x) } @@ -395,7 +407,7 @@ pareto_smooth.default <- function(x, k <- NA } else { # save time not sorting since x already sorted - if (log) { + if (log_weights) { draws_tail <- exp(draws_tail) cutoff <- exp(cutoff) } @@ -405,7 +417,7 @@ pareto_smooth.default <- function(x, if (is.finite(k) && smooth_draws) { p <- (seq_len(ndraws_tail) - 0.5) / ndraws_tail smoothed <- qgeneralized_pareto(p = p, mu = cutoff, k = k, sigma = sigma) - if (log) { + if (log_weights) { smoothed <- log(smoothed) } } else { @@ -467,11 +479,11 @@ pareto_smooth.default <- function(x, #' @noRd ps_min_ss <- function(k, ...) { if (k < 1) { - out <- 10^(1 / (1 - max(0, k))) + out <- 10^(1 / (1 - max(0, k))) } else { - out <- Inf + out <- Inf } - out + out } diff --git a/man-roxygen/args-pareto.R b/man-roxygen/args-pareto.R index a406dbde..070ded54 100644 --- a/man-roxygen/args-pareto.R +++ b/man-roxygen/args-pareto.R @@ -11,10 +11,14 @@ #' @param ndraws_tail (numeric) number of draws for the tail. If #' `ndraws_tail` is not specified, it will be calculated as #' ceiling(3 * sqrt(length(x) / r_eff)) if length(x) > 225 and -#' length(x) / 5 otherwise (see Appendix H in Vehtari et al. (2022)). +#' length(x) / 5 otherwise (see Appendix H in Vehtari et +#' al. (2022)). #' @param r_eff (numeric) relative effective sample size estimate. If #' `r_eff` is omitted, it will be calculated assuming the draws are #' from MCMC. #' @param verbose (logical) Should diagnostic messages be printed? If #' `TRUE`, messages related to Pareto diagnostics will be #' printed. Default is `FALSE`. +#' @param log_weights (logical) Are the draws log weights? Default is +#' `FALSE`. If `TRUE` computation will take into account that the +#' draws are log weights, and only right tail will be smoothed. diff --git a/man/pareto_diags.Rd b/man/pareto_diags.Rd index 38e55b5a..b47cab5d 100644 --- a/man/pareto_diags.Rd +++ b/man/pareto_diags.Rd @@ -14,6 +14,7 @@ pareto_diags(x, ...) r_eff = NULL, ndraws_tail = NULL, verbose = FALSE, + log_weights = FALSE, ... ) @@ -45,11 +46,16 @@ from MCMC.} \item{ndraws_tail}{(numeric) number of draws for the tail. If \code{ndraws_tail} is not specified, it will be calculated as ceiling(3 * sqrt(length(x) / r_eff)) if length(x) > 225 and -length(x) / 5 otherwise (see Appendix H in Vehtari et al. (2022)).} +length(x) / 5 otherwise (see Appendix H in Vehtari et +al. (2022)).} \item{verbose}{(logical) Should diagnostic messages be printed? If \code{TRUE}, messages related to Pareto diagnostics will be printed. Default is \code{FALSE}.} + +\item{log_weights}{(logical) Are the draws log weights? Default is +\code{FALSE}. If \code{TRUE} computation will take into account that the +draws are log weights, and only right tail will be smoothed.} } \value{ List of Pareto smoothing diagnostics: diff --git a/man/pareto_khat.Rd b/man/pareto_khat.Rd index 1c0d0e92..9094fd14 100644 --- a/man/pareto_khat.Rd +++ b/man/pareto_khat.Rd @@ -14,6 +14,7 @@ pareto_khat(x, ...) r_eff = NULL, ndraws_tail = NULL, verbose = FALSE, + log_weights = FALSE, ... ) @@ -45,11 +46,16 @@ from MCMC.} \item{ndraws_tail}{(numeric) number of draws for the tail. If \code{ndraws_tail} is not specified, it will be calculated as ceiling(3 * sqrt(length(x) / r_eff)) if length(x) > 225 and -length(x) / 5 otherwise (see Appendix H in Vehtari et al. (2022)).} +length(x) / 5 otherwise (see Appendix H in Vehtari et +al. (2022)).} \item{verbose}{(logical) Should diagnostic messages be printed? If \code{TRUE}, messages related to Pareto diagnostics will be printed. Default is \code{FALSE}.} + +\item{log_weights}{(logical) Are the draws log weights? Default is +\code{FALSE}. If \code{TRUE} computation will take into account that the +draws are log weights, and only right tail will be smoothed.} } \value{ \code{khat} estimated Generalized Pareto Distribution shape parameter k diff --git a/man/pareto_smooth.Rd b/man/pareto_smooth.Rd index 4a62a547..cd335da0 100644 --- a/man/pareto_smooth.Rd +++ b/man/pareto_smooth.Rd @@ -18,6 +18,7 @@ pareto_smooth(x, ...) return_k = TRUE, extra_diags = FALSE, verbose = FALSE, + log_weights = FALSE, ... ) } @@ -56,11 +57,16 @@ from MCMC.} \item{ndraws_tail}{(numeric) number of draws for the tail. If \code{ndraws_tail} is not specified, it will be calculated as ceiling(3 * sqrt(length(x) / r_eff)) if length(x) > 225 and -length(x) / 5 otherwise (see Appendix H in Vehtari et al. (2022)).} +length(x) / 5 otherwise (see Appendix H in Vehtari et +al. (2022)).} \item{verbose}{(logical) Should diagnostic messages be printed? If \code{TRUE}, messages related to Pareto diagnostics will be printed. Default is \code{FALSE}.} + +\item{log_weights}{(logical) Are the draws log weights? Default is +\code{FALSE}. If \code{TRUE} computation will take into account that the +draws are log weights, and only right tail will be smoothed.} } \value{ Either a vector \code{x} of smoothed values or a named list From 1575011c99787f4d8ef7d07a11a5ba23b344cf27 Mon Sep 17 00:00:00 2001 From: n-kall Date: Wed, 15 Nov 2023 13:38:01 +0200 Subject: [PATCH 04/14] specify version number of psis preprint --- man-roxygen/ref-vehtari-paretosmooth-2022.R | 2 +- man/pareto_diags.Rd | 2 +- man/pareto_khat.Rd | 2 +- man/pareto_smooth.Rd | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/man-roxygen/ref-vehtari-paretosmooth-2022.R b/man-roxygen/ref-vehtari-paretosmooth-2022.R index 267ec29a..30f526fc 100644 --- a/man-roxygen/ref-vehtari-paretosmooth-2022.R +++ b/man-roxygen/ref-vehtari-paretosmooth-2022.R @@ -1,4 +1,4 @@ #' @references #' Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao and #' Jonah Gabry (2022). Pareto Smoothed Importance Sampling. -#' arxiv:arXiv:1507.02646 +#' arxiv:arXiv:1507.02646 (version 8) diff --git a/man/pareto_diags.Rd b/man/pareto_diags.Rd index b47cab5d..ad1f011a 100644 --- a/man/pareto_diags.Rd +++ b/man/pareto_diags.Rd @@ -107,7 +107,7 @@ pareto_diags(d$Sigma) \references{ Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao and Jonah Gabry (2022). Pareto Smoothed Importance Sampling. -arxiv:arXiv:1507.02646 +arxiv:arXiv:1507.02646 (version 8) } \seealso{ \code{\link{pareto_khat}} for only calculating khat, and diff --git a/man/pareto_khat.Rd b/man/pareto_khat.Rd index 9094fd14..3be8d501 100644 --- a/man/pareto_khat.Rd +++ b/man/pareto_khat.Rd @@ -76,7 +76,7 @@ pareto_khat(d$Sigma) \references{ Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao and Jonah Gabry (2022). Pareto Smoothed Importance Sampling. -arxiv:arXiv:1507.02646 +arxiv:arXiv:1507.02646 (version 8) } \seealso{ \code{\link{pareto_diags}} for additional related diagnostics, and diff --git a/man/pareto_smooth.Rd b/man/pareto_smooth.Rd index cd335da0..88f8dc63 100644 --- a/man/pareto_smooth.Rd +++ b/man/pareto_smooth.Rd @@ -97,7 +97,7 @@ pareto_smooth(d$Sigma) \references{ Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao and Jonah Gabry (2022). Pareto Smoothed Importance Sampling. -arxiv:arXiv:1507.02646 +arxiv:arXiv:1507.02646 (version 8) } \seealso{ \code{\link{pareto_khat}} for only calculating khat, and From d431e5afb58b0076920ea8aa48f0765187cc5108 Mon Sep 17 00:00:00 2001 From: n-kall Date: Wed, 15 Nov 2023 18:28:48 +0200 Subject: [PATCH 05/14] validate pareto_smooth flag in weight_draws --- R/weight_draws.R | 28 +++++++++++++++++++++++----- 1 file changed, 23 insertions(+), 5 deletions(-) diff --git a/R/weight_draws.R b/R/weight_draws.R index 69debc1c..64624d6d 100644 --- a/R/weight_draws.R +++ b/R/weight_draws.R @@ -45,6 +45,9 @@ #' head(weights(x)) #' head(weights(x, log=TRUE, normalize = FALSE)) # recover original log_wts #' +#' # add weights on log scale and Pareto smooth them +#' x <- weight_draws(x, weights = log_weights, log = TRUE, pareto_smooth = TRUE) +#' #' @export weight_draws <- function(x, weights, ...) { UseMethod("weight_draws") @@ -53,11 +56,14 @@ weight_draws <- function(x, weights, ...) { #' @rdname weight_draws #' @export weight_draws.draws_matrix <- function(x, weights, log = FALSE, pareto_smooth = FALSE, ...) { + + checkmate::assert_logical(pareto_smooth) + log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { log_weights <- pareto_smooth(log_weights, tail = "right", return_k = FALSE, log = TRUE) - } + } if (".log_weight" %in% variables(x, reserved = TRUE)) { # overwrite existing weights x[, ".log_weight"] <- log_weights @@ -72,11 +78,14 @@ weight_draws.draws_matrix <- function(x, weights, log = FALSE, pareto_smooth = F #' @rdname weight_draws #' @export weight_draws.draws_array <- function(x, weights, log = FALSE, pareto_smooth = FALSE,...) { + + checkmate::assert_logical(pareto_smooth) + log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { log_weights <- pareto_smooth(log_weights, tail = "right", return_k = FALSE, log = TRUE) - } + } if (".log_weight" %in% variables(x, reserved = TRUE)) { # overwrite existing weights x[, , ".log_weight"] <- log_weights @@ -91,11 +100,14 @@ weight_draws.draws_array <- function(x, weights, log = FALSE, pareto_smooth = FA #' @rdname weight_draws #' @export weight_draws.draws_df <- function(x, weights, log = FALSE, pareto_smooth = FALSE, ...) { + + checkmate::assert_logical(pareto_smooth) + log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { log_weights <- pareto_smooth(log_weights, tail = "right", return_k = FALSE, log = TRUE) - } + } x$.log_weight <- log_weights x } @@ -103,11 +115,14 @@ weight_draws.draws_df <- function(x, weights, log = FALSE, pareto_smooth = FALSE #' @rdname weight_draws #' @export weight_draws.draws_list <- function(x, weights, log = FALSE, pareto_smooth = FALSE, ...) { + + checkmate::assert_logical(pareto_smooth) + log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { log_weights <- pareto_smooth(log_weights, tail = "right", return_k = FALSE, log = TRUE) - } + } niterations <- niterations(x) for (i in seq_len(nchains(x))) { sel <- (1 + (i - 1) * niterations):(i * niterations) @@ -119,11 +134,14 @@ weight_draws.draws_list <- function(x, weights, log = FALSE, pareto_smooth = FAL #' @rdname weight_draws #' @export weight_draws.draws_rvars <- function(x, weights, log = FALSE, pareto_smooth = FALSE, ...) { + + checkmate::assert_logical(pareto_smooth) + log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { log_weights <- pareto_smooth(log_weights, tail = "right", return_k = FALSE, log = TRUE) - } + } x$.log_weight <- rvar(log_weights) x } From 3c67dd79e9a4aedea08998f27d75da43d573938b Mon Sep 17 00:00:00 2001 From: n-kall Date: Wed, 15 Nov 2023 18:29:31 +0200 Subject: [PATCH 06/14] ad pareto smoothing weights example --- R/weight_draws.R | 2 +- man/weight_draws.Rd | 3 +++ man/weights.draws.Rd | 3 +++ 3 files changed, 7 insertions(+), 1 deletion(-) diff --git a/R/weight_draws.R b/R/weight_draws.R index 64624d6d..3bd1068e 100644 --- a/R/weight_draws.R +++ b/R/weight_draws.R @@ -46,7 +46,7 @@ #' head(weights(x, log=TRUE, normalize = FALSE)) # recover original log_wts #' #' # add weights on log scale and Pareto smooth them -#' x <- weight_draws(x, weights = log_weights, log = TRUE, pareto_smooth = TRUE) +#' x <- weight_draws(x, weights = log_wts, log = TRUE, pareto_smooth = TRUE) #' #' @export weight_draws <- function(x, weights, ...) { diff --git a/man/weight_draws.Rd b/man/weight_draws.Rd index dfc5d2ba..d866d466 100644 --- a/man/weight_draws.Rd +++ b/man/weight_draws.Rd @@ -73,6 +73,9 @@ x <- weight_draws(x, weights = log_wts, log = TRUE) head(weights(x)) head(weights(x, log=TRUE, normalize = FALSE)) # recover original log_wts +# add weights on log scale and Pareto smooth them +x <- weight_draws(x, weights = log_wts, log = TRUE, pareto_smooth = TRUE) + } \seealso{ \code{\link[=weights.draws]{weights.draws()}}, \code{\link[=resample_draws]{resample_draws()}} diff --git a/man/weights.draws.Rd b/man/weights.draws.Rd index 6b2a46ff..1a47788e 100644 --- a/man/weights.draws.Rd +++ b/man/weights.draws.Rd @@ -48,6 +48,9 @@ x <- weight_draws(x, weights = log_wts, log = TRUE) head(weights(x)) head(weights(x, log=TRUE, normalize = FALSE)) # recover original log_wts +# add weights on log scale and Pareto smooth them +x <- weight_draws(x, weights = log_wts, log = TRUE, pareto_smooth = TRUE) + } \seealso{ \code{\link{weight_draws}}, \code{\link{resample_draws}} From 64e755894e3bd2aab1e40ba25ef5405a79451b4b Mon Sep 17 00:00:00 2001 From: n-kall Date: Wed, 15 Nov 2023 18:29:45 +0200 Subject: [PATCH 07/14] fix output structure for pareto_smooth when NA --- R/pareto_smooth.R | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 27b6895a..8281fa3d 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -268,8 +268,13 @@ pareto_smooth.default <- function(x, # check for infinite or na values if (should_return_NA(x)) { - warning_no_call("Input contains infinite or NA values, Pareto smoothing not performed.") - return(list(x = x, diagnostics = NA_real_)) + warning_no_call("Input contains infinite or NA values, or is constant. Fitting of generalized Pareto distribution not performed.") + if (!return_k) { + out <- x + } else { + out <- list(x = x, diagnostics = NA_real_) + } + return(out) } if (log_weights) { From da23af6b98a2bd0075b1bc72cd05c0fcac245960 Mon Sep 17 00:00:00 2001 From: n-kall Date: Wed, 15 Nov 2023 18:40:17 +0200 Subject: [PATCH 08/14] change default for r_eff to 1 in pareto_smooth (and others) --- R/pareto_smooth.R | 2 +- man-roxygen/args-pareto.R | 2 +- man/pareto_diags.Rd | 2 +- man/pareto_khat.Rd | 2 +- man/pareto_smooth.Rd | 4 ++-- 5 files changed, 6 insertions(+), 6 deletions(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 8281fa3d..7111821b 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -251,7 +251,7 @@ pareto_smooth.rvar <- function(x, return_k = TRUE, extra_diags = FALSE, ...) { #' @export pareto_smooth.default <- function(x, tail = c("both", "right", "left"), - r_eff = NULL, + r_eff = 1, ndraws_tail = NULL, return_k = TRUE, extra_diags = FALSE, diff --git a/man-roxygen/args-pareto.R b/man-roxygen/args-pareto.R index 070ded54..98940a7e 100644 --- a/man-roxygen/args-pareto.R +++ b/man-roxygen/args-pareto.R @@ -15,7 +15,7 @@ #' al. (2022)). #' @param r_eff (numeric) relative effective sample size estimate. If #' `r_eff` is omitted, it will be calculated assuming the draws are -#' from MCMC. +#' from MCMC. Default is 1. #' @param verbose (logical) Should diagnostic messages be printed? If #' `TRUE`, messages related to Pareto diagnostics will be #' printed. Default is `FALSE`. diff --git a/man/pareto_diags.Rd b/man/pareto_diags.Rd index ad1f011a..7f850ac2 100644 --- a/man/pareto_diags.Rd +++ b/man/pareto_diags.Rd @@ -41,7 +41,7 @@ The default is \code{"both"}.} \item{r_eff}{(numeric) relative effective sample size estimate. If \code{r_eff} is omitted, it will be calculated assuming the draws are -from MCMC.} +from MCMC. Default is 1.} \item{ndraws_tail}{(numeric) number of draws for the tail. If \code{ndraws_tail} is not specified, it will be calculated as diff --git a/man/pareto_khat.Rd b/man/pareto_khat.Rd index 3be8d501..cd5c11a9 100644 --- a/man/pareto_khat.Rd +++ b/man/pareto_khat.Rd @@ -41,7 +41,7 @@ The default is \code{"both"}.} \item{r_eff}{(numeric) relative effective sample size estimate. If \code{r_eff} is omitted, it will be calculated assuming the draws are -from MCMC.} +from MCMC. Default is 1.} \item{ndraws_tail}{(numeric) number of draws for the tail. If \code{ndraws_tail} is not specified, it will be calculated as diff --git a/man/pareto_smooth.Rd b/man/pareto_smooth.Rd index 88f8dc63..f2bd5e34 100644 --- a/man/pareto_smooth.Rd +++ b/man/pareto_smooth.Rd @@ -13,7 +13,7 @@ pareto_smooth(x, ...) \method{pareto_smooth}{default}( x, tail = c("both", "right", "left"), - r_eff = NULL, + r_eff = 1, ndraws_tail = NULL, return_k = TRUE, extra_diags = FALSE, @@ -52,7 +52,7 @@ The default is \code{"both"}.} \item{r_eff}{(numeric) relative effective sample size estimate. If \code{r_eff} is omitted, it will be calculated assuming the draws are -from MCMC.} +from MCMC. Default is 1.} \item{ndraws_tail}{(numeric) number of draws for the tail. If \code{ndraws_tail} is not specified, it will be calculated as From 8adfe47fba2c6b37a7f82b56c57016b46decdaeb Mon Sep 17 00:00:00 2001 From: n-kall Date: Wed, 15 Nov 2023 18:43:55 +0200 Subject: [PATCH 09/14] change default for r_eff in pareto_smooth --- man-roxygen/args-pareto.R | 2 +- man/pareto_diags.Rd | 2 +- man/pareto_khat.Rd | 2 +- man/pareto_smooth.Rd | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/man-roxygen/args-pareto.R b/man-roxygen/args-pareto.R index 98940a7e..db89f30c 100644 --- a/man-roxygen/args-pareto.R +++ b/man-roxygen/args-pareto.R @@ -14,7 +14,7 @@ #' length(x) / 5 otherwise (see Appendix H in Vehtari et #' al. (2022)). #' @param r_eff (numeric) relative effective sample size estimate. If -#' `r_eff` is omitted, it will be calculated assuming the draws are +#' `r_eff` is NULL, it will be calculated assuming the draws are #' from MCMC. Default is 1. #' @param verbose (logical) Should diagnostic messages be printed? If #' `TRUE`, messages related to Pareto diagnostics will be diff --git a/man/pareto_diags.Rd b/man/pareto_diags.Rd index 7f850ac2..187a1237 100644 --- a/man/pareto_diags.Rd +++ b/man/pareto_diags.Rd @@ -40,7 +40,7 @@ pareto_diags(x, ...) The default is \code{"both"}.} \item{r_eff}{(numeric) relative effective sample size estimate. If -\code{r_eff} is omitted, it will be calculated assuming the draws are +\code{r_eff} is NULL, it will be calculated assuming the draws are from MCMC. Default is 1.} \item{ndraws_tail}{(numeric) number of draws for the tail. If diff --git a/man/pareto_khat.Rd b/man/pareto_khat.Rd index cd5c11a9..0d844580 100644 --- a/man/pareto_khat.Rd +++ b/man/pareto_khat.Rd @@ -40,7 +40,7 @@ pareto_khat(x, ...) The default is \code{"both"}.} \item{r_eff}{(numeric) relative effective sample size estimate. If -\code{r_eff} is omitted, it will be calculated assuming the draws are +\code{r_eff} is NULL, it will be calculated assuming the draws are from MCMC. Default is 1.} \item{ndraws_tail}{(numeric) number of draws for the tail. If diff --git a/man/pareto_smooth.Rd b/man/pareto_smooth.Rd index f2bd5e34..3e5e444d 100644 --- a/man/pareto_smooth.Rd +++ b/man/pareto_smooth.Rd @@ -51,7 +51,7 @@ returned. Default is \code{FALSE}.} The default is \code{"both"}.} \item{r_eff}{(numeric) relative effective sample size estimate. If -\code{r_eff} is omitted, it will be calculated assuming the draws are +\code{r_eff} is NULL, it will be calculated assuming the draws are from MCMC. Default is 1.} \item{ndraws_tail}{(numeric) number of draws for the tail. If From b39fdc2b1297da395d54663e64b9503aea72cdf4 Mon Sep 17 00:00:00 2001 From: n-kall Date: Fri, 17 Nov 2023 14:28:18 +0200 Subject: [PATCH 10/14] change pareto_smooth flag check in weight_draws --- R/weight_draws.R | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/R/weight_draws.R b/R/weight_draws.R index 3bd1068e..9270aa14 100644 --- a/R/weight_draws.R +++ b/R/weight_draws.R @@ -57,8 +57,8 @@ weight_draws <- function(x, weights, ...) { #' @export weight_draws.draws_matrix <- function(x, weights, log = FALSE, pareto_smooth = FALSE, ...) { - checkmate::assert_logical(pareto_smooth) + pareto_smooth <- as_one_logical(pareto_smooth) log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { @@ -79,8 +79,7 @@ weight_draws.draws_matrix <- function(x, weights, log = FALSE, pareto_smooth = F #' @export weight_draws.draws_array <- function(x, weights, log = FALSE, pareto_smooth = FALSE,...) { - checkmate::assert_logical(pareto_smooth) - + pareto_smooth <- as_one_logical(pareto_smooth) log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { @@ -101,8 +100,7 @@ weight_draws.draws_array <- function(x, weights, log = FALSE, pareto_smooth = FA #' @export weight_draws.draws_df <- function(x, weights, log = FALSE, pareto_smooth = FALSE, ...) { - checkmate::assert_logical(pareto_smooth) - + pareto_smooth <- as_one_logical(pareto_smooth) log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { @@ -116,8 +114,7 @@ weight_draws.draws_df <- function(x, weights, log = FALSE, pareto_smooth = FALSE #' @export weight_draws.draws_list <- function(x, weights, log = FALSE, pareto_smooth = FALSE, ...) { - checkmate::assert_logical(pareto_smooth) - + pareto_smooth <- as_one_logical(pareto_smooth) log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { @@ -135,8 +132,7 @@ weight_draws.draws_list <- function(x, weights, log = FALSE, pareto_smooth = FAL #' @export weight_draws.draws_rvars <- function(x, weights, log = FALSE, pareto_smooth = FALSE, ...) { - checkmate::assert_logical(pareto_smooth) - + pareto_smooth <- as_one_logical(pareto_smooth) log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { From 2bfd8182fc8a4f5c70bfc5d6fe63d21c1a70a287 Mon Sep 17 00:00:00 2001 From: n-kall Date: Fri, 17 Nov 2023 16:20:04 +0200 Subject: [PATCH 11/14] improve pareto diagnostic warning messages, add log pareto test --- R/pareto_smooth.R | 77 +++++++++++++++++------------ R/weight_draws.R | 39 ++++++++++++--- man-roxygen/args-pareto.R | 2 +- man/pareto_diags.Rd | 4 +- man/pareto_khat.Rd | 4 +- man/pareto_smooth.Rd | 4 +- tests/testthat/test-pareto_smooth.R | 15 +++++- 7 files changed, 99 insertions(+), 46 deletions(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 7111821b..1d90067d 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -29,7 +29,7 @@ pareto_khat.default <- function(x, r_eff = NULL, ndraws_tail = NULL, verbose = FALSE, - log_weights = FALSE, + are_log_weights = FALSE, ...) { smoothed <- pareto_smooth.default( x, @@ -39,7 +39,7 @@ pareto_khat.default <- function(x, verbose = verbose, return_k = TRUE, smooth_draws = FALSE, - log_weights = log_weights, + are_log_weights = are_log_weights, ...) return(smoothed$diagnostics) } @@ -126,7 +126,7 @@ pareto_diags.default <- function(x, r_eff = NULL, ndraws_tail = NULL, verbose = FALSE, - log_weights = FALSE, + are_log_weights = FALSE, ...) { smoothed <- pareto_smooth.default( @@ -138,7 +138,7 @@ pareto_diags.default <- function(x, extra_diags = TRUE, verbose = verbose, smooth_draws = FALSE, - log_weights = FALSE, + are_log_weights = FALSE, ...) return(smoothed$diagnostics) @@ -256,15 +256,15 @@ pareto_smooth.default <- function(x, return_k = TRUE, extra_diags = FALSE, verbose = FALSE, - log_weights = FALSE, + are_log_weights = FALSE, ...) { - checkmate::assert_number(ndraws_tail, null.ok = TRUE) - checkmate::assert_number(r_eff, null.ok = TRUE) - checkmate::assert_logical(extra_diags) - checkmate::assert_logical(return_k) - checkmate::assert_logical(verbose) - checkmate::assert_logical(log_weights) + checkmate::expect_numeric(ndraws_tail, null.ok = TRUE) + checkmate::expect_numeric(r_eff, null.ok = TRUE) + extra_diags <- as_one_logical(extra_diags) + return_k <- as_one_logical(return_k) + verbose <- as_one_logical(verbose) + are_log_weights <- as_one_logical(are_log_weights) # check for infinite or na values if (should_return_NA(x)) { @@ -277,8 +277,8 @@ pareto_smooth.default <- function(x, return(out) } - if (log_weights) { - tail = "right" + if (are_log_weights) { + tail <- "right" } tail <- match.arg(tail) @@ -314,7 +314,7 @@ pareto_smooth.default <- function(x, x, ndraws_tail = ndraws_tail, tail = "left", - log_weights = log_weights, + are_log_weights = are_log_weights, ... ) left_k <- smoothed$k @@ -324,13 +324,14 @@ pareto_smooth.default <- function(x, x = smoothed$x, ndraws_tail = ndraws_tail, tail = "right", - log_weights = log_weights, + are_log_weights = are_log_weights, ... ) right_k <- smoothed$k k <- max(left_k, right_k) x <- smoothed$x + } else { smoothed <- .pareto_smooth_tail( @@ -352,10 +353,11 @@ pareto_smooth.default <- function(x, if (verbose) { if (!extra_diags) { - diags_list <- .pareto_smooth_extra_diags(diags_list$khat, length(x)) + diags_list <- c(diags_list, .pareto_smooth_extra_diags(diags_list$khat, length(x))) } pareto_k_diagmsg( - diags = diags_list + diags = diags_list, + are_weights = are_log_weights ) } @@ -375,11 +377,11 @@ pareto_smooth.default <- function(x, ndraws_tail, smooth_draws = TRUE, tail = c("right", "left"), - log_weights = FALSE, + are_log_weights = FALSE, ... ) { - if (log_weights) { + if (are_log_weights) { # shift log values for safe exponentiation x <- x - max(x) } @@ -412,7 +414,7 @@ pareto_smooth.default <- function(x, k <- NA } else { # save time not sorting since x already sorted - if (log_weights) { + if (are_log_weights) { draws_tail <- exp(draws_tail) cutoff <- exp(cutoff) } @@ -422,7 +424,7 @@ pareto_smooth.default <- function(x, if (is.finite(k) && smooth_draws) { p <- (seq_len(ndraws_tail) - 0.5) / ndraws_tail smoothed <- qgeneralized_pareto(p = p, mu = cutoff, k = k, sigma = sigma) - if (log_weights) { + if (are_log_weights) { smoothed <- log(smoothed) } } else { @@ -545,27 +547,38 @@ ps_tail_length <- function(S, r_eff, ...) { #' #' Given S and scalar and k, form a diagnostic message string #' @param diags (numeric) named vector of diagnostic values +#' @param are_weights (logical) are the diagnostics for weights #' @param ... unused #' @return diagnostic message #' @noRd -pareto_k_diagmsg <- function(diags, ...) { +pareto_k_diagmsg <- function(diags, are_weights = FALSE, ...) { khat <- diags$khat min_ss <- diags$min_ss khat_threshold <- diags$khat_threshold convergence_rate <- diags$convergence_rate msg <- NULL - if (khat > 1) { - msg <- paste0(msg,'All estimates are unreliable. If the distribution of ratios is bounded,\n', - 'further draws may improve the estimates, but it is not possible to predict\n', - 'whether any feasible sample size is sufficient.') - } else { - if (khat > khat_threshold) { - msg <- paste0(msg, 'S is too small, and sample size larger than ', round(min_ss, 0), ' is needed for reliable results.\n') + + if (!are_weights) { + + if (khat > 1) { + msg <- paste0(msg,'All estimates are unreliable. If the distribution of draws is bounded,\n', + 'further draws may improve the estimates, but it is not possible to predict\n', + 'whether any feasible sample size is sufficient.') } else { - msg <- paste0(msg, 'To halve the RMSE, approximately ', round(2^(2/convergence_rate),1), ' times bigger S is needed.\n') + if (khat > khat_threshold) { + msg <- paste0(msg, 'S is too small, and sample size larger than ', round(min_ss, 0), ' is needed for reliable results.\n') + } else { + msg <- paste0(msg, 'To halve the RMSE, approximately ', round(2^(2/convergence_rate),1), ' times bigger S is needed.\n') + } + if (khat > 0.7) { + msg <- paste0(msg, 'Bias dominates RMSE, and the variance based MCSE is underestimated.\n') + } } - if (khat > 0.7) { - msg <- paste0(msg, 'Bias dominates RMSE, and the variance based MCSE is underestimated.\n') + + } else { + + if (khat > khat_threshold || khat > 0.7) { + msg <- paste0(msg, 'Pareto khat for weights is high (', round(khat, 1) ,'). This indicates a single or few weights dominate.\n', 'Inference based on weighted draws will be unreliable.\n') } } message(msg) diff --git a/R/weight_draws.R b/R/weight_draws.R index 9270aa14..816075e2 100644 --- a/R/weight_draws.R +++ b/R/weight_draws.R @@ -62,7 +62,12 @@ weight_draws.draws_matrix <- function(x, weights, log = FALSE, pareto_smooth = F log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { - log_weights <- pareto_smooth(log_weights, tail = "right", return_k = FALSE, log = TRUE) + log_weights <- pareto_smooth( + log_weights, + tail = "right", + are_log_weights = TRUE, + extra_diags = TRUE + ) } if (".log_weight" %in% variables(x, reserved = TRUE)) { # overwrite existing weights @@ -77,13 +82,18 @@ weight_draws.draws_matrix <- function(x, weights, log = FALSE, pareto_smooth = F #' @rdname weight_draws #' @export -weight_draws.draws_array <- function(x, weights, log = FALSE, pareto_smooth = FALSE,...) { +weight_draws.draws_array <- function(x, weights, log = FALSE, pareto_smooth = FALSE, ...) { pareto_smooth <- as_one_logical(pareto_smooth) log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { - log_weights <- pareto_smooth(log_weights, tail = "right", return_k = FALSE, log = TRUE) + log_weights <- pareto_smooth( + log_weights, + tail = "right", + are_log_weights = TRUE, + verbose = TRUE + )$x } if (".log_weight" %in% variables(x, reserved = TRUE)) { # overwrite existing weights @@ -104,7 +114,12 @@ weight_draws.draws_df <- function(x, weights, log = FALSE, pareto_smooth = FALSE log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { - log_weights <- pareto_smooth(log_weights, tail = "right", return_k = FALSE, log = TRUE) + log_weights <- pareto_smooth( + log_weights, + tail = "right", + are_log_weights = TRUE, + verbose = TRUE + ) } x$.log_weight <- log_weights x @@ -118,7 +133,13 @@ weight_draws.draws_list <- function(x, weights, log = FALSE, pareto_smooth = FAL log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { - log_weights <- pareto_smooth(log_weights, tail = "right", return_k = FALSE, log = TRUE) + log_weights <- pareto_smooth( + log_weights, + tail = "right", + return_k = TRUE, + are_log_weights = TRUE, + extra_diags = TRUE + ) } niterations <- niterations(x) for (i in seq_len(nchains(x))) { @@ -136,7 +157,13 @@ weight_draws.draws_rvars <- function(x, weights, log = FALSE, pareto_smooth = FA log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { - log_weights <- pareto_smooth(log_weights, tail = "right", return_k = FALSE, log = TRUE) + log_weights <- pareto_smooth( + log_weights, + tail = "right", + return_k = TRUE, + are_log_weights = TRUE, + extra_diags = TRUE + ) } x$.log_weight <- rvar(log_weights) x diff --git a/man-roxygen/args-pareto.R b/man-roxygen/args-pareto.R index db89f30c..8a4d92b4 100644 --- a/man-roxygen/args-pareto.R +++ b/man-roxygen/args-pareto.R @@ -19,6 +19,6 @@ #' @param verbose (logical) Should diagnostic messages be printed? If #' `TRUE`, messages related to Pareto diagnostics will be #' printed. Default is `FALSE`. -#' @param log_weights (logical) Are the draws log weights? Default is +#' @param are_log_weights (logical) Are the draws log weights? Default is #' `FALSE`. If `TRUE` computation will take into account that the #' draws are log weights, and only right tail will be smoothed. diff --git a/man/pareto_diags.Rd b/man/pareto_diags.Rd index 187a1237..9a1d5776 100644 --- a/man/pareto_diags.Rd +++ b/man/pareto_diags.Rd @@ -14,7 +14,7 @@ pareto_diags(x, ...) r_eff = NULL, ndraws_tail = NULL, verbose = FALSE, - log_weights = FALSE, + are_log_weights = FALSE, ... ) @@ -53,7 +53,7 @@ al. (2022)).} \code{TRUE}, messages related to Pareto diagnostics will be printed. Default is \code{FALSE}.} -\item{log_weights}{(logical) Are the draws log weights? Default is +\item{are_log_weights}{(logical) Are the draws log weights? Default is \code{FALSE}. If \code{TRUE} computation will take into account that the draws are log weights, and only right tail will be smoothed.} } diff --git a/man/pareto_khat.Rd b/man/pareto_khat.Rd index 0d844580..a4f91707 100644 --- a/man/pareto_khat.Rd +++ b/man/pareto_khat.Rd @@ -14,7 +14,7 @@ pareto_khat(x, ...) r_eff = NULL, ndraws_tail = NULL, verbose = FALSE, - log_weights = FALSE, + are_log_weights = FALSE, ... ) @@ -53,7 +53,7 @@ al. (2022)).} \code{TRUE}, messages related to Pareto diagnostics will be printed. Default is \code{FALSE}.} -\item{log_weights}{(logical) Are the draws log weights? Default is +\item{are_log_weights}{(logical) Are the draws log weights? Default is \code{FALSE}. If \code{TRUE} computation will take into account that the draws are log weights, and only right tail will be smoothed.} } diff --git a/man/pareto_smooth.Rd b/man/pareto_smooth.Rd index 3e5e444d..c0e6f017 100644 --- a/man/pareto_smooth.Rd +++ b/man/pareto_smooth.Rd @@ -18,7 +18,7 @@ pareto_smooth(x, ...) return_k = TRUE, extra_diags = FALSE, verbose = FALSE, - log_weights = FALSE, + are_log_weights = FALSE, ... ) } @@ -64,7 +64,7 @@ al. (2022)).} \code{TRUE}, messages related to Pareto diagnostics will be printed. Default is \code{FALSE}.} -\item{log_weights}{(logical) Are the draws log weights? Default is +\item{are_log_weights}{(logical) Are the draws log weights? Default is \code{FALSE}. If \code{TRUE} computation will take into account that the draws are log weights, and only right tail will be smoothed.} } diff --git a/tests/testthat/test-pareto_smooth.R b/tests/testthat/test-pareto_smooth.R index dff22ee1..6b67d2b0 100644 --- a/tests/testthat/test-pareto_smooth.R +++ b/tests/testthat/test-pareto_smooth.R @@ -73,7 +73,7 @@ test_that("pareto_khat diagnostics messages are as expected", { diags$khat <- 1.1 expect_message(pareto_k_diagmsg(diags), - paste0('All estimates are unreliable. If the distribution of ratios is bounded,\n', + paste0('All estimates are unreliable. If the distribution of draws is bounded,\n', 'further draws may improve the estimates, but it is not possible to predict\n', 'whether any feasible sample size is sufficient.')) @@ -192,3 +192,16 @@ test_that("pareto_smooth returns x with smoothed tail", { expect_false(isTRUE(all.equal(sort(tau), sort(tau_smoothed)))) }) + +test_that("pareto_smooth works for log_weights", { + w <- c(1:25, 1e3, 1e3, 1e3) + lw <- log(w) + + ps <- pareto_smooth(lw, are_log_weights = TRUE, verbose = FALSE, ndraws_tail = 10) + + # only right tail is smoothed + expect_equal(ps$x[1:15], lw[1:15]) + + expect_true(ps$diagnostics$khat > 0.7) + +}) From ef4cc2a604772f2b7ae7c15853a2d9f64190d140 Mon Sep 17 00:00:00 2001 From: n-kall Date: Thu, 23 Nov 2023 14:58:10 +0200 Subject: [PATCH 12/14] cleanup and test for pareto_smooth in weight_draws --- R/weight_draws.R | 48 ++++++++++-------------------- tests/testthat/test-weight_draws.R | 11 ++++++- 2 files changed, 26 insertions(+), 33 deletions(-) diff --git a/R/weight_draws.R b/R/weight_draws.R index 816075e2..34494820 100644 --- a/R/weight_draws.R +++ b/R/weight_draws.R @@ -62,12 +62,7 @@ weight_draws.draws_matrix <- function(x, weights, log = FALSE, pareto_smooth = F log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { - log_weights <- pareto_smooth( - log_weights, - tail = "right", - are_log_weights = TRUE, - extra_diags = TRUE - ) + log_weights <- pareto_smooth_log_weights(log_weights) } if (".log_weight" %in% variables(x, reserved = TRUE)) { # overwrite existing weights @@ -88,12 +83,7 @@ weight_draws.draws_array <- function(x, weights, log = FALSE, pareto_smooth = FA log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { - log_weights <- pareto_smooth( - log_weights, - tail = "right", - are_log_weights = TRUE, - verbose = TRUE - )$x + log_weights <- pareto_smooth_log_weights(log_weights) } if (".log_weight" %in% variables(x, reserved = TRUE)) { # overwrite existing weights @@ -114,12 +104,7 @@ weight_draws.draws_df <- function(x, weights, log = FALSE, pareto_smooth = FALSE log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { - log_weights <- pareto_smooth( - log_weights, - tail = "right", - are_log_weights = TRUE, - verbose = TRUE - ) + log_weights <- pareto_smooth_log_weights(log_weights) } x$.log_weight <- log_weights x @@ -133,13 +118,7 @@ weight_draws.draws_list <- function(x, weights, log = FALSE, pareto_smooth = FAL log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { - log_weights <- pareto_smooth( - log_weights, - tail = "right", - return_k = TRUE, - are_log_weights = TRUE, - extra_diags = TRUE - ) + log_weights <- pareto_smooth_log_weights(log_weights) } niterations <- niterations(x) for (i in seq_len(nchains(x))) { @@ -157,13 +136,7 @@ weight_draws.draws_rvars <- function(x, weights, log = FALSE, pareto_smooth = FA log <- as_one_logical(log) log_weights <- validate_weights(weights, x, log = log) if (pareto_smooth) { - log_weights <- pareto_smooth( - log_weights, - tail = "right", - return_k = TRUE, - are_log_weights = TRUE, - extra_diags = TRUE - ) + log_weights <- pareto_smooth_log_weights(log_weights) } x$.log_weight <- rvar(log_weights) x @@ -219,3 +192,14 @@ validate_weights <- function(weights, draws, log = FALSE) { } weights } + + +pareto_smooth_log_weights <- function(log_weights) { + pareto_smooth( + log_weights, + tail = "right", + return_k = TRUE, + are_log_weights = TRUE, + extra_diags = TRUE + )$x +} diff --git a/tests/testthat/test-weight_draws.R b/tests/testthat/test-weight_draws.R index 87337439..fb6e6cc3 100644 --- a/tests/testthat/test-weight_draws.R +++ b/tests/testthat/test-weight_draws.R @@ -63,7 +63,6 @@ test_that("weight_draws works on draws_rvars", { expect_equal(weights2, weights) }) - # conversion preserves weights -------------------------------------------- test_that("conversion between formats preserves weights", { @@ -88,3 +87,13 @@ test_that("conversion between formats preserves weights", { expect_equal(as_draws_rvars(draws[[!!type]]), draws$rvars) } }) + +# pareto smoothing ---------------- + +test_that("pareto smoothing smooths weights in weight_draws", { + x <- example_draws() + lw <- sort(log(abs(rt(ndraws(x), 1)))) + weighted <- weight_draws(x, lw, pareto_smooth = FALSE, log = TRUE) + smoothed <- weight_draws(x, lw, pareto_smooth = TRUE, log = TRUE) + expect_false(all(weights(weighted) == weights(smoothed))) +}) From 0f3cf670c1daba91eb5cfc943e1964229d26d34b Mon Sep 17 00:00:00 2001 From: n-kall Date: Thu, 23 Nov 2023 15:26:35 +0200 Subject: [PATCH 13/14] formatting tweak for pareto messages --- R/pareto_smooth.R | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 1d90067d..52ec6bba 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -561,24 +561,24 @@ pareto_k_diagmsg <- function(diags, are_weights = FALSE, ...) { if (!are_weights) { if (khat > 1) { - msg <- paste0(msg,'All estimates are unreliable. If the distribution of draws is bounded,\n', - 'further draws may improve the estimates, but it is not possible to predict\n', - 'whether any feasible sample size is sufficient.') + msg <- paste0(msg, "All estimates are unreliable. If the distribution of draws is bounded,\n", + "further draws may improve the estimates, but it is not possible to predict\n", + "whether any feasible sample size is sufficient.") } else { if (khat > khat_threshold) { - msg <- paste0(msg, 'S is too small, and sample size larger than ', round(min_ss, 0), ' is needed for reliable results.\n') + msg <- paste0(msg, "S is too small, and sample size larger than ", round(min_ss, 0), " is needed for reliable results.\n") } else { - msg <- paste0(msg, 'To halve the RMSE, approximately ', round(2^(2/convergence_rate),1), ' times bigger S is needed.\n') + msg <- paste0(msg, "To halve the RMSE, approximately ", round(2^(2 / convergence_rate), 1), " times bigger S is needed.\n") } if (khat > 0.7) { - msg <- paste0(msg, 'Bias dominates RMSE, and the variance based MCSE is underestimated.\n') + msg <- paste0(msg, "Bias dominates RMSE, and the variance based MCSE is underestimated.\n") } } } else { if (khat > khat_threshold || khat > 0.7) { - msg <- paste0(msg, 'Pareto khat for weights is high (', round(khat, 1) ,'). This indicates a single or few weights dominate.\n', 'Inference based on weighted draws will be unreliable.\n') + msg <- paste0(msg, "Pareto khat for weights is high (", round(khat, 1) ,"). This indicates a single or few weights dominate.\n", "Inference based on weighted draws will be unreliable.\n") } } message(msg) From b062ce4f3062bb38c901806fe55022e3db3b90de Mon Sep 17 00:00:00 2001 From: n-kall Date: Thu, 23 Nov 2023 16:37:43 +0200 Subject: [PATCH 14/14] update news --- NEWS.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 7b9af2b6..a6909143 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,10 +2,11 @@ ### Enhancements +* Add `pareto_smooth` option to `weight_draws`, to Pareto smooth + weights before adding to a draws object. * Matrix multiplication of `rvar`s can now be done with the base matrix multiplication operator (`%*%`) instead of `%**%` in R >= 4.3. - # posterior 1.5.0 ### Enhancements