From 5f38bee761c9df466ca4a54738414dca5084eab4 Mon Sep 17 00:00:00 2001 From: Ozan147 Date: Mon, 22 Aug 2022 08:54:40 +0300 Subject: [PATCH 01/30] added pareto-k diagnostics --- R/gdp.R | 194 ++++++++++++++++++++++++++++++++++++++++++++++ R/pareto_smooth.R | 123 +++++++++++++++++++++++++++++ 2 files changed, 317 insertions(+) create mode 100644 R/gdp.R create mode 100644 R/pareto_smooth.R diff --git a/R/gdp.R b/R/gdp.R new file mode 100644 index 00000000..1442e594 --- /dev/null +++ b/R/gdp.R @@ -0,0 +1,194 @@ +#' The Generalized Pareto Distribution +#' +#' Density, distribution function, quantile function and random generation +#' for the generalized Pareto distribution with location equal to \code{mu}, +#' scale equal to \code{sigma}, and shape equal to \code{k}. +#' +#' @param x,q vector of quantiles. +#' @param p vector of probabilities. +#' @param n number of observations. If \code{length(n) > 1}, the length is taken +#' to be the number required. +#' @param mu scalar location parameter +#' @param sigma scalar, positive scale parameter +#' @param k scalar shape parameter +#' @param log,log.p logical; if TRUE, probabilities p are given as log(p). +#' @param lower.tail logical; if TRUE (default), probabilities are \eqn{P[X \le x]} +#' otherwise, \eqn{P[X > x]}. +#' +#' @name GPD + + +#' @rdname GPD +#' @export +dgpd <- function(x, mu = 0, sigma = 1, k = 0, log = FALSE) { + stopifnot(length(mu) == 1 && length(sigma) == 1 && length(k) == 1) + if (is.na(sigma) || sigma <= 0) { + return(rep(NaN, length(x))) + } + d <- (x - mu) / sigma + ind <- (d > 0) & ((1 + k * d) > 0) + ind[is.na(ind)] <- FALSE + if (k == 0) { + d[ind] <- -d[ind] - log(sigma) + } else { + d[ind] <- log1p(k * d[ind]) * -(1 / k + 1) - log(sigma) + } + d[!ind] <- -Inf + if (!log) { + d <- exp(d) + } + d +} + + +#' @rdname GPD +#' @export +pgpd <- function(q, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, log.p = FALSE) { + stopifnot(length(mu) == 1 && length(sigma) == 1 && length(k) == 1) + if (is.na(sigma) || sigma <= 0) { + return(rep(NaN, length(q))) + } + q <- pmax(q - mu, 0) / sigma + if (k == 0) { + p <- 1 - exp(-q) + } else { + p <- -expm1(log(pmax(1 + k * q, 0)) * -(1 / k)) + } + if (!lower.tail) { + p <- 1 - p + } + if (log.p) { + p <- log(p) + } + p +} + +#' @rdname GPD +#' @export +qgpd <- function(p, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, log.p = FALSE) { + stopifnot(length(mu) == 1 && length(sigma) == 1 && length(k) == 1) + if (is.na(sigma) || sigma <= 0) { + return(rep(NaN, length(p))) + } + if (log.p) { + p <- exp(p) + } + if (!lower.tail) { + p <- 1 - p + } + if (k == 0) { + q <- mu - sigma * log1p(-p) + } else { + q <- mu + sigma * expm1(-k * log1p(-p)) / k + } + q +} + +#' @rdname GPD +#' @export +rgpd <- function(n, mu = 0, sigma = 1, k = 0) { + stopifnot( + length(n) == 1 && length(mu) == 1 && length(sigma) == 1 && length(k) == 1 + ) + if (is.na(sigma) || sigma <= 0) { + return(rep(NaN, n)) + } + if (k == 0) { + r <- mu + sigma * rexp(n) + } else { + r <- mu + sigma * expm1(-k * log(runif(n))) / k + } + r +} + + +#' Estimate parameters of the Generalized Pareto distribution +#' +#' Given a sample \eqn{x}, Estimate the parameters \eqn{k} and \eqn{\sigma} of +#' the generalized Pareto distribution (GPD), assuming the location parameter is +#' 0. By default the fit uses a prior for \eqn{k}, which will stabilize +#' estimates for very small sample sizes (and low effective sample sizes in the +#' case of MCMC samples). The weakly informative prior is a Gaussian prior +#' centered at 0.5. +#' +#' @export +#' @param x A numeric vector. The sample from which to estimate the parameters. +#' @param wip Logical indicating whether to adjust \eqn{k} based on a weakly +#' informative Gaussian prior centered on 0.5. Defaults to `TRUE`. +#' @param min_grid_pts The minimum number of grid points used in the fitting +#' algorithm. The actual number used is `min_grid_pts + floor(sqrt(length(x)))`. +#' @param sort_x If `TRUE` (the default), the first step in the fitting +#' algorithm is to sort the elements of `x`. If `x` is already +#' sorted in ascending order then `sort_x` can be set to `FALSE` to +#' skip the initial sorting step. +#' @return A named list with components `k` and `sigma`. +#' +#' @details Here the parameter \eqn{k} is the negative of \eqn{k} in Zhang & +#' Stephens (2009). +#' +#' @seealso [psis()], [pareto-k-diagnostic] +#' +#' @references +#' Zhang, J., and Stephens, M. A. (2009). A new and efficient estimation method +#' for the generalized Pareto distribution. *Technometrics* **51**, 316-325. +#' +gpdfit <- function(x, wip = TRUE, min_grid_pts = 30, sort_x = TRUE) { + # see section 4 of Zhang and Stephens (2009) + if (sort_x) { + x <- sort.int(x) + } + N <- length(x) + prior <- 3 + M <- min_grid_pts + floor(sqrt(N)) + jj <- seq_len(M) + xstar <- x[floor(N / 4 + 0.5)] # first quartile of sample + theta <- 1 / x[N] + (1 - sqrt(M / (jj - 0.5))) / prior / xstar + k <- matrixStats::rowMeans2(log1p(-theta %o% x)) + l_theta <- N * (log(-theta / k) - k - 1) # profile log-lik + w_theta <- exp(l_theta - matrixStats::logSumExp(l_theta)) # normalize + theta_hat <- sum(theta * w_theta) + k_hat <- mean.default(log1p(-theta_hat * x)) + sigma_hat <- -k_hat / theta_hat + + # adjust k_hat based on weakly informative prior, Gaussian centered on 0.5. + # this stabilizes estimates for very small Monte Carlo sample sizes and low neff + if (wip) { + k_hat <- (k_hat * N + 0.5 * 10) / (N + 10) + } + + if (is.na(k_hat)) { + k_hat <- Inf + sigma_hat <- NaN + } + + list("k" = k_hat, "sigma" = sigma_hat) +} + +# Below is a separate function for the marginal posterior. +# We can extend to have computationally more intensive but more accurate methods. +gpdpost <- function(x, min_grid_pts = 30, sort_x = TRUE) { + + if (sort_x) { + x <- sort.int(x) + } + N <- length(x) + prior <- 3 + M <- min_grid_pts + floor(sqrt(N)) + jj <- seq_len(M) + xstar <- x[floor(N / 4 + 0.5)] # first quartile of sample + theta <- 1 / x[N] + (1 - sqrt(M / (jj - 0.5))) / prior / xstar + k <- matrixStats::rowMeans2(log1p(-theta %o% x)) + l_theta <- N * (log(-theta / k) - k - 1) # profile log-lik + w_theta <- exp(l_theta - matrixStats::logSumExp(l_theta)) # normalize + + # quadrature weights for k are same as for theta + k_w <- w_theta + # quadrature weights are just the normalized likelihoods + # we get the unnormalized posterior by multiplying these by the prior + k_d <- k_w * dgpd(-theta, mu = -1 / x[N], sigma = 1 / prior / xstar, k = 0.5) + # normalize using the trapezoidal rule + Z <- sum((k_d[-M] + k_d[-1]) * (k[-M] - k[-1])) / 2 + k_d <- k_d / Z + + list("k" = k, "k_w" = k_w, "k_d" = k_d) +} diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R new file mode 100644 index 00000000..340b3108 --- /dev/null +++ b/R/pareto_smooth.R @@ -0,0 +1,123 @@ +#' Given Pareto-k computes the minimum sample size for reliable Pareto +#' smoothed estimate (to have small probability of large error) +ps_min_ss <- function(k) { + # Eq (11) in PSIS paper + 10^(1 / (1 - max(0, k))) +} + +#' Given sample size S computes khat threshold for reliable Pareto +#' smoothed estimate (to have small probability of large error) +ps_khat_threshold <- function(S) { + 1 - 1 / log10(S) +} + +#' Given S and scalar or array of k's, compute the relative +#' convergence rate of PSIS estimate RMSE +ps_convergence_rate <- function(k, S) { + # allow array of k's + rate <- numeric(length(k)) + # k<0 bounded distribution + rate[k < 0] <- 1 + # k>0 non-finite mean + rate[k > 1] <- 0 + # limit value at k=1/2 + rate[k == 0.5] <- 1 - 1 / log(S) + # smooth approximation for the rest (see Appendix of PSIS paper) + ki <- (k > 0 & k < 1 & k != 0.5) + kk <- k[ki] + rate[ki] <- pmax( + 0, + (2 * (kk - 1) * S^(2 * kk + 1) + (1 - 2 * kk) * S^(2 * kk) + S^2) / + ((S - 1) * (S - S^(2 * kk))) + ) + rate +} + +#' Given S and scalar and k, form a diagnostic message string +pareto_k_diagmsg <- function(k, S) { + msg <- paste0('With k=', round(k,2), ' and S=', round(S,0), ':\n') + if (k > 1) { + msg <- paste0(msg,'All estimates are unreliable. If the distribution of ratios is bounded,\nfurther draws may improve the estimates, but it is not possible to predict\nwhether any feasible sample size is sufficient.') + } else { + if (k > ps_khat_threshold(S)) { + msg <- paste0(msg, 'S is too small, and sample size larger than ', round(ps_min_ss(k),0), ' is neeeded for reliable results.\n') + } else { + msg <- paste0(msg, 'To half the RMSE, approximately ', round(2^(2/ps_convergence_rate(k,S)),1), ' times bigger S is needed.\n') + } + } + if (k > 0.7 & k < 1) { + msg <- paste0(msg, 'Bias dominates RMSE, and the variance based MCSE is underestimated.\n') + } + message(msg) +} + +# To do ---- +# - use r_eff (scalar, vector, or NULL to compute r_eff from draws) +# - support "left" and both "tails" + +# Question: +# - if draws is a draws matrix or similar object with multiple parameters +# and r_eff is a vector, then how to structure the output? +pareto_khat <- function(draws, + ndraws_tail = ceiling(3 * sqrt(length(draws))), + tail = c("right", "left", "both"), + r_eff = 1, + diagmsg = TRUE) { + + tail <- match.arg(tail) + S <- length(draws) + + if (ndraws_tail >= 5) { + ord <- sort.int(draws, index.return = TRUE) + tail_ids <- seq(S - ndraws_tail + 1, S) + draws_tail <- ord$x[tail_ids] + if (abs(max(draws_tail) - min(draws_tail)) < .Machine$double.eps / 100) { + warning( + "Can't fit generalized Pareto distribution ", + "because all tail values are the same.", + call. = FALSE + ) + } else { + cutoff <- ord$x[min(tail_ids) - 1] # largest value smaller than tail values + # save time not sorting since x already sorted + fit <- gpdfit(draws_tail - cutoff, sort_x = FALSE) + k <- fit$k + sigma <- fit$sigma + if (is.finite(k)) { + p <- (seq_len(ndraws_tail) - 0.5) / ndraws_tail + tail <- qgpd(p = p, mu = cutoff, k = k, sigma = sigma) + } else { + tail <- draws_tail + } + draws[ord$ix[tail_ids]] <- tail + } + } else { + warning( + "Can't fit generalized Pareto distribution ", + "because ndraws_tail is less than 5.", + call. = FALSE + ) + } + + # truncate at max of raw draws + draws[draws > ord$x[S]] <- ord$x[S] + + min_ss <- ps_min_ss(k) + + khat_threshold <- ps_khat_threshold(S) + + convergence_rate <- ps_convergence_rate(k, S) + + if (diagmsg) { + pareto_k_diagmsg(k, S) + } + + list( + "k" = k, + "sigma" = sigma, + "min_ss" = min_ss, + "khat_threshold" = khat_threshold, + "convergence_rate" = convergence_rate + ) + +} From 9ca1c1a03ffc7c5686b13a6a9817345e7dec010f Mon Sep 17 00:00:00 2001 From: n-kall Date: Thu, 23 Mar 2023 13:27:53 +0200 Subject: [PATCH 02/30] Begin cleanup of Pareto-smoothing functions --- R/{gdp.R => gpd.R} | 0 R/pareto_smooth.R | 85 ++++++++++++++++++++++++++++++++++++---------- 2 files changed, 68 insertions(+), 17 deletions(-) rename R/{gdp.R => gpd.R} (100%) diff --git a/R/gdp.R b/R/gpd.R similarity index 100% rename from R/gdp.R rename to R/gpd.R diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 340b3108..e36d3146 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -1,18 +1,35 @@ +##' Pareto-smoothing minimum sample-size +##' #' Given Pareto-k computes the minimum sample size for reliable Pareto #' smoothed estimate (to have small probability of large error) +##' @param k +##' @return minimum sample size ps_min_ss <- function(k) { # Eq (11) in PSIS paper 10^(1 / (1 - max(0, k))) } -#' Given sample size S computes khat threshold for reliable Pareto -#' smoothed estimate (to have small probability of large error) + + +##' Pareto-smoothing k-hat threshold +##' +##' Given sample size S computes khat threshold for reliable Pareto +##' smoothed estimate (to have small probability of large error) +##' @param S sample size +##' @return threshold ps_khat_threshold <- function(S) { 1 - 1 / log10(S) } -#' Given S and scalar or array of k's, compute the relative -#' convergence rate of PSIS estimate RMSE + + +##' Pareto-smoothing convergence rate +##' +##' Given S and scalar or array of k's, compute the relative +##' convergence rate of PSIS estimate RMSE +##' @param k pareto-k values +##' @param S sample size +##' @return convergence rate ps_convergence_rate <- function(k, S) { # allow array of k's rate <- numeric(length(k)) @@ -33,7 +50,13 @@ ps_convergence_rate <- function(k, S) { rate } -#' Given S and scalar and k, form a diagnostic message string + +##' Pareto-k diagnostic message +##' +##' Given S and scalar and k, form a diagnostic message string +##' @param k pareto-k values +##' @param S sample size +##' @return diagnostic message pareto_k_diagmsg <- function(k, S) { msg <- paste0('With k=', round(k,2), ' and S=', round(S,0), ':\n') if (k > 1) { @@ -42,10 +65,10 @@ pareto_k_diagmsg <- function(k, S) { if (k > ps_khat_threshold(S)) { msg <- paste0(msg, 'S is too small, and sample size larger than ', round(ps_min_ss(k),0), ' is neeeded for reliable results.\n') } else { - msg <- paste0(msg, 'To half the RMSE, approximately ', round(2^(2/ps_convergence_rate(k,S)),1), ' times bigger S is needed.\n') + msg <- paste0(msg, 'To halve the RMSE, approximately ', round(2^(2/ps_convergence_rate(k,S)),1), ' times bigger S is needed.\n') } } - if (k > 0.7 & k < 1) { + if (k > 0.7 && k < 1) { msg <- paste0(msg, 'Bias dominates RMSE, and the variance based MCSE is underestimated.\n') } message(msg) @@ -58,17 +81,46 @@ pareto_k_diagmsg <- function(k, S) { # Question: # - if draws is a draws matrix or similar object with multiple parameters # and r_eff is a vector, then how to structure the output? -pareto_khat <- function(draws, - ndraws_tail = ceiling(3 * sqrt(length(draws))), - tail = c("right", "left", "both"), - r_eff = 1, - diagmsg = TRUE) { + +##' Pareto-khat +##' +##' Calculate pareto-khat value given draws +##' @param draws +##' @param ndraws_tail number of draws for the tail. Default is max(ceiling(3 * sqrt(length(draws))), S / 5) +##' @param tail which tail +##' @param r_eff relative effective. Default is "auto" +##' @param diagmsg display diagnostic message? +##' @return list +pareto_khat <- function(x, ...) { + UseMethod("pareto_khat") +} + +pareto_khat.default <- function(x, ...) { + x <- as_draws(x) + pareto_khat(x, ...) +} + + +pareto_khat.draws <- function(x, + ndraws_tail = "default", + tail = c("right", "left", "both"), + r_eff = NULL, + diagmsg = TRUE) { + + if (is.null(r_eff) { + r_eff <- 1 # TODO: calculate this + } + + S <- ndraws(x) + + if (ndraws_tail == "default") { + ndraws_tail <- max(ceiling(3 * sqrt(length(x))), S / 5) + } tail <- match.arg(tail) - S <- length(draws) if (ndraws_tail >= 5) { - ord <- sort.int(draws, index.return = TRUE) + ord <- sort.int(x, index.return = TRUE) tail_ids <- seq(S - ndraws_tail + 1, S) draws_tail <- ord$x[tail_ids] if (abs(max(draws_tail) - min(draws_tail)) < .Machine$double.eps / 100) { @@ -89,7 +141,7 @@ pareto_khat <- function(draws, } else { tail <- draws_tail } - draws[ord$ix[tail_ids]] <- tail + x[ord$ix[tail_ids]] <- tail } } else { warning( @@ -100,7 +152,7 @@ pareto_khat <- function(draws, } # truncate at max of raw draws - draws[draws > ord$x[S]] <- ord$x[S] + x[x > ord$x[S]] <- ord$x[S] min_ss <- ps_min_ss(k) @@ -119,5 +171,4 @@ pareto_khat <- function(draws, "khat_threshold" = khat_threshold, "convergence_rate" = convergence_rate ) - } From 73380a9bf02ec49852ec34f523fc66eff23110eb Mon Sep 17 00:00:00 2001 From: n-kall Date: Thu, 23 Mar 2023 13:45:06 +0200 Subject: [PATCH 03/30] Minor changes to pareto-smooth functions --- R/pareto_smooth.R | 34 ++++++++++++++++++---------------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index e36d3146..8b221f31 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -74,27 +74,29 @@ pareto_k_diagmsg <- function(k, S) { message(msg) } -# To do ---- -# - use r_eff (scalar, vector, or NULL to compute r_eff from draws) -# - support "left" and both "tails" - -# Question: -# - if draws is a draws matrix or similar object with multiple parameters -# and r_eff is a vector, then how to structure the output? - ##' Pareto-khat ##' ##' Calculate pareto-khat value given draws -##' @param draws -##' @param ndraws_tail number of draws for the tail. Default is max(ceiling(3 * sqrt(length(draws))), S / 5) +##' @template args-method-x +##' @param ndraws_tail (numeric) number of draws for the tail. Default is max(ceiling(3 * sqrt(length(draws))), S / 5) ##' @param tail which tail ##' @param r_eff relative effective. Default is "auto" -##' @param diagmsg display diagnostic message? -##' @return list +##' @param verbose (logical) Should a diagnostic message be displayed? Default is `TRUE`. +##' @template args-methods-dots +##' @return List of Pareto-smoothing diagnostics pareto_khat <- function(x, ...) { UseMethod("pareto_khat") } +# To do ---- +# - use r_eff (scalar, vector, or NULL to compute r_eff from draws) +# - support "left" and both "tails" + +# Question: +# - if draws is a draws matrix or similar object with multiple parameters +# and r_eff is a vector, then how to structure the output? + + pareto_khat.default <- function(x, ...) { x <- as_draws(x) pareto_khat(x, ...) @@ -105,13 +107,13 @@ pareto_khat.draws <- function(x, ndraws_tail = "default", tail = c("right", "left", "both"), r_eff = NULL, - diagmsg = TRUE) { + verbose = TRUE) { - if (is.null(r_eff) { + if (is.null(r_eff)) { r_eff <- 1 # TODO: calculate this } - S <- ndraws(x) + S <- length(x) if (ndraws_tail == "default") { ndraws_tail <- max(ceiling(3 * sqrt(length(x))), S / 5) @@ -160,7 +162,7 @@ pareto_khat.draws <- function(x, convergence_rate <- ps_convergence_rate(k, S) - if (diagmsg) { + if (verbose) { pareto_k_diagmsg(k, S) } From 11c40687ee15f3852a7e983302dbda3e4e0af090 Mon Sep 17 00:00:00 2001 From: n-kall Date: Fri, 31 Mar 2023 10:59:05 +0300 Subject: [PATCH 04/30] add left tail pareto-k calculation --- NAMESPACE | 5 ++ R/pareto_smooth.R | 120 ++++++++++++++++++++++++++++++++++------------ 2 files changed, 95 insertions(+), 30 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 709acd43..299d9327 100755 --- a/NAMESPACE +++ b/NAMESPACE @@ -319,6 +319,7 @@ export(chain_ids) export(default_convergence_measures) export(default_mcse_measures) export(default_summary_measures) +export(dgpd) export(diag) export(draw_ids) export(draws_array) @@ -339,6 +340,7 @@ export(example_draws) export(extract_variable) export(extract_variable_matrix) export(for_each_draw) +export(gpdfit) export(is_draws) export(is_draws_array) export(is_draws_df) @@ -359,6 +361,8 @@ export(ndraws) export(niterations) export(nvariables) export(order_draws) +export(pgpd) +export(qgpd) export(quantile2) export(r_scale) export(rdo) @@ -367,6 +371,7 @@ export(repair_draws) export(resample_draws) export(reserved_variables) export(rfun) +export(rgpd) export(rhat) export(rhat_basic) export(rstar) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 8b221f31..c7f56df5 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -3,8 +3,9 @@ #' Given Pareto-k computes the minimum sample size for reliable Pareto #' smoothed estimate (to have small probability of large error) ##' @param k +##' @param ... unused ##' @return minimum sample size -ps_min_ss <- function(k) { +ps_min_ss <- function(k, ...) { # Eq (11) in PSIS paper 10^(1 / (1 - max(0, k))) } @@ -14,10 +15,12 @@ ps_min_ss <- function(k) { ##' Pareto-smoothing k-hat threshold ##' ##' Given sample size S computes khat threshold for reliable Pareto -##' smoothed estimate (to have small probability of large error) +##' smoothed estimate (to have small probability of large error). See +##' section 3.2.4, equation (13). ##' @param S sample size +##' @param ... unused ##' @return threshold -ps_khat_threshold <- function(S) { +ps_khat_threshold <- function(S, ...) { 1 - 1 / log10(S) } @@ -29,8 +32,9 @@ ps_khat_threshold <- function(S) { ##' convergence rate of PSIS estimate RMSE ##' @param k pareto-k values ##' @param S sample size +##' @param ... unused ##' @return convergence rate -ps_convergence_rate <- function(k, S) { +ps_convergence_rate <- function(k, S, ...) { # allow array of k's rate <- numeric(length(k)) # k<0 bounded distribution @@ -56,8 +60,9 @@ ps_convergence_rate <- function(k, S) { ##' Given S and scalar and k, form a diagnostic message string ##' @param k pareto-k values ##' @param S sample size +##' @param ... unused ##' @return diagnostic message -pareto_k_diagmsg <- function(k, S) { +pareto_k_diagmsg <- function(k, S, ...) { msg <- paste0('With k=', round(k,2), ' and S=', round(S,0), ':\n') if (k > 1) { msg <- paste0(msg,'All estimates are unreliable. If the distribution of ratios is bounded,\nfurther draws may improve the estimates, but it is not possible to predict\nwhether any feasible sample size is sufficient.') @@ -77,8 +82,9 @@ pareto_k_diagmsg <- function(k, S) { ##' Pareto-khat ##' ##' Calculate pareto-khat value given draws -##' @template args-method-x -##' @param ndraws_tail (numeric) number of draws for the tail. Default is max(ceiling(3 * sqrt(length(draws))), S / 5) +##' @template args-methods-x +##' @param x object for which method is defined +##' @param ndraws_tail (numeric) number of draws for the tail. Default is max(ceiling(3 * sqrt(length(draws))), S / 5) (Supplementary H) ##' @param tail which tail ##' @param r_eff relative effective. Default is "auto" ##' @param verbose (logical) Should a diagnostic message be displayed? Default is `TRUE`. @@ -98,16 +104,46 @@ pareto_khat <- function(x, ...) { pareto_khat.default <- function(x, ...) { - x <- as_draws(x) - pareto_khat(x, ...) + # x <- as_draws(x) + pareto_khat.draws(x, ...) } -pareto_khat.draws <- function(x, - ndraws_tail = "default", - tail = c("right", "left", "both"), - r_eff = NULL, - verbose = TRUE) { +pareto_khat.default <- function(x, + ndraws_tail = "default", + tail = c("right", "left"), + r_eff = NULL, + verbose = FALSE, + extra_diags = FALSE) { + + tail <- match.arg(tail) + + # flip sign to do left tail + if (tail == "left") { + original_x <- x + x <- -x + } + + out <- .pareto_khat( + x, + ndraws_tail = ndraws_tail, + r_eff = r_eff, + verbose = verbose, + extra_diags = extra_diags + ) + + + names(out) <- paste(names(out), tail, sep = "_") + out + +} + + +.pareto_khat <- function(x, + ndraws_tail = "default", + r_eff = NULL, + verbose = FALSE, + extra_diags = FALSE) { if (is.null(r_eff)) { r_eff <- 1 # TODO: calculate this @@ -119,8 +155,6 @@ pareto_khat.draws <- function(x, ndraws_tail <- max(ceiling(3 * sqrt(length(x))), S / 5) } - tail <- match.arg(tail) - if (ndraws_tail >= 5) { ord <- sort.int(x, index.return = TRUE) tail_ids <- seq(S - ndraws_tail + 1, S) @@ -139,11 +173,9 @@ pareto_khat.draws <- function(x, sigma <- fit$sigma if (is.finite(k)) { p <- (seq_len(ndraws_tail) - 0.5) / ndraws_tail - tail <- qgpd(p = p, mu = cutoff, k = k, sigma = sigma) - } else { - tail <- draws_tail + draws_tail <- qgpd(p = p, mu = cutoff, k = k, sigma = sigma) } - x[ord$ix[tail_ids]] <- tail + x[ord$ix[tail_ids]] <- draws_tail } } else { warning( @@ -156,21 +188,49 @@ pareto_khat.draws <- function(x, # truncate at max of raw draws x[x > ord$x[S]] <- ord$x[S] - min_ss <- ps_min_ss(k) + other_diags <- list() + + if (extra_diags) { - khat_threshold <- ps_khat_threshold(S) + min_ss <- ps_min_ss(k) - convergence_rate <- ps_convergence_rate(k, S) + khat_threshold <- ps_khat_threshold(S) + + convergence_rate <- ps_convergence_rate(k, S) + + other_diags <- list( + "sigma" = sigma, + "min_ss" = min_ss, + "khat_threshold" = khat_threshold, + "convergence_rate" = convergence_rate + ) + } if (verbose) { pareto_k_diagmsg(k, S) } - - list( - "k" = k, - "sigma" = sigma, - "min_ss" = min_ss, - "khat_threshold" = khat_threshold, - "convergence_rate" = convergence_rate + c( + list( + "k" = k + ), + other_diags ) } + +pareto_smooth_tail <- function(x, cutoff) { + + len <- length(x) + + # save time not sorting since x already sorted + fit <- gpdfit(x - cutoff, sort_x = FALSE) + k <- fit$k + sigma <- fit$sigma + if (is.finite(k)) { + p <- (seq_len(len) - 0.5) / len + qq <- qgpd(p, k, sigma) + cutoff + tail <- log(qq) + } else { + tail <- x + } + list(tail = tail, k = k) +} From 9290c215db4f029c90e25f21484fa1de4d209965 Mon Sep 17 00:00:00 2001 From: n-kall Date: Fri, 31 Mar 2023 11:02:46 +0300 Subject: [PATCH 05/30] fix pareto_khat.default to not convert to draws --- R/pareto_smooth.R | 7 ------- 1 file changed, 7 deletions(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index c7f56df5..6f35d61e 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -102,13 +102,6 @@ pareto_khat <- function(x, ...) { # - if draws is a draws matrix or similar object with multiple parameters # and r_eff is a vector, then how to structure the output? - -pareto_khat.default <- function(x, ...) { - # x <- as_draws(x) - pareto_khat.draws(x, ...) -} - - pareto_khat.default <- function(x, ndraws_tail = "default", tail = c("right", "left"), From c51bb534f7ab7d0060851268d94f44d5d1f5e4d2 Mon Sep 17 00:00:00 2001 From: n-kall Date: Sat, 1 Apr 2023 07:47:04 +0300 Subject: [PATCH 06/30] add pareto-khat calculation for "both" tails --- R/pareto_smooth.R | 69 ++++++++++++++++++++++++++++++++++++----------- 1 file changed, 53 insertions(+), 16 deletions(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 6f35d61e..b0f21b64 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -55,6 +55,7 @@ ps_convergence_rate <- function(k, S, ...) { } + ##' Pareto-k diagnostic message ##' ##' Given S and scalar and k, form a diagnostic message string @@ -65,7 +66,9 @@ ps_convergence_rate <- function(k, S, ...) { pareto_k_diagmsg <- function(k, S, ...) { msg <- paste0('With k=', round(k,2), ' and S=', round(S,0), ':\n') if (k > 1) { - msg <- paste0(msg,'All estimates are unreliable. If the distribution of ratios is bounded,\nfurther draws may improve the estimates, but it is not possible to predict\nwhether any feasible sample size is sufficient.') + 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 (k > ps_khat_threshold(S)) { msg <- paste0(msg, 'S is too small, and sample size larger than ', round(ps_min_ss(k),0), ' is neeeded for reliable results.\n') @@ -84,10 +87,13 @@ pareto_k_diagmsg <- function(k, S, ...) { ##' Calculate pareto-khat value given draws ##' @template args-methods-x ##' @param x object for which method is defined -##' @param ndraws_tail (numeric) number of draws for the tail. Default is max(ceiling(3 * sqrt(length(draws))), S / 5) (Supplementary H) +##' @param ndraws_tail (numeric) number of draws for the tail. Default +##' is max(ceiling(3 * sqrt(length(draws))), S / 5) (Supplementary +##' H) ##' @param tail which tail ##' @param r_eff relative effective. Default is "auto" -##' @param verbose (logical) Should a diagnostic message be displayed? Default is `TRUE`. +##' @param verbose (logical) Should a diagnostic message be displayed? +##' Default is `TRUE`. ##' @template args-methods-dots ##' @return List of Pareto-smoothing diagnostics pareto_khat <- function(x, ...) { @@ -104,29 +110,60 @@ pareto_khat <- function(x, ...) { pareto_khat.default <- function(x, ndraws_tail = "default", - tail = c("right", "left"), + tail = c("both", "right", "left"), r_eff = NULL, verbose = FALSE, extra_diags = FALSE) { tail <- match.arg(tail) - # flip sign to do left tail - if (tail == "left") { + if (tail == "both") { + # left first original_x <- x - x <- -x - } + x <- -x + + khat_left <- .pareto_khat( + x, + ndraws_tail = ndraws_tail, + r_eff = r_eff, + verbose = verbose, + extra_diags = extra_diags + ) - out <- .pareto_khat( - x, - ndraws_tail = ndraws_tail, - r_eff = r_eff, - verbose = verbose, - extra_diags = extra_diags - ) + khat_right <- .pareto_khat( + original_x, + ndraws_tail = ndraws_tail, + r_eff = r_eff, + verbose = verbose, + extra_diags = extra_diags + ) + + # take max of khats and corresponding diagnostics + if (khat_left$k >= khat_right$k) { + out <- khat_left + } else { + out <- khat_right + } + + } else { + if (tail == "left") { + # flip sign to do left tail + + original_x <- x + x <- -x + } - names(out) <- paste(names(out), tail, sep = "_") + out <- .pareto_khat( + x, + ndraws_tail = ndraws_tail, + r_eff = r_eff, + verbose = verbose, + extra_diags = extra_diags + ) + names(out) <- paste(names(out), tail, sep = "_") + } + out } From 4f958348a54a37a38a686526995bea1b6d84b53f Mon Sep 17 00:00:00 2001 From: n-kall Date: Sat, 1 Apr 2023 09:08:12 +0300 Subject: [PATCH 07/30] improve pareto-k diagnostic message handling --- R/pareto_smooth.R | 98 ++++++++++++++++++++++++----------------------- 1 file changed, 51 insertions(+), 47 deletions(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index b0f21b64..98b0a817 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -63,17 +63,22 @@ ps_convergence_rate <- function(k, S, ...) { ##' @param S sample size ##' @param ... unused ##' @return diagnostic message -pareto_k_diagmsg <- function(k, S, ...) { - msg <- paste0('With k=', round(k,2), ' and S=', round(S,0), ':\n') +pareto_k_diagmsg <- function(diags, ...) { + k <- diags$k + sigma <- diags$sigma + min_ss <- diags$min_ss + khat_threshold <- diags$khat_threshold + convergence_rate <- diags$convergence_rate + msg <- NULL if (k > 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 (k > ps_khat_threshold(S)) { - msg <- paste0(msg, 'S is too small, and sample size larger than ', round(ps_min_ss(k),0), ' is neeeded for reliable results.\n') + if (k > khat_threshold) { + msg <- paste0(msg, 'S is too small, and sample size larger than ', round(min_ss, 0), ' is neeeded for reliable results.\n') } else { - msg <- paste0(msg, 'To halve the RMSE, approximately ', round(2^(2/ps_convergence_rate(k,S)),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 (k > 0.7 && k < 1) { @@ -94,6 +99,7 @@ pareto_k_diagmsg <- function(k, S, ...) { ##' @param r_eff relative effective. Default is "auto" ##' @param verbose (logical) Should a diagnostic message be displayed? ##' Default is `TRUE`. +##' @param extra_diags include extra pareto-k diagnostics? ##' @template args-methods-dots ##' @return List of Pareto-smoothing diagnostics pareto_khat <- function(x, ...) { @@ -113,7 +119,7 @@ pareto_khat.default <- function(x, tail = c("both", "right", "left"), r_eff = NULL, verbose = FALSE, - extra_diags = FALSE) { + extra_diags = FALSE, ...) { tail <- match.arg(tail) @@ -126,7 +132,6 @@ pareto_khat.default <- function(x, x, ndraws_tail = ndraws_tail, r_eff = r_eff, - verbose = verbose, extra_diags = extra_diags ) @@ -134,16 +139,11 @@ pareto_khat.default <- function(x, original_x, ndraws_tail = ndraws_tail, r_eff = r_eff, - verbose = verbose, extra_diags = extra_diags ) # take max of khats and corresponding diagnostics - if (khat_left$k >= khat_right$k) { - out <- khat_left - } else { - out <- khat_right - } + k <- max(khat_right, khat_left) } else { @@ -154,25 +154,55 @@ pareto_khat.default <- function(x, x <- -x } - out <- .pareto_khat( + k <- .pareto_khat( x, ndraws_tail = ndraws_tail, - r_eff = r_eff, - verbose = verbose, - extra_diags = extra_diags + r_eff = r_eff ) - names(out) <- paste(names(out), tail, sep = "_") + + } + + out <- list(pareto_k = k) + + if (extra_diags) { + diags <- .pareto_k_extra_diags(k, length(x)) + out <- c(out, diags) } - + + if (verbose) { + if (!extra_diags) { + diags <- .pareto_k_extra_diags(k, length(x)) + } + pareto_k_diagmsg( + k = k, + diags = diags + ) + } + + names(out) <- paste(names(out), substring(tail, 1, 1), sep = "_") out } +.pareto_k_extra_diags <- function(k, S, ...) { + + min_ss <- ps_min_ss(k) + + khat_threshold <- ps_khat_threshold(S) + + convergence_rate <- ps_convergence_rate(k, S) + + other_diags <- list( +# "sigma" = sigma, + "min_ss" = min_ss, + "khat_threshold" = khat_threshold, + "convergence_rate" = convergence_rate + ) +} .pareto_khat <- function(x, ndraws_tail = "default", r_eff = NULL, - verbose = FALSE, extra_diags = FALSE) { if (is.null(r_eff)) { @@ -218,33 +248,7 @@ pareto_khat.default <- function(x, # truncate at max of raw draws x[x > ord$x[S]] <- ord$x[S] - other_diags <- list() - - if (extra_diags) { - - min_ss <- ps_min_ss(k) - - khat_threshold <- ps_khat_threshold(S) - - convergence_rate <- ps_convergence_rate(k, S) - - other_diags <- list( - "sigma" = sigma, - "min_ss" = min_ss, - "khat_threshold" = khat_threshold, - "convergence_rate" = convergence_rate - ) - } - - if (verbose) { - pareto_k_diagmsg(k, S) - } - c( - list( - "k" = k - ), - other_diags - ) + k } pareto_smooth_tail <- function(x, cutoff) { From a1d1486b375773e1f92fa1c8cc0b157b91511fb4 Mon Sep 17 00:00:00 2001 From: n-kall Date: Tue, 4 Apr 2023 14:06:40 +0300 Subject: [PATCH 08/30] calculate ndraws_tail for pareto_k in accordance with psis paper --- R/pareto_smooth.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 98b0a817..6c569c26 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -179,7 +179,6 @@ pareto_khat.default <- function(x, ) } - names(out) <- paste(names(out), substring(tail, 1, 1), sep = "_") out } @@ -212,7 +211,8 @@ pareto_khat.default <- function(x, S <- length(x) if (ndraws_tail == "default") { - ndraws_tail <- max(ceiling(3 * sqrt(length(x))), S / 5) + # Appendix H in PSIS paper + ndraws_tail <- ifelse(S > 225, ceiling(3 * sqrt(S/r_eff)), S / 5) } if (ndraws_tail >= 5) { From 500fc60cde529e8b3099763b903e54344038f104 Mon Sep 17 00:00:00 2001 From: n-kall Date: Fri, 14 Apr 2023 08:31:08 +0300 Subject: [PATCH 09/30] refactor pareto-smoothing functions --- R/pareto_smooth.R | 428 ++++++++++++++++++++++++++-------------------- 1 file changed, 239 insertions(+), 189 deletions(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 6c569c26..0c00ddb1 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -1,8 +1,242 @@ +##' Pareto khat diagnostic +##' +##' Calculate pareto khat value from pareto smoothing the tail of +##' x +##' +##' @family diagnostics +##' @template args-methods-x +##' @param ndraws_tail (numeric) number of draws for the tail. Default +##' is max(ceiling(3 * sqrt(length(draws))), S / 5) (Supplementary +##' H) +##' @param tail Which tail to smooth? (one of "right", "left", or "both"). +##' @param r_eff relative effeciency. If "mcmc", it will be calculated +##' automatically. +##' @param verbose (logical) Should diagnostic messages be printed? +##' Default is `FALSE`. +##' @param extra_diags (logical) Should extra pareto smoothing +##' diagnostics be included in output? Default is `FALSE`. +##' @template args-methods-dots +##' @return numeric vector of pareto-smoothing diagnostics +##' +##' @examples +##' mu <- extract_variable_matrix(example_draws(), "mu") +##' pareto_khat(mu) +##' +##' d <- as_draws_rvars(example_draws("multi_normal")) +##' pareto_khat(d$Sigma) +##' +##' @export +pareto_khat <- function(x, ...) UseMethod("pareto_khat") + +##' @rdname pareto_khat +##' @export +pareto_khat.default <- function(x, + ...) { + pareto_smooth.default(x, return_smoothed = FALSE, ...) +} + +##' @rdname pareto_smooth +##' @export +pareto_khat.rvar <- function(x, ...) { + summarise_rvar_by_element(x, pareto_khat.default, ...) +} + +##' Pareto smoothing +##' +##' Smooth the tail of the input with pareto-smoothing. +##' +##' @template args-methods-x +##' @param ndraws_tail (numeric) number of draws for the tail. If +##' `NULL` , it will be calculated as max(ceiling(3 * +##' sqrt(length(x))), length(x) / 5) (Supplementary H) +##' @param tail Which tail to smooth? (one of "right", "left", or +##' "both"). +##' @param r_eff relative effeciency. If "mcmc", it will be calculated +##' automatically. +##' @param verbose (logical) Should diagnostic messages be printed? +##' Default is `FALSE`. +##' @param extra_diags (logical) Should extra pareto smoothing +##' diagnostics be included in output? Default is `FALSE`. +##' @template args-methods-dots +##' @return List of smoothed values and numeric pareto-smoothing +##' diagnostics. +##' +##' @examples +##' mu <- extract_variable_matrix(example_draws(), "mu") +##' pareto_smooth(mu) +##' +##' d <- as_draws_rvars(example_draws("multi_normal")) +##' pareto_smooth(d$Sigma) +pareto_smooth <- function(x, ...) UseMethod("pareto_smooth") + +##' @rdname pareto_smooth +##' @export +pareto_smooth.rvar <- function(x, ...) { + summarise_rvar_by_element(x, pareto_smooth, ...) +} + + +##' @rdname pareto_smooth +##' @export +pareto_smooth.default <- function(x, + tail = c("both", "right", "left"), + r_eff = "mcmc", + ndraws_tail = NULL, + extra_diags = FALSE, + verbose = FALSE, + return_smoothed = TRUE, + ...) { + + tail <- match.arg(tail) + S <- length(x) + + # automatically calculate relative efficiency + if (r_eff == "mcmc") { + r_eff <- ess_basic(x) / S + } + + # automatically calculate tail length + if (is.null(ndraws_tail)) { + ndraws_tail <- ps_tail_length(S, r_eff) + } + + if (tail == "both") { + + if (ndraws_tail > S / 2) { + warning("Number of tail draws cannot be more than half ", + "the length of the draws if both tails are fit, ", + "changing to ", S / 2, ".") + ndraws_tail <- S / 2 + } + + # left tail + smoothed <- .pareto_smooth_tail(x, ndraws_tail = ndraws_tail, tail = "left") + left_k <- smoothed$k + + # right tail + smoothed <-.pareto_smooth_tail(x = smoothed$x, ndraws_tail = ndraws_tail, tail = "right") + right_k <- smoothed$k + + k <- max(left_k, right_k) + x <- smoothed$x + } else { + + smoothed <- .pareto_smooth_tail( + x, + ndraws_tail = ndraws_tail, + tail = tail + ) + k <- smoothed$k + x <- smoothed$x + } + + diags <- c(k = k) + + if (extra_diags) { + ext_diags <- .pareto_smooth_extra_diags(k, S) + diags <- c(diags, ext_diags) + } + + if (verbose) { + if (!extra_diags) { + diags <- .pareto_smooth_extra_diags(out$k, length(x)) + } + pareto_k_diagmsg( + diags = diags + ) + } + if (return_smoothed) { + out <- list(x = x, diagnostics = diags) + } else { + out <- diags + } + return(out) +} + +.pareto_smooth_tail <- function(x, + ndraws_tail, + tail = c("right", "left"), + ... + ) { + + 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 + + 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) { + warning( + "Can't fit generalized Pareto distribution ", + "because all tail values are the same.", + call. = FALSE + ) + } else { + # save time not sorting since x already sorted + fit <- gpdfit(draws_tail - cutoff, sort_x = FALSE) + k <- fit$k + sigma <- fit$sigma + if (is.finite(k)) { + p <- (seq_len(ndraws_tail) - 0.5) / ndraws_tail + smoothed <- qgpd(p = p, mu = cutoff, k = k, sigma = sigma) + } + } + } else { + warning( + "Can't fit generalized Pareto distribution ", + "because ndraws_tail is less than 5.", + call. = FALSE + ) + } + + # truncate at max of raw draws + smoothed[smoothed > max_tail] <- max_tail + + x[ord$ix[tail_ids]] <- smoothed + + + if (tail == "left") { + x <- -x + } + + out <- list(x = x, k = k) + + return(out) +} + +.pareto_smooth_extra_diags <- function(k, S, ...) { + + min_ss <- ps_min_ss(k) + + khat_threshold <- ps_khat_threshold(S) + + convergence_rate <- ps_convergence_rate(k, S) + + other_diags <- list( + "min_ss" = min_ss, + "khat_threshold" = khat_threshold, + "convergence_rate" = convergence_rate + ) +} + ##' Pareto-smoothing minimum sample-size ##' #' Given Pareto-k computes the minimum sample size for reliable Pareto #' smoothed estimate (to have small probability of large error) -##' @param k +##' @param k pareto k value ##' @param ... unused ##' @return minimum sample size ps_min_ss <- function(k, ...) { @@ -10,8 +244,6 @@ ps_min_ss <- function(k, ...) { 10^(1 / (1 - max(0, k))) } - - ##' Pareto-smoothing k-hat threshold ##' ##' Given sample size S computes khat threshold for reliable Pareto @@ -24,8 +256,6 @@ ps_khat_threshold <- function(S, ...) { 1 - 1 / log10(S) } - - ##' Pareto-smoothing convergence rate ##' ##' Given S and scalar or array of k's, compute the relative @@ -54,7 +284,10 @@ ps_convergence_rate <- function(k, S, ...) { rate } - +ps_tail_length <- function(S, r_eff, ...) { + # Appendix H in PSIS paper + ifelse(S > 225, ceiling(3 * sqrt(S / r_eff)), S / 5) +} ##' Pareto-k diagnostic message ##' @@ -65,7 +298,6 @@ ps_convergence_rate <- function(k, S, ...) { ##' @return diagnostic message pareto_k_diagmsg <- function(diags, ...) { k <- diags$k - sigma <- diags$sigma min_ss <- diags$min_ss khat_threshold <- diags$khat_threshold convergence_rate <- diags$convergence_rate @@ -86,185 +318,3 @@ pareto_k_diagmsg <- function(diags, ...) { } message(msg) } - -##' Pareto-khat -##' -##' Calculate pareto-khat value given draws -##' @template args-methods-x -##' @param x object for which method is defined -##' @param ndraws_tail (numeric) number of draws for the tail. Default -##' is max(ceiling(3 * sqrt(length(draws))), S / 5) (Supplementary -##' H) -##' @param tail which tail -##' @param r_eff relative effective. Default is "auto" -##' @param verbose (logical) Should a diagnostic message be displayed? -##' Default is `TRUE`. -##' @param extra_diags include extra pareto-k diagnostics? -##' @template args-methods-dots -##' @return List of Pareto-smoothing diagnostics -pareto_khat <- function(x, ...) { - UseMethod("pareto_khat") -} - -# To do ---- -# - use r_eff (scalar, vector, or NULL to compute r_eff from draws) -# - support "left" and both "tails" - -# Question: -# - if draws is a draws matrix or similar object with multiple parameters -# and r_eff is a vector, then how to structure the output? - -pareto_khat.default <- function(x, - ndraws_tail = "default", - tail = c("both", "right", "left"), - r_eff = NULL, - verbose = FALSE, - extra_diags = FALSE, ...) { - - tail <- match.arg(tail) - - if (tail == "both") { - # left first - original_x <- x - x <- -x - - khat_left <- .pareto_khat( - x, - ndraws_tail = ndraws_tail, - r_eff = r_eff, - extra_diags = extra_diags - ) - - khat_right <- .pareto_khat( - original_x, - ndraws_tail = ndraws_tail, - r_eff = r_eff, - extra_diags = extra_diags - ) - - # take max of khats and corresponding diagnostics - k <- max(khat_right, khat_left) - - } else { - - if (tail == "left") { - # flip sign to do left tail - - original_x <- x - x <- -x - } - - k <- .pareto_khat( - x, - ndraws_tail = ndraws_tail, - r_eff = r_eff - ) - - } - - out <- list(pareto_k = k) - - if (extra_diags) { - diags <- .pareto_k_extra_diags(k, length(x)) - out <- c(out, diags) - } - - if (verbose) { - if (!extra_diags) { - diags <- .pareto_k_extra_diags(k, length(x)) - } - pareto_k_diagmsg( - k = k, - diags = diags - ) - } - - out - -} - -.pareto_k_extra_diags <- function(k, S, ...) { - - min_ss <- ps_min_ss(k) - - khat_threshold <- ps_khat_threshold(S) - - convergence_rate <- ps_convergence_rate(k, S) - - other_diags <- list( -# "sigma" = sigma, - "min_ss" = min_ss, - "khat_threshold" = khat_threshold, - "convergence_rate" = convergence_rate - ) -} - -.pareto_khat <- function(x, - ndraws_tail = "default", - r_eff = NULL, - extra_diags = FALSE) { - - if (is.null(r_eff)) { - r_eff <- 1 # TODO: calculate this - } - - S <- length(x) - - if (ndraws_tail == "default") { - # Appendix H in PSIS paper - ndraws_tail <- ifelse(S > 225, ceiling(3 * sqrt(S/r_eff)), S / 5) - } - - if (ndraws_tail >= 5) { - ord <- sort.int(x, index.return = TRUE) - tail_ids <- seq(S - ndraws_tail + 1, S) - draws_tail <- ord$x[tail_ids] - if (abs(max(draws_tail) - min(draws_tail)) < .Machine$double.eps / 100) { - warning( - "Can't fit generalized Pareto distribution ", - "because all tail values are the same.", - call. = FALSE - ) - } else { - cutoff <- ord$x[min(tail_ids) - 1] # largest value smaller than tail values - # save time not sorting since x already sorted - fit <- gpdfit(draws_tail - cutoff, sort_x = FALSE) - k <- fit$k - sigma <- fit$sigma - if (is.finite(k)) { - p <- (seq_len(ndraws_tail) - 0.5) / ndraws_tail - draws_tail <- qgpd(p = p, mu = cutoff, k = k, sigma = sigma) - } - x[ord$ix[tail_ids]] <- draws_tail - } - } else { - warning( - "Can't fit generalized Pareto distribution ", - "because ndraws_tail is less than 5.", - call. = FALSE - ) - } - - # truncate at max of raw draws - x[x > ord$x[S]] <- ord$x[S] - - k -} - -pareto_smooth_tail <- function(x, cutoff) { - - len <- length(x) - - # save time not sorting since x already sorted - fit <- gpdfit(x - cutoff, sort_x = FALSE) - k <- fit$k - sigma <- fit$sigma - if (is.finite(k)) { - p <- (seq_len(len) - 0.5) / len - qq <- qgpd(p, k, sigma) + cutoff - tail <- log(qq) - } else { - tail <- x - } - list(tail = tail, k = k) -} From 76cfe413ded00cc75264e986d73820abad6ba589 Mon Sep 17 00:00:00 2001 From: n-kall Date: Wed, 26 Apr 2023 13:54:54 +0300 Subject: [PATCH 10/30] Improve documentation for pareto smoothing functions --- NAMESPACE | 6 + R/gpd.R | 1 - R/pareto_smooth.R | 190 ++++++++++---------- man-roxygen/args-pareto.R | 23 +++ man-roxygen/ref-vehtari-paretosmooth-2022.R | 4 + man/GPD.Rd | 42 +++++ man/gpdfit.Rd | 41 +++++ man/pareto_k_diagmsg.Rd | 19 ++ man/pareto_khat.Rd | 65 +++++++ man/pareto_smooth.Rd | 78 ++++++++ man/ps_convergence_rate.Rd | 22 +++ man/ps_khat_threshold.Rd | 21 +++ man/ps_min_ss.Rd | 20 +++ 13 files changed, 431 insertions(+), 101 deletions(-) create mode 100644 man-roxygen/args-pareto.R create mode 100644 man-roxygen/ref-vehtari-paretosmooth-2022.R create mode 100644 man/GPD.Rd create mode 100644 man/gpdfit.Rd create mode 100644 man/pareto_k_diagmsg.Rd create mode 100644 man/pareto_khat.Rd create mode 100644 man/pareto_smooth.Rd create mode 100644 man/ps_convergence_rate.Rd create mode 100644 man/ps_khat_threshold.Rd create mode 100644 man/ps_min_ss.Rd diff --git a/NAMESPACE b/NAMESPACE index 299d9327..a996f2eb 100755 --- a/NAMESPACE +++ b/NAMESPACE @@ -199,6 +199,10 @@ S3method(order_draws,draws_list) S3method(order_draws,draws_matrix) S3method(order_draws,draws_rvars) S3method(order_draws,rvar) +S3method(pareto_khat,default) +S3method(pareto_khat,rvar) +S3method(pareto_smooth,default) +S3method(pareto_smooth,rvar) S3method(pillar_shaft,rvar) S3method(print,draws_array) S3method(print,draws_df) @@ -361,6 +365,8 @@ export(ndraws) export(niterations) export(nvariables) export(order_draws) +export(pareto_khat) +export(pareto_smooth) export(pgpd) export(qgpd) export(quantile2) diff --git a/R/gpd.R b/R/gpd.R index 1442e594..a4f0f927 100644 --- a/R/gpd.R +++ b/R/gpd.R @@ -126,7 +126,6 @@ rgpd <- function(n, mu = 0, sigma = 1, k = 0) { #' @details Here the parameter \eqn{k} is the negative of \eqn{k} in Zhang & #' Stephens (2009). #' -#' @seealso [psis()], [pareto-k-diagnostic] #' #' @references #' Zhang, J., and Stephens, M. A. (2009). A new and efficient estimation method diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 0c00ddb1..2cafd230 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -1,86 +1,69 @@ -##' Pareto khat diagnostic -##' -##' Calculate pareto khat value from pareto smoothing the tail of -##' x -##' -##' @family diagnostics -##' @template args-methods-x -##' @param ndraws_tail (numeric) number of draws for the tail. Default -##' is max(ceiling(3 * sqrt(length(draws))), S / 5) (Supplementary -##' H) -##' @param tail Which tail to smooth? (one of "right", "left", or "both"). -##' @param r_eff relative effeciency. If "mcmc", it will be calculated -##' automatically. -##' @param verbose (logical) Should diagnostic messages be printed? -##' Default is `FALSE`. -##' @param extra_diags (logical) Should extra pareto smoothing -##' diagnostics be included in output? Default is `FALSE`. -##' @template args-methods-dots -##' @return numeric vector of pareto-smoothing diagnostics -##' -##' @examples -##' mu <- extract_variable_matrix(example_draws(), "mu") -##' pareto_khat(mu) -##' -##' d <- as_draws_rvars(example_draws("multi_normal")) -##' pareto_khat(d$Sigma) -##' -##' @export +#' Pareto khat diagnostic +#' +#' Estimate Pareto k value from Pareto smoothing the tail(s) of x. For +#' further details see Vehtari et al. (2022). +#' +#' @template args-pareto +#' @template args-methods-dots +#' @return numeric vector of Pareto smoothing diagnostics +#' +#' @examples +#' mu <- extract_variable_matrix(example_draws(), "mu") +#' pareto_khat(mu) +#' +#' d <- as_draws_rvars(example_draws("multi_normal")) +#' pareto_khat(d$Sigma) +#' +#' @export pareto_khat <- function(x, ...) UseMethod("pareto_khat") -##' @rdname pareto_khat -##' @export +#' @rdname pareto_khat +#' @export pareto_khat.default <- function(x, ...) { pareto_smooth.default(x, return_smoothed = FALSE, ...) } -##' @rdname pareto_smooth -##' @export +#' @rdname pareto_khat +#' @export pareto_khat.rvar <- function(x, ...) { summarise_rvar_by_element(x, pareto_khat.default, ...) } -##' Pareto smoothing -##' -##' Smooth the tail of the input with pareto-smoothing. -##' -##' @template args-methods-x -##' @param ndraws_tail (numeric) number of draws for the tail. If -##' `NULL` , it will be calculated as max(ceiling(3 * -##' sqrt(length(x))), length(x) / 5) (Supplementary H) -##' @param tail Which tail to smooth? (one of "right", "left", or -##' "both"). -##' @param r_eff relative effeciency. If "mcmc", it will be calculated -##' automatically. -##' @param verbose (logical) Should diagnostic messages be printed? -##' Default is `FALSE`. -##' @param extra_diags (logical) Should extra pareto smoothing -##' diagnostics be included in output? Default is `FALSE`. -##' @template args-methods-dots -##' @return List of smoothed values and numeric pareto-smoothing -##' diagnostics. -##' -##' @examples -##' mu <- extract_variable_matrix(example_draws(), "mu") -##' pareto_smooth(mu) -##' -##' d <- as_draws_rvars(example_draws("multi_normal")) -##' pareto_smooth(d$Sigma) +#' Pareto smoothing +#' +#' Smooth the tail(s) of x by fitting a generalized Pareto +#' distribution. For further details see Vehtari et al. (2022). +#' +#' @template args-pareto +#' @template args-methods-dots +#' @param return_smoothed (logical) Should the smoothed x be +#' returned? If `TRUE`, returns smoothed x. If `FALSE`, acts the +#' same as `pareto_khat`. Default is `TRUE`. +#' @return List containing vector of smoothed values and vector of +#' pareto smoothing diagnostics. +#' +#' @examples +#' mu <- extract_variable_matrix(example_draws(), "mu") +#' pareto_smooth(mu) +#' +#' d <- as_draws_rvars(example_draws("multi_normal")) +#' pareto_smooth(d$Sigma) +#' @export pareto_smooth <- function(x, ...) UseMethod("pareto_smooth") -##' @rdname pareto_smooth -##' @export +#' @rdname pareto_smooth +#' @export pareto_smooth.rvar <- function(x, ...) { summarise_rvar_by_element(x, pareto_smooth, ...) } -##' @rdname pareto_smooth -##' @export +#' @rdname pareto_smooth +#' @export pareto_smooth.default <- function(x, tail = c("both", "right", "left"), - r_eff = "mcmc", + r_eff = NULL, ndraws_tail = NULL, extra_diags = FALSE, verbose = FALSE, @@ -91,7 +74,7 @@ pareto_smooth.default <- function(x, S <- length(x) # automatically calculate relative efficiency - if (r_eff == "mcmc") { + if (is.null(r_eff)) { r_eff <- ess_basic(x) / S } @@ -110,11 +93,19 @@ pareto_smooth.default <- function(x, } # left tail - smoothed <- .pareto_smooth_tail(x, ndraws_tail = ndraws_tail, tail = "left") + smoothed <- .pareto_smooth_tail( + x, + ndraws_tail = ndraws_tail, + tail = "left" + ) left_k <- smoothed$k # right tail - smoothed <-.pareto_smooth_tail(x = smoothed$x, ndraws_tail = ndraws_tail, tail = "right") + smoothed <-.pareto_smooth_tail( + x = smoothed$x, + ndraws_tail = ndraws_tail, + tail = "right" + ) right_k <- smoothed$k k <- max(left_k, right_k) @@ -145,7 +136,7 @@ pareto_smooth.default <- function(x, diags = diags ) } - if (return_smoothed) { + if (return_smoothed) { out <- list(x = x, diagnostics = diags) } else { out <- diags @@ -225,45 +216,45 @@ pareto_smooth.default <- function(x, convergence_rate <- ps_convergence_rate(k, S) - other_diags <- list( - "min_ss" = min_ss, - "khat_threshold" = khat_threshold, - "convergence_rate" = convergence_rate - ) + other_diags <- list( + "min_ss" = min_ss, + "khat_threshold" = khat_threshold, + "convergence_rate" = convergence_rate + ) } -##' Pareto-smoothing minimum sample-size -##' +#' Pareto-smoothing minimum sample-size +#' #' Given Pareto-k computes the minimum sample size for reliable Pareto #' smoothed estimate (to have small probability of large error) -##' @param k pareto k value -##' @param ... unused -##' @return minimum sample size +#' @param k pareto k value +#' @param ... unused +#' @return minimum sample size ps_min_ss <- function(k, ...) { # Eq (11) in PSIS paper 10^(1 / (1 - max(0, k))) } -##' Pareto-smoothing k-hat threshold -##' -##' Given sample size S computes khat threshold for reliable Pareto -##' smoothed estimate (to have small probability of large error). See -##' section 3.2.4, equation (13). -##' @param S sample size -##' @param ... unused -##' @return threshold +#' Pareto-smoothing k-hat threshold +#' +#' Given sample size S computes khat threshold for reliable Pareto +#' smoothed estimate (to have small probability of large error). See +#' section 3.2.4, equation (13). +#' @param S sample size +#' @param ... unused +#' @return threshold ps_khat_threshold <- function(S, ...) { 1 - 1 / log10(S) } -##' Pareto-smoothing convergence rate -##' -##' Given S and scalar or array of k's, compute the relative -##' convergence rate of PSIS estimate RMSE -##' @param k pareto-k values -##' @param S sample size -##' @param ... unused -##' @return convergence rate +#' Pareto-smoothing convergence rate +#' +#' Given S and scalar or array of k's, compute the relative +#' convergence rate of PSIS estimate RMSE +#' @param k pareto-k values +#' @param S sample size +#' @param ... unused +#' @return convergence rate ps_convergence_rate <- function(k, S, ...) { # allow array of k's rate <- numeric(length(k)) @@ -289,13 +280,12 @@ ps_tail_length <- function(S, r_eff, ...) { ifelse(S > 225, ceiling(3 * sqrt(S / r_eff)), S / 5) } -##' Pareto-k diagnostic message -##' -##' Given S and scalar and k, form a diagnostic message string -##' @param k pareto-k values -##' @param S sample size -##' @param ... unused -##' @return diagnostic message +#' Pareto-k diagnostic message +#' +#' Given S and scalar and k, form a diagnostic message string +#' @param diags (numeric) named vector of diagnostic values +#' @param ... unused +#' @return diagnostic message pareto_k_diagmsg <- function(diags, ...) { k <- diags$k min_ss <- diags$min_ss diff --git a/man-roxygen/args-pareto.R b/man-roxygen/args-pareto.R new file mode 100644 index 00000000..ca6eabc7 --- /dev/null +++ b/man-roxygen/args-pareto.R @@ -0,0 +1,23 @@ +#' @param x (multiple options) One of: +#' - A matrix of draws for a single variable (iterations x chains). See +#' [extract_variable_matrix()]. +#' - An [`rvar`]. +#' @param tail (string) The tail to smooth: +#' * `"right"`: smooth only the right (upper) tail +#' * `"left"`: smooth only the left (lower) tail +#' * `"both"`: smooth both tails and return the maximum k value +#' +#' The default is `"both"`. +#' @param ndraws_tail (numeric) number of draws for the tail. If +#' `ndraws_tail` is not specified, it will be calculated as +#' max(ceiling(3 * sqrt(length(x))), length(x) / 5). +#' @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 smoothing diagnostics will be +#' printed. Default is `FALSE`. +#' @param extra_diags (logical) Should extra Pareto smoothing +#' diagnostics be included in output? If `TRUE`, `min_ss`, +#' `khat_threshold` and `convergence_rate` for the calculated k +#' value will be returned. Default is `FALSE`. diff --git a/man-roxygen/ref-vehtari-paretosmooth-2022.R b/man-roxygen/ref-vehtari-paretosmooth-2022.R new file mode 100644 index 00000000..267ec29a --- /dev/null +++ b/man-roxygen/ref-vehtari-paretosmooth-2022.R @@ -0,0 +1,4 @@ +#' @references +#' Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao and +#' Jonah Gabry (2022). Pareto Smoothed Importance Sampling. +#' arxiv:arXiv:1507.02646 diff --git a/man/GPD.Rd b/man/GPD.Rd new file mode 100644 index 00000000..9f6e5542 --- /dev/null +++ b/man/GPD.Rd @@ -0,0 +1,42 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gpd.R +\name{GPD} +\alias{GPD} +\alias{dgpd} +\alias{pgpd} +\alias{qgpd} +\alias{rgpd} +\title{The Generalized Pareto Distribution} +\usage{ +dgpd(x, mu = 0, sigma = 1, k = 0, log = FALSE) + +pgpd(q, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, log.p = FALSE) + +qgpd(p, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, log.p = FALSE) + +rgpd(n, mu = 0, sigma = 1, k = 0) +} +\arguments{ +\item{x, q}{vector of quantiles.} + +\item{mu}{scalar location parameter} + +\item{sigma}{scalar, positive scale parameter} + +\item{k}{scalar shape parameter} + +\item{log, log.p}{logical; if TRUE, probabilities p are given as log(p).} + +\item{lower.tail}{logical; if TRUE (default), probabilities are \eqn{P[X \le x]} +otherwise, \eqn{P[X > x]}.} + +\item{p}{vector of probabilities.} + +\item{n}{number of observations. If \code{length(n) > 1}, the length is taken +to be the number required.} +} +\description{ +Density, distribution function, quantile function and random generation +for the generalized Pareto distribution with location equal to \code{mu}, +scale equal to \code{sigma}, and shape equal to \code{k}. +} diff --git a/man/gpdfit.Rd b/man/gpdfit.Rd new file mode 100644 index 00000000..a496f4ad --- /dev/null +++ b/man/gpdfit.Rd @@ -0,0 +1,41 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gpd.R +\name{gpdfit} +\alias{gpdfit} +\title{Estimate parameters of the Generalized Pareto distribution} +\usage{ +gpdfit(x, wip = TRUE, min_grid_pts = 30, sort_x = TRUE) +} +\arguments{ +\item{x}{A numeric vector. The sample from which to estimate the parameters.} + +\item{wip}{Logical indicating whether to adjust \eqn{k} based on a weakly +informative Gaussian prior centered on 0.5. Defaults to \code{TRUE}.} + +\item{min_grid_pts}{The minimum number of grid points used in the fitting +algorithm. The actual number used is \code{min_grid_pts + floor(sqrt(length(x)))}.} + +\item{sort_x}{If \code{TRUE} (the default), the first step in the fitting +algorithm is to sort the elements of \code{x}. If \code{x} is already +sorted in ascending order then \code{sort_x} can be set to \code{FALSE} to +skip the initial sorting step.} +} +\value{ +A named list with components \code{k} and \code{sigma}. +} +\description{ +Given a sample \eqn{x}, Estimate the parameters \eqn{k} and \eqn{\sigma} of +the generalized Pareto distribution (GPD), assuming the location parameter is +0. By default the fit uses a prior for \eqn{k}, which will stabilize +estimates for very small sample sizes (and low effective sample sizes in the +case of MCMC samples). The weakly informative prior is a Gaussian prior +centered at 0.5. +} +\details{ +Here the parameter \eqn{k} is the negative of \eqn{k} in Zhang & +Stephens (2009). +} +\references{ +Zhang, J., and Stephens, M. A. (2009). A new and efficient estimation method +for the generalized Pareto distribution. \emph{Technometrics} \strong{51}, 316-325. +} diff --git a/man/pareto_k_diagmsg.Rd b/man/pareto_k_diagmsg.Rd new file mode 100644 index 00000000..cf000caa --- /dev/null +++ b/man/pareto_k_diagmsg.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pareto_smooth.R +\name{pareto_k_diagmsg} +\alias{pareto_k_diagmsg} +\title{Pareto-k diagnostic message} +\usage{ +pareto_k_diagmsg(diags, ...) +} +\arguments{ +\item{diags}{(numeric) named vector of diagnostic values} + +\item{...}{unused} +} +\value{ +diagnostic message +} +\description{ +Given S and scalar and k, form a diagnostic message string +} diff --git a/man/pareto_khat.Rd b/man/pareto_khat.Rd new file mode 100644 index 00000000..397c9cb4 --- /dev/null +++ b/man/pareto_khat.Rd @@ -0,0 +1,65 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pareto_smooth.R +\name{pareto_khat} +\alias{pareto_khat} +\alias{pareto_khat.default} +\alias{pareto_khat.rvar} +\title{Pareto khat diagnostic} +\usage{ +pareto_khat(x, ...) + +\method{pareto_khat}{default}(x, ...) + +\method{pareto_khat}{rvar}(x, ...) +} +\arguments{ +\item{x}{(multiple options) One of: +\itemize{ +\item A matrix of draws for a single variable (iterations x chains). See +\code{\link[=extract_variable_matrix]{extract_variable_matrix()}}. +\item An \code{\link{rvar}}. +}} + +\item{...}{Arguments passed to individual methods (if applicable).} + +\item{tail}{(string) The tail to smooth: +\itemize{ +\item \code{"right"}: smooth only the right (upper) tail +\item \code{"left"}: smooth only the left (lower) tail +\item \code{"both"}: smooth both tails and return the maximum k value +} + +The default is \code{"both"}.} + +\item{ndraws_tail}{(numeric) number of draws for the tail. If +\code{ndraws_tail} is not specified, it will be calculated as +max(ceiling(3 * sqrt(length(x))), length(x) / 5).} + +\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.} + +\item{verbose}{(logical) Should diagnostic messages be printed? If +\code{TRUE}, messages related to Pareto smoothing diagnostics will be +printed. Default is \code{FALSE}.} + +\item{extra_diags}{(logical) Should extra Pareto smoothing +diagnostics be included in output? If \code{TRUE}, \code{min_ss}, +\code{khat_threshold} and \code{convergence_rate} for the calculated k +value will be returned. Default is \code{FALSE}.} +} +\value{ +numeric vector of Pareto smoothing diagnostics +} +\description{ +Estimate Pareto k value from Pareto smoothing the tail(s) of x. For +further details see Vehtari et al. (2022). +} +\examples{ +mu <- extract_variable_matrix(example_draws(), "mu") +pareto_khat(mu) + +d <- as_draws_rvars(example_draws("multi_normal")) +pareto_khat(d$Sigma) + +} diff --git a/man/pareto_smooth.Rd b/man/pareto_smooth.Rd new file mode 100644 index 00000000..869c0479 --- /dev/null +++ b/man/pareto_smooth.Rd @@ -0,0 +1,78 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pareto_smooth.R +\name{pareto_smooth} +\alias{pareto_smooth} +\alias{pareto_smooth.rvar} +\alias{pareto_smooth.default} +\title{Pareto smoothing} +\usage{ +pareto_smooth(x, ...) + +\method{pareto_smooth}{rvar}(x, ...) + +\method{pareto_smooth}{default}( + x, + tail = c("both", "right", "left"), + r_eff = NULL, + ndraws_tail = NULL, + extra_diags = FALSE, + verbose = FALSE, + return_smoothed = TRUE, + ... +) +} +\arguments{ +\item{x}{(multiple options) One of: +\itemize{ +\item A matrix of draws for a single variable (iterations x chains). See +\code{\link[=extract_variable_matrix]{extract_variable_matrix()}}. +\item An \code{\link{rvar}}. +}} + +\item{...}{Arguments passed to individual methods (if applicable).} + +\item{tail}{(string) The tail to smooth: +\itemize{ +\item \code{"right"}: smooth only the right (upper) tail +\item \code{"left"}: smooth only the left (lower) tail +\item \code{"both"}: smooth both tails and return the maximum k value +} + +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.} + +\item{ndraws_tail}{(numeric) number of draws for the tail. If +\code{ndraws_tail} is not specified, it will be calculated as +max(ceiling(3 * sqrt(length(x))), length(x) / 5).} + +\item{extra_diags}{(logical) Should extra Pareto smoothing +diagnostics be included in output? If \code{TRUE}, \code{min_ss}, +\code{khat_threshold} and \code{convergence_rate} for the calculated k +value will be returned. Default is \code{FALSE}.} + +\item{verbose}{(logical) Should diagnostic messages be printed? If +\code{TRUE}, messages related to Pareto smoothing diagnostics will be +printed. Default is \code{FALSE}.} + +\item{return_smoothed}{(logical) Should the smoothed x be +returned? If \code{TRUE}, returns smoothed x. If \code{FALSE}, acts the +same as \code{pareto_khat}. Default is \code{TRUE}.} +} +\value{ +List containing vector of smoothed values and vector of +pareto smoothing diagnostics. +} +\description{ +Smooth the tail(s) of x by fitting a generalized Pareto +distribution. For further details see Vehtari et al. (2022). +} +\examples{ +mu <- extract_variable_matrix(example_draws(), "mu") +pareto_smooth(mu) + +d <- as_draws_rvars(example_draws("multi_normal")) +pareto_smooth(d$Sigma) +} diff --git a/man/ps_convergence_rate.Rd b/man/ps_convergence_rate.Rd new file mode 100644 index 00000000..40948fc4 --- /dev/null +++ b/man/ps_convergence_rate.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pareto_smooth.R +\name{ps_convergence_rate} +\alias{ps_convergence_rate} +\title{Pareto-smoothing convergence rate} +\usage{ +ps_convergence_rate(k, S, ...) +} +\arguments{ +\item{k}{pareto-k values} + +\item{S}{sample size} + +\item{...}{unused} +} +\value{ +convergence rate +} +\description{ +Given S and scalar or array of k's, compute the relative +convergence rate of PSIS estimate RMSE +} diff --git a/man/ps_khat_threshold.Rd b/man/ps_khat_threshold.Rd new file mode 100644 index 00000000..1cdfd459 --- /dev/null +++ b/man/ps_khat_threshold.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pareto_smooth.R +\name{ps_khat_threshold} +\alias{ps_khat_threshold} +\title{Pareto-smoothing k-hat threshold} +\usage{ +ps_khat_threshold(S, ...) +} +\arguments{ +\item{S}{sample size} + +\item{...}{unused} +} +\value{ +threshold +} +\description{ +Given sample size S computes khat threshold for reliable Pareto +smoothed estimate (to have small probability of large error). See +section 3.2.4, equation (13). +} diff --git a/man/ps_min_ss.Rd b/man/ps_min_ss.Rd new file mode 100644 index 00000000..2b48a38e --- /dev/null +++ b/man/ps_min_ss.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pareto_smooth.R +\name{ps_min_ss} +\alias{ps_min_ss} +\title{Pareto-smoothing minimum sample-size} +\usage{ +ps_min_ss(k, ...) +} +\arguments{ +\item{k}{pareto k value} + +\item{...}{unused} +} +\value{ +minimum sample size +} +\description{ +Given Pareto-k computes the minimum sample size for reliable Pareto +smoothed estimate (to have small probability of large error) +} From cae2893c2d68653177cc9182c4f8de039ae78e9b Mon Sep 17 00:00:00 2001 From: n-kall Date: Thu, 27 Apr 2023 10:42:11 +0300 Subject: [PATCH 11/30] fix documentation for ndraws_tail in pareto smoothing functions --- man-roxygen/args-pareto.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/man-roxygen/args-pareto.R b/man-roxygen/args-pareto.R index ca6eabc7..daff8977 100644 --- a/man-roxygen/args-pareto.R +++ b/man-roxygen/args-pareto.R @@ -10,7 +10,8 @@ #' The default is `"both"`. #' @param ndraws_tail (numeric) number of draws for the tail. If #' `ndraws_tail` is not specified, it will be calculated as -#' max(ceiling(3 * sqrt(length(x))), length(x) / 5). +#' ceiling(3 * sqrt(length(x) / r_eff)) if length(x) > 225 and +#' 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. From 30dff2e8d7e8b8b0f7e3a954518610139b1e9b73 Mon Sep 17 00:00:00 2001 From: n-kall Date: Thu, 27 Apr 2023 10:57:20 +0300 Subject: [PATCH 12/30] update documentation for pareto smoothing --- man/pareto_khat.Rd | 3 ++- man/pareto_smooth.Rd | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/man/pareto_khat.Rd b/man/pareto_khat.Rd index 397c9cb4..842bb3d0 100644 --- a/man/pareto_khat.Rd +++ b/man/pareto_khat.Rd @@ -33,7 +33,8 @@ The default is \code{"both"}.} \item{ndraws_tail}{(numeric) number of draws for the tail. If \code{ndraws_tail} is not specified, it will be calculated as -max(ceiling(3 * sqrt(length(x))), length(x) / 5).} +ceiling(3 * sqrt(length(x) / r_eff)) if length(x) > 225 and +length(x) / 5 otherwise (see Appendix H in Vehtari et al. (2022)).} \item{r_eff}{(numeric) relative effective sample size estimate. If \code{r_eff} is omitted, it will be calculated assuming the draws are diff --git a/man/pareto_smooth.Rd b/man/pareto_smooth.Rd index 869c0479..5d4ebd62 100644 --- a/man/pareto_smooth.Rd +++ b/man/pareto_smooth.Rd @@ -46,7 +46,8 @@ from MCMC.} \item{ndraws_tail}{(numeric) number of draws for the tail. If \code{ndraws_tail} is not specified, it will be calculated as -max(ceiling(3 * sqrt(length(x))), length(x) / 5).} +ceiling(3 * sqrt(length(x) / r_eff)) if length(x) > 225 and +length(x) / 5 otherwise (see Appendix H in Vehtari et al. (2022)).} \item{extra_diags}{(logical) Should extra Pareto smoothing diagnostics be included in output? If \code{TRUE}, \code{min_ss}, From 0d1acf97dc7dada8153aac0edabc61421ba67d48 Mon Sep 17 00:00:00 2001 From: n-kall Date: Thu, 27 Apr 2023 10:57:31 +0300 Subject: [PATCH 13/30] return pareto diagnostics as list --- R/pareto_smooth.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 2cafd230..56e0b56c 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -121,7 +121,7 @@ pareto_smooth.default <- function(x, x <- smoothed$x } - diags <- c(k = k) + diags <- list(k = k) if (extra_diags) { ext_diags <- .pareto_smooth_extra_diags(k, S) From 049d8fed5f92c18b2a4632451ac066f50a224718 Mon Sep 17 00:00:00 2001 From: n-kall Date: Wed, 3 May 2023 14:01:34 +0300 Subject: [PATCH 14/30] Add tests and cleanup pareto diagnostics --- R/pareto_smooth.R | 54 +++++---- tests/testthat/test-pareto_smooth.R | 171 ++++++++++++++++++++++++++++ 2 files changed, 206 insertions(+), 19 deletions(-) create mode 100644 tests/testthat/test-pareto_smooth.R diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 56e0b56c..1554bafd 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -13,7 +13,6 @@ #' #' d <- as_draws_rvars(example_draws("multi_normal")) #' pareto_khat(d$Sigma) -#' #' @export pareto_khat <- function(x, ...) UseMethod("pareto_khat") @@ -27,7 +26,10 @@ pareto_khat.default <- function(x, #' @rdname pareto_khat #' @export pareto_khat.rvar <- function(x, ...) { - summarise_rvar_by_element(x, pareto_khat.default, ...) + unlist( + summarise_rvar_by_element_with_chains(x, pareto_khat.default, ...), + recursive = FALSE + ) } #' Pareto smoothing @@ -37,8 +39,8 @@ pareto_khat.rvar <- function(x, ...) { #' #' @template args-pareto #' @template args-methods-dots -#' @param return_smoothed (logical) Should the smoothed x be -#' returned? If `TRUE`, returns smoothed x. If `FALSE`, acts the +#' @param return_smoothed (logical) Should the smoothed x be returned? +#' If `TRUE`, returns x with smoothed tail(s). If `FALSE`, acts the #' same as `pareto_khat`. Default is `TRUE`. #' @return List containing vector of smoothed values and vector of #' pareto smoothing diagnostics. @@ -55,10 +57,9 @@ pareto_smooth <- function(x, ...) UseMethod("pareto_smooth") #' @rdname pareto_smooth #' @export pareto_smooth.rvar <- function(x, ...) { - summarise_rvar_by_element(x, pareto_smooth, ...) + unlist(summarise_rvar_by_element_with_chains(x, pareto_smooth.default, ...), recursive = FALSE) } - #' @rdname pareto_smooth #' @export pareto_smooth.default <- function(x, @@ -70,9 +71,14 @@ pareto_smooth.default <- function(x, return_smoothed = TRUE, ...) { + # check for infinite or na values + if (should_return_NA(x)) { + return(NA_real_) + } + tail <- match.arg(tail) S <- length(x) - + # automatically calculate relative efficiency if (is.null(r_eff)) { r_eff <- ess_basic(x) / S @@ -92,6 +98,12 @@ pareto_smooth.default <- function(x, ndraws_tail <- S / 2 } + if (ndraws_tail < 5) { + warning("Number of tail draws cannot be less than 5. ", + "Changing to ", 5, ".") + ndraws_tail <- 5 + } + # left tail smoothed <- .pareto_smooth_tail( x, @@ -121,7 +133,7 @@ pareto_smooth.default <- function(x, x <- smoothed$x } - diags <- list(k = k) + diags <- list(khat = k) if (extra_diags) { ext_diags <- .pareto_smooth_extra_diags(k, S) @@ -130,7 +142,7 @@ pareto_smooth.default <- function(x, if (verbose) { if (!extra_diags) { - diags <- .pareto_smooth_extra_diags(out$k, length(x)) + diags <- .pareto_smooth_extra_diags(diags$khat, length(x)) } pareto_k_diagmsg( diags = diags @@ -175,6 +187,8 @@ pareto_smooth.default <- function(x, "because all tail values are the same.", call. = FALSE ) + smoothed <- NULL + k <- NA } else { # save time not sorting since x already sorted fit <- gpdfit(draws_tail - cutoff, sort_x = FALSE) @@ -191,13 +205,15 @@ pareto_smooth.default <- function(x, "because ndraws_tail is less than 5.", call. = FALSE ) + smoothed <- NULL + k <- NA } # truncate at max of raw draws - smoothed[smoothed > max_tail] <- max_tail - - x[ord$ix[tail_ids]] <- smoothed - + if (!is.null(smoothed)) { + smoothed[smoothed > max_tail] <- max_tail + x[ord$ix[tail_ids]] <- smoothed + } if (tail == "left") { x <- -x @@ -287,24 +303,24 @@ ps_tail_length <- function(S, r_eff, ...) { #' @param ... unused #' @return diagnostic message pareto_k_diagmsg <- function(diags, ...) { - k <- diags$k + khat <- diags$khat min_ss <- diags$min_ss khat_threshold <- diags$khat_threshold convergence_rate <- diags$convergence_rate msg <- NULL - if (k > 1) { + 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 (k > khat_threshold) { + if (khat > khat_threshold) { msg <- paste0(msg, 'S is too small, and sample size larger than ', round(min_ss, 0), ' is neeeded 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 (k > 0.7 && k < 1) { - 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') + } } message(msg) } diff --git a/tests/testthat/test-pareto_smooth.R b/tests/testthat/test-pareto_smooth.R new file mode 100644 index 00000000..f71c3ce4 --- /dev/null +++ b/tests/testthat/test-pareto_smooth.R @@ -0,0 +1,171 @@ +test_that("pareto_khat returns expected reasonable values", { + tau <- extract_variable_matrix(example_draws(), "tau") + + pk <- pareto_khat(tau) + expect_true(names(pk) == "khat") + +}) + +test_that("pareto_khat handles tail argument", { + + # as tau is bounded (0, Inf) the left pareto k should be lower than + # right + tau <- extract_variable_matrix(example_draws(), "tau") + pkl <- pareto_khat(tau, tail = "left") + pkr <- pareto_khat(tau, tail = "right") + pkb <- pareto_khat(tau) + expect_true(pkl$khat < pkr$khat) + expect_equal(pkr$khat, pkb$khat) +}) + +test_that("pareto_khat handles ndraws_tail argument", { + + tau <- extract_variable_matrix(example_draws(), "tau") + pk10 <- pareto_khat(tau, tail = "right", ndraws_tail = 10) + pk25 <- pareto_khat(tau, tail = "right", ndraws_tail = 25) + expect_true(pk10$khat > pk25$khat) + + expect_warning(pareto_khat(tau, tail = "both", ndraws_tail = 201), + "Number of tail draws cannot be more than half ", + "the length of the draws if both tails are fit, ", + "changing to 200.") + + expect_warning(pareto_khat(tau, tail = "both", ndraws_tail = 4), + "Number of tail draws cannot be less than 5. Changing to 5.") + +}) + + +test_that("pareto_khat handles r_eff argument", { + + tau <- extract_variable_matrix(example_draws(), "tau") + pk1 <- pareto_khat(tau, r_eff = 1) + pk0.6 <- pareto_khat(tau, r_eff = 0.6) + expect_true(pk1$khat < pk0.6$khat) + +}) + + +test_that("pareto_khat diagnostics messages are as expected", { + + diags <- list( + khat = 0.5, + min_ss = 10, + khat_threshold = 0.55, + convergence_rate = 0.99 + ) + + expect_message(pareto_k_diagmsg(diags), + paste0('To halve the RMSE, approximately 4.1 times bigger S is needed.')) + + diags$khat <- 0.6 + + expect_message(pareto_k_diagmsg(diags), + paste0('S is too small, and sample size larger than 10 is neeeded for reliable results.\n')) + + diags$khat <- 0.71 + diags$khat_threshold <- 0.8 + + expect_message(pareto_k_diagmsg(diags), + paste0('To halve the RMSE, approximately 4.1 times bigger S is needed.\n', 'Bias dominates RMSE, and the variance based MCSE is underestimated.\n')) + + + diags$khat <- 1.1 + + expect_message(pareto_k_diagmsg(diags), + paste0('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.')) + +}) + +test_that("pareto_khat extra diags are returned", { + tau <- extract_variable_matrix(example_draws(), "tau") + + pk <- pareto_khat(tau, extra_diag = TRUE) + expect_true(all(names(pk) == c("khat", "min_ss", "khat_threshold", "convergence_rate"))) +}) + + +test_that("pareto_khat diagnostics handles verbose argument", { + tau <- extract_variable_matrix(example_draws(), "tau") + + expect_message(pareto_khat(tau, extra_diag = TRUE, verbose = TRUE)) + expect_message(pareto_khat(tau, extra_diag = FALSE, verbose = TRUE)) + expect_no_message(pareto_khat(tau, extra_diag = TRUE, verbose = FALSE)) + +}) + +test_that("pareto_khat diagnostics accept vectors as input", { + set.seed(1234) + x <- rnorm(1000) + + pk <- pareto_khat(x) + expect_true(pk < 1) +}) + +test_that("pareto_khat diagnostics handle special cases correctly", { + set.seed(1234) + + x <- c(rnorm(50), NA) + expect_true(is.na(pareto_khat(x))) + + x <- c(rnorm(50), Inf) + expect_true(is.na(pareto_khat(x))) + + x <- rep(1, 50) + expect_true(is.na(pareto_khat(x))) + +}) + + +test_that("pareto_khat functions work with matrix with chains", { + + tau_chains <- extract_variable_matrix(example_draws(), "tau") + tau_nochains <- extract_variable(example_draws(), "tau") + + expect_equal(pareto_khat(tau_chains, ndraws_tail = 20), + pareto_khat(tau_nochains, ndraws_tail = 20)) + + ps_chains <- pareto_smooth(tau_chains, ndraws_tail = 20) + ps_nochains <- pareto_smooth(tau_nochains, ndraws_tail = 20) + + expect_equal(as.numeric(ps_chains$x), as.numeric(ps_nochains$x)) + +}) + +test_that("pareto_khat functions work with rvars with and without chains", { + + tau_chains <- extract_variable_matrix(example_draws(), "tau") + tau_nochains <- extract_variable(example_draws(), "tau") + + tau_rvar_chains <- rvar(tau_chains, with_chains = TRUE) + tau_rvar_nochains <- rvar(tau_nochains) + + expect_equal(pareto_khat(tau_chains, ndraws_tail = 20), + pareto_khat(tau_rvar_chains, ndraws_tail = 20)) + + expect_equal(pareto_khat(tau_rvar_chains, ndraws_tail = 20), + pareto_khat(tau_rvar_nochains, ndraws_tail = 20)) + + ps_chains <- pareto_smooth(tau_chains, ndraws_tail = 20) + ps_rvar_chains <- pareto_smooth(tau_rvar_chains, ndraws_tail = 20) + + ps_nochains <- pareto_smooth(tau_nochains, ndraws_tail = 20) + ps_rvar_nochains <- pareto_smooth(tau_rvar_nochains, ndraws_tail = 20) + + expect_equal(as.numeric(ps_rvar_chains$x), as.numeric(ps_chains$x)) + expect_equal(as.numeric(ps_rvar_nochains$x), as.numeric(ps_rvar_chains$x)) + +}) + +test_that("pareto_smooth returns x with smoothed tail", { + tau <- extract_variable_matrix(example_draws(), "tau") + + tau_smoothed <- pareto_smooth(tau, ndraws_tail = 10, tail = "right")$x + + expect_equal(sort(tau)[1:390], sort(tau_smoothed)[1:390]) + + expect_false(isTRUE(all.equal(sort(tau), sort(tau_smoothed)))) + +}) From 3ae4af4d19cf3fe6a9096b6484d7c3bdb1770316 Mon Sep 17 00:00:00 2001 From: n-kall Date: Tue, 9 May 2023 11:36:09 +0300 Subject: [PATCH 15/30] improvements to pareto smooth functions (fix rvar method, simplify) --- R/pareto_smooth.R | 84 ++++++++++++++++++----------- tests/testthat/test-pareto_smooth.R | 5 +- 2 files changed, 57 insertions(+), 32 deletions(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 1554bafd..45b5e48b 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -20,15 +20,18 @@ pareto_khat <- function(x, ...) UseMethod("pareto_khat") #' @export pareto_khat.default <- function(x, ...) { - pareto_smooth.default(x, return_smoothed = FALSE, ...) + out <- pareto_smooth.default(x, ...) + return(out$diagnostics) } #' @rdname pareto_khat #' @export pareto_khat.rvar <- function(x, ...) { - unlist( - summarise_rvar_by_element_with_chains(x, pareto_khat.default, ...), - recursive = FALSE + draws_diags <- summarise_rvar_by_element_with_chains(x, pareto_smooth.default, ...) + dim(draws_diags) <- dim(draws_diags) %||% length(draws_diags) + margins <- seq_along(dim(draws_diags)) + list( + khat = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat) ) } @@ -39,11 +42,9 @@ pareto_khat.rvar <- function(x, ...) { #' #' @template args-pareto #' @template args-methods-dots -#' @param return_smoothed (logical) Should the smoothed x be returned? -#' If `TRUE`, returns x with smoothed tail(s). If `FALSE`, acts the -#' same as `pareto_khat`. Default is `TRUE`. -#' @return List containing vector of smoothed values and vector of -#' pareto smoothing diagnostics. +#' @return List containing vector `x` of smoothed values and list +#' `diagnostics` of Pareto smoothing diagnostics (`khat`, and +#' optionally `min_ss`, `khat_threshold`, and `convergence_rate`). #' #' @examples #' mu <- extract_variable_matrix(example_draws(), "mu") @@ -56,8 +57,29 @@ pareto_smooth <- function(x, ...) UseMethod("pareto_smooth") #' @rdname pareto_smooth #' @export -pareto_smooth.rvar <- function(x, ...) { - unlist(summarise_rvar_by_element_with_chains(x, pareto_smooth.default, ...), recursive = FALSE) +pareto_smooth.rvar <- function(x, extra_diags = FALSE, ...) { + draws_diags <- summarise_rvar_by_element_with_chains(x, pareto_smooth.default, extra_diags = extra_diags, ...) + dim(draws_diags) <- dim(draws_diags) %||% length(draws_diags) + margins <- seq_along(dim(draws_diags)) + + if (extra_diags) { + + diags <- list( + khat = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat), + min_ss = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$min_ss), + khat_threshold = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat_threshold), + convergence_rage = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$convergence_rate) + ) + } else { + diags <- list( + khat = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat) + ) + } + + list( + x = rvar(apply(draws_diags, margins, function(x) x[[1]]$x), nchains = nchains(x)), + diagnostics = diags + ) } #' @rdname pareto_smooth @@ -68,17 +90,16 @@ pareto_smooth.default <- function(x, ndraws_tail = NULL, extra_diags = FALSE, verbose = FALSE, - return_smoothed = TRUE, ...) { # check for infinite or na values if (should_return_NA(x)) { - return(NA_real_) + return(list(x = x, diagnostics = NA_real_)) } tail <- match.arg(tail) S <- length(x) - + # automatically calculate relative efficiency if (is.null(r_eff)) { r_eff <- ess_basic(x) / S @@ -88,9 +109,9 @@ pareto_smooth.default <- function(x, if (is.null(ndraws_tail)) { ndraws_tail <- ps_tail_length(S, r_eff) } - + if (tail == "both") { - + if (ndraws_tail > S / 2) { warning("Number of tail draws cannot be more than half ", "the length of the draws if both tails are fit, ", @@ -123,7 +144,7 @@ pareto_smooth.default <- function(x, k <- max(left_k, right_k) x <- smoothed$x } else { - + smoothed <- .pareto_smooth_tail( x, ndraws_tail = ndraws_tail, @@ -148,14 +169,13 @@ pareto_smooth.default <- function(x, diags = diags ) } - if (return_smoothed) { - out <- list(x = x, diagnostics = diags) - } else { - out <- diags - } + + out <- list(x = x, diagnostics = diags) + return(out) } +# internal function to pareto smooth the tail of a vector .pareto_smooth_tail <- function(x, ndraws_tail, tail = c("right", "left"), @@ -166,8 +186,8 @@ pareto_smooth.default <- function(x, S <- length(x) tail_ids <- seq(S - ndraws_tail + 1, S) - - + + if (tail == "left") { x <- -x } @@ -175,13 +195,14 @@ pareto_smooth.default <- function(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 - + 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) { + print(paste(max_tail, min_tail)) warning( "Can't fit generalized Pareto distribution ", "because all tail values are the same.", @@ -214,16 +235,18 @@ pareto_smooth.default <- function(x, smoothed[smoothed > max_tail] <- max_tail x[ord$ix[tail_ids]] <- smoothed } - + if (tail == "left") { x <- -x } - + out <- list(x = x, k = k) return(out) } +# internal function to calculate the extra diagnostics or a given +# pareto k and sample size S .pareto_smooth_extra_diags <- function(k, S, ...) { min_ss <- ps_min_ss(k) @@ -233,9 +256,9 @@ pareto_smooth.default <- function(x, convergence_rate <- ps_convergence_rate(k, S) other_diags <- list( - "min_ss" = min_ss, - "khat_threshold" = khat_threshold, - "convergence_rate" = convergence_rate + min_ss = min_ss, + khat_threshold = khat_threshold, + convergence_rate = convergence_rate ) } @@ -323,4 +346,5 @@ pareto_k_diagmsg <- function(diags, ...) { } } message(msg) + invisible(diags) } diff --git a/tests/testthat/test-pareto_smooth.R b/tests/testthat/test-pareto_smooth.R index f71c3ce4..e2bddb60 100644 --- a/tests/testthat/test-pareto_smooth.R +++ b/tests/testthat/test-pareto_smooth.R @@ -154,8 +154,9 @@ test_that("pareto_khat functions work with rvars with and without chains", { ps_nochains <- pareto_smooth(tau_nochains, ndraws_tail = 20) ps_rvar_nochains <- pareto_smooth(tau_rvar_nochains, ndraws_tail = 20) - expect_equal(as.numeric(ps_rvar_chains$x), as.numeric(ps_chains$x)) - expect_equal(as.numeric(ps_rvar_nochains$x), as.numeric(ps_rvar_chains$x)) + expect_equal(ps_rvar_chains$x, rvar(ps_chains$x, with_chains = TRUE)) + + expect_equal(ps_rvar_nochains$x, rvar(ps_nochains$x)) }) From 6e232c33d68df5787ab1554ce5b0f6130200ef66 Mon Sep 17 00:00:00 2001 From: n-kall Date: Tue, 9 May 2023 11:39:40 +0300 Subject: [PATCH 16/30] cleanup and documentation --- R/pareto_smooth.R | 1 - man/pareto_khat.Rd | 1 - man/pareto_smooth.Rd | 22 +++++++++------------- 3 files changed, 9 insertions(+), 15 deletions(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 45b5e48b..67143ce2 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -202,7 +202,6 @@ pareto_smooth.default <- function(x, if (ndraws_tail >= 5) { ord <- sort.int(x, index.return = TRUE) if (abs(max_tail - min_tail) < .Machine$double.eps / 100) { - print(paste(max_tail, min_tail)) warning( "Can't fit generalized Pareto distribution ", "because all tail values are the same.", diff --git a/man/pareto_khat.Rd b/man/pareto_khat.Rd index 842bb3d0..a580aa15 100644 --- a/man/pareto_khat.Rd +++ b/man/pareto_khat.Rd @@ -62,5 +62,4 @@ pareto_khat(mu) d <- as_draws_rvars(example_draws("multi_normal")) pareto_khat(d$Sigma) - } diff --git a/man/pareto_smooth.Rd b/man/pareto_smooth.Rd index 5d4ebd62..0c3a0d17 100644 --- a/man/pareto_smooth.Rd +++ b/man/pareto_smooth.Rd @@ -8,7 +8,7 @@ \usage{ pareto_smooth(x, ...) -\method{pareto_smooth}{rvar}(x, ...) +\method{pareto_smooth}{rvar}(x, extra_diags = FALSE, ...) \method{pareto_smooth}{default}( x, @@ -17,7 +17,6 @@ pareto_smooth(x, ...) ndraws_tail = NULL, extra_diags = FALSE, verbose = FALSE, - return_smoothed = TRUE, ... ) } @@ -31,6 +30,11 @@ pareto_smooth(x, ...) \item{...}{Arguments passed to individual methods (if applicable).} +\item{extra_diags}{(logical) Should extra Pareto smoothing +diagnostics be included in output? If \code{TRUE}, \code{min_ss}, +\code{khat_threshold} and \code{convergence_rate} for the calculated k +value will be returned. Default is \code{FALSE}.} + \item{tail}{(string) The tail to smooth: \itemize{ \item \code{"right"}: smooth only the right (upper) tail @@ -49,22 +53,14 @@ from MCMC.} ceiling(3 * sqrt(length(x) / r_eff)) if length(x) > 225 and length(x) / 5 otherwise (see Appendix H in Vehtari et al. (2022)).} -\item{extra_diags}{(logical) Should extra Pareto smoothing -diagnostics be included in output? If \code{TRUE}, \code{min_ss}, -\code{khat_threshold} and \code{convergence_rate} for the calculated k -value will be returned. Default is \code{FALSE}.} - \item{verbose}{(logical) Should diagnostic messages be printed? If \code{TRUE}, messages related to Pareto smoothing diagnostics will be printed. Default is \code{FALSE}.} - -\item{return_smoothed}{(logical) Should the smoothed x be -returned? If \code{TRUE}, returns smoothed x. If \code{FALSE}, acts the -same as \code{pareto_khat}. Default is \code{TRUE}.} } \value{ -List containing vector of smoothed values and vector of -pareto smoothing diagnostics. +List containing vector \code{x} of smoothed values and list +\code{diagnostics} of Pareto smoothing diagnostics (\code{khat}, and +optionally \code{min_ss}, \code{khat_threshold}, and \code{convergence_rate}). } \description{ Smooth the tail(s) of x by fitting a generalized Pareto From ef9a2b3703f0d9edb5804225b03f09517d568714 Mon Sep 17 00:00:00 2001 From: n-kall Date: Tue, 9 May 2023 11:54:23 +0300 Subject: [PATCH 17/30] remove currently unused genralized pareto distribution functions --- NAMESPACE | 5 -- R/gpd.R | 120 ++-------------------------------------------- R/pareto_smooth.R | 2 +- man/GPD.Rd | 42 ---------------- 4 files changed, 4 insertions(+), 165 deletions(-) delete mode 100644 man/GPD.Rd diff --git a/NAMESPACE b/NAMESPACE index 78079afc..ad92fad7 100755 --- a/NAMESPACE +++ b/NAMESPACE @@ -400,7 +400,6 @@ export(chain_ids) export(default_convergence_measures) export(default_mcse_measures) export(default_summary_measures) -export(dgpd) export(diag) export(dissent) export(draw_ids) @@ -423,7 +422,6 @@ export(example_draws) export(extract_variable) export(extract_variable_matrix) export(for_each_draw) -export(gpdfit) export(is_draws) export(is_draws_array) export(is_draws_df) @@ -450,8 +448,6 @@ export(nvariables) export(order_draws) export(pareto_khat) export(pareto_smooth) -export(pgpd) -export(qgpd) export(quantile2) export(r_scale) export(rdo) @@ -460,7 +456,6 @@ export(repair_draws) export(resample_draws) export(reserved_variables) export(rfun) -export(rgpd) export(rhat) export(rhat_basic) export(rstar) diff --git a/R/gpd.R b/R/gpd.R index a4f0f927..6a751dd0 100644 --- a/R/gpd.R +++ b/R/gpd.R @@ -1,71 +1,5 @@ -#' The Generalized Pareto Distribution -#' -#' Density, distribution function, quantile function and random generation -#' for the generalized Pareto distribution with location equal to \code{mu}, -#' scale equal to \code{sigma}, and shape equal to \code{k}. -#' -#' @param x,q vector of quantiles. -#' @param p vector of probabilities. -#' @param n number of observations. If \code{length(n) > 1}, the length is taken -#' to be the number required. -#' @param mu scalar location parameter -#' @param sigma scalar, positive scale parameter -#' @param k scalar shape parameter -#' @param log,log.p logical; if TRUE, probabilities p are given as log(p). -#' @param lower.tail logical; if TRUE (default), probabilities are \eqn{P[X \le x]} -#' otherwise, \eqn{P[X > x]}. -#' -#' @name GPD - - -#' @rdname GPD -#' @export -dgpd <- function(x, mu = 0, sigma = 1, k = 0, log = FALSE) { - stopifnot(length(mu) == 1 && length(sigma) == 1 && length(k) == 1) - if (is.na(sigma) || sigma <= 0) { - return(rep(NaN, length(x))) - } - d <- (x - mu) / sigma - ind <- (d > 0) & ((1 + k * d) > 0) - ind[is.na(ind)] <- FALSE - if (k == 0) { - d[ind] <- -d[ind] - log(sigma) - } else { - d[ind] <- log1p(k * d[ind]) * -(1 / k + 1) - log(sigma) - } - d[!ind] <- -Inf - if (!log) { - d <- exp(d) - } - d -} - - -#' @rdname GPD -#' @export -pgpd <- function(q, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, log.p = FALSE) { - stopifnot(length(mu) == 1 && length(sigma) == 1 && length(k) == 1) - if (is.na(sigma) || sigma <= 0) { - return(rep(NaN, length(q))) - } - q <- pmax(q - mu, 0) / sigma - if (k == 0) { - p <- 1 - exp(-q) - } else { - p <- -expm1(log(pmax(1 + k * q, 0)) * -(1 / k)) - } - if (!lower.tail) { - p <- 1 - p - } - if (log.p) { - p <- log(p) - } - p -} - -#' @rdname GPD -#' @export -qgpd <- function(p, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, log.p = FALSE) { +# quantile function for generalized pareto +qgeneralized_pareto <- function(p, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, log.p = FALSE) { stopifnot(length(mu) == 1 && length(sigma) == 1 && length(k) == 1) if (is.na(sigma) || sigma <= 0) { return(rep(NaN, length(p))) @@ -84,24 +18,6 @@ qgpd <- function(p, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, log.p = FALSE) q } -#' @rdname GPD -#' @export -rgpd <- function(n, mu = 0, sigma = 1, k = 0) { - stopifnot( - length(n) == 1 && length(mu) == 1 && length(sigma) == 1 && length(k) == 1 - ) - if (is.na(sigma) || sigma <= 0) { - return(rep(NaN, n)) - } - if (k == 0) { - r <- mu + sigma * rexp(n) - } else { - r <- mu + sigma * expm1(-k * log(runif(n))) / k - } - r -} - - #' Estimate parameters of the Generalized Pareto distribution #' #' Given a sample \eqn{x}, Estimate the parameters \eqn{k} and \eqn{\sigma} of @@ -111,7 +27,6 @@ rgpd <- function(n, mu = 0, sigma = 1, k = 0) { #' case of MCMC samples). The weakly informative prior is a Gaussian prior #' centered at 0.5. #' -#' @export #' @param x A numeric vector. The sample from which to estimate the parameters. #' @param wip Logical indicating whether to adjust \eqn{k} based on a weakly #' informative Gaussian prior centered on 0.5. Defaults to `TRUE`. @@ -160,34 +75,5 @@ gpdfit <- function(x, wip = TRUE, min_grid_pts = 30, sort_x = TRUE) { sigma_hat <- NaN } - list("k" = k_hat, "sigma" = sigma_hat) -} - -# Below is a separate function for the marginal posterior. -# We can extend to have computationally more intensive but more accurate methods. -gpdpost <- function(x, min_grid_pts = 30, sort_x = TRUE) { - - if (sort_x) { - x <- sort.int(x) - } - N <- length(x) - prior <- 3 - M <- min_grid_pts + floor(sqrt(N)) - jj <- seq_len(M) - xstar <- x[floor(N / 4 + 0.5)] # first quartile of sample - theta <- 1 / x[N] + (1 - sqrt(M / (jj - 0.5))) / prior / xstar - k <- matrixStats::rowMeans2(log1p(-theta %o% x)) - l_theta <- N * (log(-theta / k) - k - 1) # profile log-lik - w_theta <- exp(l_theta - matrixStats::logSumExp(l_theta)) # normalize - - # quadrature weights for k are same as for theta - k_w <- w_theta - # quadrature weights are just the normalized likelihoods - # we get the unnormalized posterior by multiplying these by the prior - k_d <- k_w * dgpd(-theta, mu = -1 / x[N], sigma = 1 / prior / xstar, k = 0.5) - # normalize using the trapezoidal rule - Z <- sum((k_d[-M] + k_d[-1]) * (k[-M] - k[-1])) / 2 - k_d <- k_d / Z - - list("k" = k, "k_w" = k_w, "k_d" = k_d) + list(k = k_hat, sigma = sigma_hat) } diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 67143ce2..86928596 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -216,7 +216,7 @@ pareto_smooth.default <- function(x, sigma <- fit$sigma if (is.finite(k)) { p <- (seq_len(ndraws_tail) - 0.5) / ndraws_tail - smoothed <- qgpd(p = p, mu = cutoff, k = k, sigma = sigma) + smoothed <- qgeneralized_pareto(p = p, mu = cutoff, k = k, sigma = sigma) } } } else { diff --git a/man/GPD.Rd b/man/GPD.Rd deleted file mode 100644 index 9f6e5542..00000000 --- a/man/GPD.Rd +++ /dev/null @@ -1,42 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/gpd.R -\name{GPD} -\alias{GPD} -\alias{dgpd} -\alias{pgpd} -\alias{qgpd} -\alias{rgpd} -\title{The Generalized Pareto Distribution} -\usage{ -dgpd(x, mu = 0, sigma = 1, k = 0, log = FALSE) - -pgpd(q, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, log.p = FALSE) - -qgpd(p, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, log.p = FALSE) - -rgpd(n, mu = 0, sigma = 1, k = 0) -} -\arguments{ -\item{x, q}{vector of quantiles.} - -\item{mu}{scalar location parameter} - -\item{sigma}{scalar, positive scale parameter} - -\item{k}{scalar shape parameter} - -\item{log, log.p}{logical; if TRUE, probabilities p are given as log(p).} - -\item{lower.tail}{logical; if TRUE (default), probabilities are \eqn{P[X \le x]} -otherwise, \eqn{P[X > x]}.} - -\item{p}{vector of probabilities.} - -\item{n}{number of observations. If \code{length(n) > 1}, the length is taken -to be the number required.} -} -\description{ -Density, distribution function, quantile function and random generation -for the generalized Pareto distribution with location equal to \code{mu}, -scale equal to \code{sigma}, and shape equal to \code{k}. -} From 077c975f61839cbe6941fe546b8ee239e8172586 Mon Sep 17 00:00:00 2001 From: n-kall Date: Tue, 9 May 2023 12:14:25 +0300 Subject: [PATCH 18/30] handle extra diagnostics in pareto_khat.rvar --- R/pareto_smooth.R | 33 ++++++++++++++++++++++----------- 1 file changed, 22 insertions(+), 11 deletions(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 86928596..39906680 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -26,13 +26,26 @@ pareto_khat.default <- function(x, #' @rdname pareto_khat #' @export -pareto_khat.rvar <- function(x, ...) { - draws_diags <- summarise_rvar_by_element_with_chains(x, pareto_smooth.default, ...) +pareto_khat.rvar <- function(x, extra_diags = FALSE, ...) { + draws_diags <- summarise_rvar_by_element_with_chains(x, pareto_smooth.default, extra_diags = extra_diags, ...) dim(draws_diags) <- dim(draws_diags) %||% length(draws_diags) margins <- seq_along(dim(draws_diags)) - list( - khat = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat) - ) + + if (extra_diags) { + + diags <- list( + khat = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat), + min_ss = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$min_ss), + khat_threshold = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat_threshold), + convergence_rage = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$convergence_rate) + ) + } else { + diags <- list( + khat = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat) + ) + } + + diags } #' Pareto smoothing @@ -202,10 +215,9 @@ pareto_smooth.default <- function(x, if (ndraws_tail >= 5) { ord <- sort.int(x, index.return = TRUE) if (abs(max_tail - min_tail) < .Machine$double.eps / 100) { - warning( + warning_no_call( "Can't fit generalized Pareto distribution ", - "because all tail values are the same.", - call. = FALSE + "because all tail values are the same." ) smoothed <- NULL k <- NA @@ -220,10 +232,9 @@ pareto_smooth.default <- function(x, } } } else { - warning( + warning_no_call( "Can't fit generalized Pareto distribution ", - "because ndraws_tail is less than 5.", - call. = FALSE + "because ndraws_tail is less than 5." ) smoothed <- NULL k <- NA From ca58e6dccfb225eee454da852862fcd4bc9377c4 Mon Sep 17 00:00:00 2001 From: n-kall Date: Tue, 9 May 2023 12:37:24 +0300 Subject: [PATCH 19/30] improve documentation and comments for pareto smoothing functions --- R/gpd.R | 4 +++- R/pareto_smooth.R | 39 ++++++++++++++++++++++++++---------- man/gpdfit.Rd | 41 -------------------------------------- man/pareto_k_diagmsg.Rd | 19 ------------------ man/pareto_khat.Rd | 20 ++++++++++++------- man/pareto_smooth.Rd | 9 +++++++-- man/ps_convergence_rate.Rd | 22 -------------------- man/ps_khat_threshold.Rd | 21 ------------------- man/ps_min_ss.Rd | 20 ------------------- 9 files changed, 52 insertions(+), 143 deletions(-) delete mode 100644 man/gpdfit.Rd delete mode 100644 man/pareto_k_diagmsg.Rd delete mode 100644 man/ps_convergence_rate.Rd delete mode 100644 man/ps_khat_threshold.Rd delete mode 100644 man/ps_min_ss.Rd diff --git a/R/gpd.R b/R/gpd.R index 6a751dd0..a893dd38 100644 --- a/R/gpd.R +++ b/R/gpd.R @@ -1,4 +1,5 @@ -# quantile function for generalized pareto +#' Quantile function for generalized pareto +#' @noRd qgeneralized_pareto <- function(p, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, log.p = FALSE) { stopifnot(length(mu) == 1 && length(sigma) == 1 && length(k) == 1) if (is.na(sigma) || sigma <= 0) { @@ -46,6 +47,7 @@ qgeneralized_pareto <- function(p, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, #' Zhang, J., and Stephens, M. A. (2009). A new and efficient estimation method #' for the generalized Pareto distribution. *Technometrics* **51**, 316-325. #' +#' @noRd gpdfit <- function(x, wip = TRUE, min_grid_pts = 30, sort_x = TRUE) { # see section 4 of Zhang and Stephens (2009) if (sort_x) { diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 39906680..7cd068b9 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -5,7 +5,12 @@ #' #' @template args-pareto #' @template args-methods-dots -#' @return numeric vector of Pareto smoothing diagnostics +#' @return #' List of Pareto smoothing diagnostics: +#' * `khat`: estimated Pareto k shape parameter, and optionally +#' * `min_ss`: minimum sample size for reliable Pareto smoothed estimate +#' * `khat_threshold`: khat-threshold for reliable Pareto smoothed estimate +#' * `convergence_rate`: Pareto smoothed estimate RMSE convergence rate. +#' #' #' @examples #' mu <- extract_variable_matrix(example_draws(), "mu") @@ -56,8 +61,11 @@ pareto_khat.rvar <- function(x, extra_diags = FALSE, ...) { #' @template args-pareto #' @template args-methods-dots #' @return List containing vector `x` of smoothed values and list -#' `diagnostics` of Pareto smoothing diagnostics (`khat`, and -#' optionally `min_ss`, `khat_threshold`, and `convergence_rate`). +#' `diagnostics` of Pareto smoothing diagnostics: +#' * `khat`: estimated Pareto k shape parameter, and optionally +#' * `min_ss`: minimum sample size for reliable Pareto smoothed estimate +#' * `khat_threshold`: khat-threshold for reliable Pareto smoothed estimate +#' * `convergence_rate`: Pareto smoothed estimate RMSE convergence rate. #' #' @examples #' mu <- extract_variable_matrix(example_draws(), "mu") @@ -103,8 +111,8 @@ pareto_smooth.default <- function(x, ndraws_tail = NULL, extra_diags = FALSE, verbose = FALSE, - ...) { - + ...) { + # check for infinite or na values if (should_return_NA(x)) { return(list(x = x, diagnostics = NA_real_)) @@ -188,7 +196,9 @@ pareto_smooth.default <- function(x, return(out) } -# internal function to pareto smooth the tail of a vector +#' Pareto smooth tail +#' internal function to pareto smooth the tail of a vector +#' @noRd .pareto_smooth_tail <- function(x, ndraws_tail, tail = c("right", "left"), @@ -255,8 +265,11 @@ pareto_smooth.default <- function(x, return(out) } -# internal function to calculate the extra diagnostics or a given -# pareto k and sample size S +#' Extra pareto-k diagnostics +#' +#' internal function to calculate the extra diagnostics or a given +#' pareto k and sample size S +#' @noRd .pareto_smooth_extra_diags <- function(k, S, ...) { min_ss <- ps_min_ss(k) @@ -276,11 +289,12 @@ pareto_smooth.default <- function(x, #' #' Given Pareto-k computes the minimum sample size for reliable Pareto #' smoothed estimate (to have small probability of large error) +#' Equqtion (11) in PSIS paper #' @param k pareto k value #' @param ... unused #' @return minimum sample size +#' @noRd ps_min_ss <- function(k, ...) { - # Eq (11) in PSIS paper 10^(1 / (1 - max(0, k))) } @@ -292,6 +306,7 @@ ps_min_ss <- function(k, ...) { #' @param S sample size #' @param ... unused #' @return threshold +#' @noRd ps_khat_threshold <- function(S, ...) { 1 - 1 / log10(S) } @@ -304,6 +319,7 @@ ps_khat_threshold <- function(S, ...) { #' @param S sample size #' @param ... unused #' @return convergence rate +#' @noRd ps_convergence_rate <- function(k, S, ...) { # allow array of k's rate <- numeric(length(k)) @@ -324,8 +340,10 @@ ps_convergence_rate <- function(k, S, ...) { rate } +#' Calculate the tail length from S and r_eff +#' Appendix H in PSIS paper +#' @noRd ps_tail_length <- function(S, r_eff, ...) { - # Appendix H in PSIS paper ifelse(S > 225, ceiling(3 * sqrt(S / r_eff)), S / 5) } @@ -335,6 +353,7 @@ ps_tail_length <- function(S, r_eff, ...) { #' @param diags (numeric) named vector of diagnostic values #' @param ... unused #' @return diagnostic message +#' @noRd pareto_k_diagmsg <- function(diags, ...) { khat <- diags$khat min_ss <- diags$min_ss diff --git a/man/gpdfit.Rd b/man/gpdfit.Rd deleted file mode 100644 index a496f4ad..00000000 --- a/man/gpdfit.Rd +++ /dev/null @@ -1,41 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/gpd.R -\name{gpdfit} -\alias{gpdfit} -\title{Estimate parameters of the Generalized Pareto distribution} -\usage{ -gpdfit(x, wip = TRUE, min_grid_pts = 30, sort_x = TRUE) -} -\arguments{ -\item{x}{A numeric vector. The sample from which to estimate the parameters.} - -\item{wip}{Logical indicating whether to adjust \eqn{k} based on a weakly -informative Gaussian prior centered on 0.5. Defaults to \code{TRUE}.} - -\item{min_grid_pts}{The minimum number of grid points used in the fitting -algorithm. The actual number used is \code{min_grid_pts + floor(sqrt(length(x)))}.} - -\item{sort_x}{If \code{TRUE} (the default), the first step in the fitting -algorithm is to sort the elements of \code{x}. If \code{x} is already -sorted in ascending order then \code{sort_x} can be set to \code{FALSE} to -skip the initial sorting step.} -} -\value{ -A named list with components \code{k} and \code{sigma}. -} -\description{ -Given a sample \eqn{x}, Estimate the parameters \eqn{k} and \eqn{\sigma} of -the generalized Pareto distribution (GPD), assuming the location parameter is -0. By default the fit uses a prior for \eqn{k}, which will stabilize -estimates for very small sample sizes (and low effective sample sizes in the -case of MCMC samples). The weakly informative prior is a Gaussian prior -centered at 0.5. -} -\details{ -Here the parameter \eqn{k} is the negative of \eqn{k} in Zhang & -Stephens (2009). -} -\references{ -Zhang, J., and Stephens, M. A. (2009). A new and efficient estimation method -for the generalized Pareto distribution. \emph{Technometrics} \strong{51}, 316-325. -} diff --git a/man/pareto_k_diagmsg.Rd b/man/pareto_k_diagmsg.Rd deleted file mode 100644 index cf000caa..00000000 --- a/man/pareto_k_diagmsg.Rd +++ /dev/null @@ -1,19 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/pareto_smooth.R -\name{pareto_k_diagmsg} -\alias{pareto_k_diagmsg} -\title{Pareto-k diagnostic message} -\usage{ -pareto_k_diagmsg(diags, ...) -} -\arguments{ -\item{diags}{(numeric) named vector of diagnostic values} - -\item{...}{unused} -} -\value{ -diagnostic message -} -\description{ -Given S and scalar and k, form a diagnostic message string -} diff --git a/man/pareto_khat.Rd b/man/pareto_khat.Rd index a580aa15..3cd9d6cb 100644 --- a/man/pareto_khat.Rd +++ b/man/pareto_khat.Rd @@ -10,7 +10,7 @@ pareto_khat(x, ...) \method{pareto_khat}{default}(x, ...) -\method{pareto_khat}{rvar}(x, ...) +\method{pareto_khat}{rvar}(x, extra_diags = FALSE, ...) } \arguments{ \item{x}{(multiple options) One of: @@ -22,6 +22,11 @@ pareto_khat(x, ...) \item{...}{Arguments passed to individual methods (if applicable).} +\item{extra_diags}{(logical) Should extra Pareto smoothing +diagnostics be included in output? If \code{TRUE}, \code{min_ss}, +\code{khat_threshold} and \code{convergence_rate} for the calculated k +value will be returned. Default is \code{FALSE}.} + \item{tail}{(string) The tail to smooth: \itemize{ \item \code{"right"}: smooth only the right (upper) tail @@ -43,14 +48,15 @@ from MCMC.} \item{verbose}{(logical) Should diagnostic messages be printed? If \code{TRUE}, messages related to Pareto smoothing diagnostics will be printed. Default is \code{FALSE}.} - -\item{extra_diags}{(logical) Should extra Pareto smoothing -diagnostics be included in output? If \code{TRUE}, \code{min_ss}, -\code{khat_threshold} and \code{convergence_rate} for the calculated k -value will be returned. Default is \code{FALSE}.} } \value{ -numeric vector of Pareto smoothing diagnostics +#' List of Pareto smoothing diagnostics: +\itemize{ +\item \code{khat}: estimated Pareto k shape parameter, and optionally +\item \code{min_ss}: minimum sample size for reliable Pareto smoothed estimate +\item \code{khat_threshold}: khat-threshold for reliable Pareto smoothed estimate +\item \code{convergence_rate}: Pareto smoothed estimate RMSE convergence rate. +} } \description{ Estimate Pareto k value from Pareto smoothing the tail(s) of x. For diff --git a/man/pareto_smooth.Rd b/man/pareto_smooth.Rd index 0c3a0d17..fe8327f6 100644 --- a/man/pareto_smooth.Rd +++ b/man/pareto_smooth.Rd @@ -59,8 +59,13 @@ printed. Default is \code{FALSE}.} } \value{ List containing vector \code{x} of smoothed values and list -\code{diagnostics} of Pareto smoothing diagnostics (\code{khat}, and -optionally \code{min_ss}, \code{khat_threshold}, and \code{convergence_rate}). +\code{diagnostics} of Pareto smoothing diagnostics: +\itemize{ +\item \code{khat}: estimated Pareto k shape parameter, and optionally +\item \code{min_ss}: minimum sample size for reliable Pareto smoothed estimate +\item \code{khat_threshold}: khat-threshold for reliable Pareto smoothed estimate +\item \code{convergence_rate}: Pareto smoothed estimate RMSE convergence rate. +} } \description{ Smooth the tail(s) of x by fitting a generalized Pareto diff --git a/man/ps_convergence_rate.Rd b/man/ps_convergence_rate.Rd deleted file mode 100644 index 40948fc4..00000000 --- a/man/ps_convergence_rate.Rd +++ /dev/null @@ -1,22 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/pareto_smooth.R -\name{ps_convergence_rate} -\alias{ps_convergence_rate} -\title{Pareto-smoothing convergence rate} -\usage{ -ps_convergence_rate(k, S, ...) -} -\arguments{ -\item{k}{pareto-k values} - -\item{S}{sample size} - -\item{...}{unused} -} -\value{ -convergence rate -} -\description{ -Given S and scalar or array of k's, compute the relative -convergence rate of PSIS estimate RMSE -} diff --git a/man/ps_khat_threshold.Rd b/man/ps_khat_threshold.Rd deleted file mode 100644 index 1cdfd459..00000000 --- a/man/ps_khat_threshold.Rd +++ /dev/null @@ -1,21 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/pareto_smooth.R -\name{ps_khat_threshold} -\alias{ps_khat_threshold} -\title{Pareto-smoothing k-hat threshold} -\usage{ -ps_khat_threshold(S, ...) -} -\arguments{ -\item{S}{sample size} - -\item{...}{unused} -} -\value{ -threshold -} -\description{ -Given sample size S computes khat threshold for reliable Pareto -smoothed estimate (to have small probability of large error). See -section 3.2.4, equation (13). -} diff --git a/man/ps_min_ss.Rd b/man/ps_min_ss.Rd deleted file mode 100644 index 2b48a38e..00000000 --- a/man/ps_min_ss.Rd +++ /dev/null @@ -1,20 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/pareto_smooth.R -\name{ps_min_ss} -\alias{ps_min_ss} -\title{Pareto-smoothing minimum sample-size} -\usage{ -ps_min_ss(k, ...) -} -\arguments{ -\item{k}{pareto k value} - -\item{...}{unused} -} -\value{ -minimum sample size -} -\description{ -Given Pareto-k computes the minimum sample size for reliable Pareto -smoothed estimate (to have small probability of large error) -} From 7da27e334e480ef3d138594c5b3c33006933a374 Mon Sep 17 00:00:00 2001 From: n-kall Date: Tue, 9 May 2023 12:52:24 +0300 Subject: [PATCH 20/30] add input checks for pareto_smooth --- R/pareto_smooth.R | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 7cd068b9..358a46dd 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -112,9 +112,15 @@ pareto_smooth.default <- function(x, extra_diags = FALSE, verbose = 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(verbose) # 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_)) } From e2330b353319f05348b1b34235206184f2674665 Mon Sep 17 00:00:00 2001 From: n-kall Date: Tue, 9 May 2023 13:26:11 +0300 Subject: [PATCH 21/30] fix documentation issue in pareto_khat --- R/pareto_smooth.R | 19 +++++++++++++++---- man/pareto_khat.Rd | 28 ++++++++++++++++++---------- 2 files changed, 33 insertions(+), 14 deletions(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 358a46dd..b16427a6 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -5,13 +5,12 @@ #' #' @template args-pareto #' @template args-methods-dots -#' @return #' List of Pareto smoothing diagnostics: +#' @return List of Pareto smoothing diagnostics: #' * `khat`: estimated Pareto k shape parameter, and optionally #' * `min_ss`: minimum sample size for reliable Pareto smoothed estimate #' * `khat_threshold`: khat-threshold for reliable Pareto smoothed estimate #' * `convergence_rate`: Pareto smoothed estimate RMSE convergence rate. #' -#' #' @examples #' mu <- extract_variable_matrix(example_draws(), "mu") #' pareto_khat(mu) @@ -24,9 +23,21 @@ pareto_khat <- function(x, ...) UseMethod("pareto_khat") #' @rdname pareto_khat #' @export pareto_khat.default <- function(x, + tail = c("both", "right", "left"), + r_eff = NULL, + ndraws_tail = NULL, + extra_diags = FALSE, + verbose = FALSE, ...) { - out <- pareto_smooth.default(x, ...) - return(out$diagnostics) + smoothed <- pareto_smooth.default( + x, + tail = tail, + r_eff = r_eff, + ndraws_tail = ndraws_tail, + extra_diags = extra_diags, + verbose = verbose, + ...) + return(smoothed$diagnostics) } #' @rdname pareto_khat diff --git a/man/pareto_khat.Rd b/man/pareto_khat.Rd index 3cd9d6cb..9eded27c 100644 --- a/man/pareto_khat.Rd +++ b/man/pareto_khat.Rd @@ -8,7 +8,15 @@ \usage{ pareto_khat(x, ...) -\method{pareto_khat}{default}(x, ...) +\method{pareto_khat}{default}( + x, + tail = c("both", "right", "left"), + r_eff = NULL, + ndraws_tail = NULL, + extra_diags = FALSE, + verbose = FALSE, + ... +) \method{pareto_khat}{rvar}(x, extra_diags = FALSE, ...) } @@ -22,11 +30,6 @@ pareto_khat(x, ...) \item{...}{Arguments passed to individual methods (if applicable).} -\item{extra_diags}{(logical) Should extra Pareto smoothing -diagnostics be included in output? If \code{TRUE}, \code{min_ss}, -\code{khat_threshold} and \code{convergence_rate} for the calculated k -value will be returned. Default is \code{FALSE}.} - \item{tail}{(string) The tail to smooth: \itemize{ \item \code{"right"}: smooth only the right (upper) tail @@ -36,21 +39,26 @@ value will be 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 +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)).} -\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.} +\item{extra_diags}{(logical) Should extra Pareto smoothing +diagnostics be included in output? If \code{TRUE}, \code{min_ss}, +\code{khat_threshold} and \code{convergence_rate} for the calculated k +value will be returned. Default is \code{FALSE}.} \item{verbose}{(logical) Should diagnostic messages be printed? If \code{TRUE}, messages related to Pareto smoothing diagnostics will be printed. Default is \code{FALSE}.} } \value{ -#' List of Pareto smoothing diagnostics: +List of Pareto smoothing diagnostics: \itemize{ \item \code{khat}: estimated Pareto k shape parameter, and optionally \item \code{min_ss}: minimum sample size for reliable Pareto smoothed estimate From 60fe39c4ce40c38cfbed2589ca93efcb63c37ec3 Mon Sep 17 00:00:00 2001 From: n-kall Date: Tue, 9 May 2023 14:44:16 +0300 Subject: [PATCH 22/30] typofix 'convergence_rage' -> 'convergence_rate' --- R/pareto_smooth.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index b16427a6..b9fc13c4 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -53,7 +53,7 @@ pareto_khat.rvar <- function(x, extra_diags = FALSE, ...) { khat = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat), min_ss = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$min_ss), khat_threshold = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat_threshold), - convergence_rage = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$convergence_rate) + convergence_rate = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$convergence_rate) ) } else { diags <- list( @@ -100,7 +100,7 @@ pareto_smooth.rvar <- function(x, extra_diags = FALSE, ...) { khat = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat), min_ss = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$min_ss), khat_threshold = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat_threshold), - convergence_rage = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$convergence_rate) + convergence_rate = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$convergence_rate) ) } else { diags <- list( From 982d408489afc31ddb0045626f1154380179b3ff Mon Sep 17 00:00:00 2001 From: n-kall Date: Fri, 12 May 2023 18:14:23 +0300 Subject: [PATCH 23/30] set minimum ss when pareto-k >=1 to infinity --- R/pareto_smooth.R | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index b9fc13c4..4131d812 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -312,7 +312,11 @@ pareto_smooth.default <- function(x, #' @return minimum sample size #' @noRd ps_min_ss <- function(k, ...) { - 10^(1 / (1 - max(0, k))) + if (k < 1) { + out <- 10^(1 / (1 - max(0, k))) + } else { + out <- Inf + } } #' Pareto-smoothing k-hat threshold From ad68f27732b78a6e01ccf081fc6024288d599611 Mon Sep 17 00:00:00 2001 From: n-kall <33577035+n-kall@users.noreply.github.com> Date: Sat, 13 May 2023 09:50:32 +0300 Subject: [PATCH 24/30] set min_ss based on pareto_k to inf if k >= 1 --- R/pareto_smooth.R | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index b9fc13c4..071c8b34 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -312,9 +312,15 @@ pareto_smooth.default <- function(x, #' @return minimum sample size #' @noRd ps_min_ss <- function(k, ...) { - 10^(1 / (1 - max(0, k))) + if (k < 1) { + out <- 10^(1 / (1 - max(0, k))) + } else { + out <- Inf + } + out } + #' Pareto-smoothing k-hat threshold #' #' Given sample size S computes khat threshold for reliable Pareto From 3e823e7e0299b68c7ccba47c810642fd0f01f3a4 Mon Sep 17 00:00:00 2001 From: n-kall Date: Tue, 16 May 2023 14:41:19 +0300 Subject: [PATCH 25/30] add pareto smoothing reference in docs --- R/pareto_smooth.R | 2 ++ man/pareto_khat.Rd | 5 +++++ man/pareto_smooth.Rd | 5 +++++ 3 files changed, 12 insertions(+) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 071c8b34..79a6187d 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -5,6 +5,7 @@ #' #' @template args-pareto #' @template args-methods-dots +#' @template ref-vehtari-paretosmooth-2022 #' @return List of Pareto smoothing diagnostics: #' * `khat`: estimated Pareto k shape parameter, and optionally #' * `min_ss`: minimum sample size for reliable Pareto smoothed estimate @@ -71,6 +72,7 @@ pareto_khat.rvar <- function(x, extra_diags = FALSE, ...) { #' #' @template args-pareto #' @template args-methods-dots +#' @template ref-vehtari-paretosmooth-2022 #' @return List containing vector `x` of smoothed values and list #' `diagnostics` of Pareto smoothing diagnostics: #' * `khat`: estimated Pareto k shape parameter, and optionally diff --git a/man/pareto_khat.Rd b/man/pareto_khat.Rd index 9eded27c..cc024b9f 100644 --- a/man/pareto_khat.Rd +++ b/man/pareto_khat.Rd @@ -77,3 +77,8 @@ pareto_khat(mu) d <- as_draws_rvars(example_draws("multi_normal")) 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 +} diff --git a/man/pareto_smooth.Rd b/man/pareto_smooth.Rd index fe8327f6..e5344d59 100644 --- a/man/pareto_smooth.Rd +++ b/man/pareto_smooth.Rd @@ -78,3 +78,8 @@ pareto_smooth(mu) d <- as_draws_rvars(example_draws("multi_normal")) 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 +} From 4c9aa75d5a5667ce127e52081ae2290129ee0bfb Mon Sep 17 00:00:00 2001 From: n-kall Date: Wed, 17 May 2023 12:24:56 +0300 Subject: [PATCH 26/30] add new function pareto_diags and adjust existing functions --- NAMESPACE | 3 + R/pareto_smooth.R | 222 ++++++++++++++++++++-------- man-roxygen/args-pareto.R | 14 +- man/pareto_khat.Rd | 28 +--- man/pareto_smooth.Rd | 47 +++--- tests/testthat/test-pareto_smooth.R | 18 ++- 6 files changed, 219 insertions(+), 113 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 6e6d7e2b..c34bf942 100755 --- a/NAMESPACE +++ b/NAMESPACE @@ -223,6 +223,8 @@ S3method(order_draws,draws_list) S3method(order_draws,draws_matrix) S3method(order_draws,draws_rvars) S3method(order_draws,rvar) +S3method(pareto_diags,default) +S3method(pareto_diags,rvar) S3method(pareto_khat,default) S3method(pareto_khat,rvar) S3method(pareto_smooth,default) @@ -448,6 +450,7 @@ export(ndraws) export(niterations) export(nvariables) export(order_draws) +export(pareto_diags) export(pareto_khat) export(pareto_smooth) export(quantile2) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 79a6187d..3418cf88 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -1,17 +1,12 @@ #' Pareto khat diagnostic #' -#' Estimate Pareto k value from Pareto smoothing the tail(s) of x. For +#' Estimate Pareto k value by fitting a Generalized Pareto Distribution to one or two tails of x. For #' further details see Vehtari et al. (2022). #' #' @template args-pareto #' @template args-methods-dots #' @template ref-vehtari-paretosmooth-2022 -#' @return List of Pareto smoothing diagnostics: -#' * `khat`: estimated Pareto k shape parameter, and optionally -#' * `min_ss`: minimum sample size for reliable Pareto smoothed estimate -#' * `khat_threshold`: khat-threshold for reliable Pareto smoothed estimate -#' * `convergence_rate`: Pareto smoothed estimate RMSE convergence rate. -#' +#' @return `khat` estimated Generalized Pareto Distribution shape parameter k #' @examples #' mu <- extract_variable_matrix(example_draws(), "mu") #' pareto_khat(mu) @@ -27,7 +22,6 @@ pareto_khat.default <- function(x, tail = c("both", "right", "left"), r_eff = NULL, ndraws_tail = NULL, - extra_diags = FALSE, verbose = FALSE, ...) { smoothed <- pareto_smooth.default( @@ -35,50 +29,139 @@ pareto_khat.default <- function(x, tail = tail, r_eff = r_eff, ndraws_tail = ndraws_tail, - extra_diags = extra_diags, verbose = verbose, + return_k = TRUE, + smooth_draws = FALSE, ...) return(smoothed$diagnostics) } #' @rdname pareto_khat #' @export -pareto_khat.rvar <- function(x, extra_diags = FALSE, ...) { - draws_diags <- summarise_rvar_by_element_with_chains(x, pareto_smooth.default, extra_diags = extra_diags, ...) +pareto_khat.rvar <- function(x, ...) { + draws_diags <- summarise_rvar_by_element_with_chains( + x, + pareto_smooth.default, + return_k = TRUE, + smooth_draws = FALSE, + ... + ) dim(draws_diags) <- dim(draws_diags) %||% length(draws_diags) margins <- seq_along(dim(draws_diags)) - if (extra_diags) { + diags <- list( + khat = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat) + ) - diags <- list( - khat = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat), - min_ss = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$min_ss), - khat_threshold = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat_threshold), - convergence_rate = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$convergence_rate) - ) - } else { - diags <- list( - khat = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat) - ) - } + diags +} + + +#' Pareto smoothing diagnostics +#' +#' Compute diagnostics for Pareto smoothing the tail draws of x by +#' replacing tail draws by order statistics of a generalized Pareto +#' distribution fit to the tail(s). For further details see Vehtari et +#' al. (2022). +#' +#' @template args-pareto +#' @template args-methods-dots +#' @template ref-vehtari-paretosmooth-2022 +#' @return List of Pareto smoothing diagnostics: +#' * `khat`: estimated Pareto k shape parameter, +#' * `min_ss`: minimum sample size for reliable Pareto smoothed estimate, +#' * `khat_threshold`: khat-threshold for reliable Pareto smoothed estimate, +#' * `convergence_rate`: Pareto smoothed estimate RMSE convergence rate. +#' +#' @examples +#' mu <- extract_variable_matrix(example_draws(), "mu") +#' pareto_diags(mu) +#' +#' d <- as_draws_rvars(example_draws("multi_normal")) +#' pareto_diags(d$Sigma) +#' @export +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, + ...) { - diags + smoothed <- pareto_smooth.default( + x, + tail = tail, + r_eff = r_eff, + ndraws_tail = ndraws_tail, + return_k = TRUE, + extra_diags = TRUE, + verbose = verbose, + smooth_draws = FALSE, + ...) + + return(smoothed$diagnostics) + +} + + +#' @rdname pareto_diags +#' @export +pareto_diags.rvar <- function(x, ...) { + draws_diags <- summarise_rvar_by_element_with_chains( + x, + pareto_smooth.default, + return_k = TRUE, + smooth_draws = FALSE, + extra_diags = TRUE, + ... + ) + + dim(draws_diags) <- dim(draws_diags) %||% length(draws_diags) + margins <- seq_along(dim(draws_diags)) + + diags <- list( + khat = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat), + min_ss = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$min_ss), + khat_threshold = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat_threshold), + convergence_rate = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$convergence_rate) + ) + + diags } + + #' Pareto smoothing #' -#' Smooth the tail(s) of x by fitting a generalized Pareto -#' distribution. For further details see Vehtari et al. (2022). +#' Smooth the tail draws of x by replacing tail draws by order +#' statistics of a generalized Pareto distribution fit to the +#' tail(s). For further details see Vehtari et al. (2022). #' #' @template args-pareto +#' @param return_k (logical) Should the Pareto khat be included in +#' output? If `TRUE`, output will be a list containing of smoothed +#' draws and diagnostics. Default is `TRUE`. +#' @param extra_diags (logical) Should extra Pareto khat diagnostics +#' be included in output? If `TRUE`, `min_ss`, `khat_threshold` and +#' `convergence_rate` for the calculated k value will be +#' returned. Default is `FALSE`. #' @template args-methods-dots #' @template ref-vehtari-paretosmooth-2022 -#' @return List containing vector `x` of smoothed values and list -#' `diagnostics` of Pareto smoothing diagnostics: -#' * `khat`: estimated Pareto k shape parameter, and optionally -#' * `min_ss`: minimum sample size for reliable Pareto smoothed estimate -#' * `khat_threshold`: khat-threshold for reliable Pareto smoothed estimate -#' * `convergence_rate`: Pareto smoothed estimate RMSE convergence rate. +#' @return Either a vector `x` of smoothed values and or a list +#' containing the vector `x` and `diagnostics` of Pareto smoothing +#' diagnostics: +#' * `khat`: estimated Pareto k shape parameter, and +#' optionally +#' * `min_ss`: minimum sample size for reliable Pareto +#' smoothed estimate +#' * `khat_threshold`: khat-threshold for reliable +#' Pareto smoothed estimat +#' * `convergence_rate`: Pareto smoothed +#' estimate RMSE convergence rate. #' #' @examples #' mu <- extract_variable_matrix(example_draws(), "mu") @@ -91,29 +174,38 @@ pareto_smooth <- function(x, ...) UseMethod("pareto_smooth") #' @rdname pareto_smooth #' @export -pareto_smooth.rvar <- function(x, extra_diags = FALSE, ...) { - draws_diags <- summarise_rvar_by_element_with_chains(x, pareto_smooth.default, extra_diags = extra_diags, ...) - dim(draws_diags) <- dim(draws_diags) %||% length(draws_diags) - margins <- seq_along(dim(draws_diags)) +pareto_smooth.rvar <- function(x, return_k = TRUE, extra_diags = FALSE, ...) { if (extra_diags) { - - diags <- list( - khat = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat), - min_ss = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$min_ss), - khat_threshold = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat_threshold), - convergence_rate = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$convergence_rate) - ) - } else { - diags <- list( - khat = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat) - ) + return_k <- TRUE } + + draws_diags <- summarise_rvar_by_element_with_chains(x, pareto_smooth.default, return_k = return_k, extra_diags = extra_diags, ...) + dim(draws_diags) <- dim(draws_diags) %||% length(draws_diags) + margins <- seq_along(dim(draws_diags)) - list( + if (return_k) { + if (extra_diags) { + + diags <- list( + khat = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat), + min_ss = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$min_ss), + khat_threshold = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat_threshold), + convergence_rate = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$convergence_rate) + ) + } else { + diags <- list( + khat = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat) + ) + } + out <- list( 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)) + } + out } #' @rdname pareto_smooth @@ -122,6 +214,7 @@ pareto_smooth.default <- function(x, tail = c("both", "right", "left"), r_eff = NULL, ndraws_tail = NULL, + return_k = TRUE, extra_diags = FALSE, verbose = FALSE, ...) { @@ -129,6 +222,7 @@ pareto_smooth.default <- function(x, 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) # check for infinite or na values @@ -169,7 +263,8 @@ pareto_smooth.default <- function(x, smoothed <- .pareto_smooth_tail( x, ndraws_tail = ndraws_tail, - tail = "left" + tail = "left", + ... ) left_k <- smoothed$k @@ -177,7 +272,8 @@ pareto_smooth.default <- function(x, smoothed <-.pareto_smooth_tail( x = smoothed$x, ndraws_tail = ndraws_tail, - tail = "right" + tail = "right", + ... ) right_k <- smoothed$k @@ -188,29 +284,34 @@ pareto_smooth.default <- function(x, smoothed <- .pareto_smooth_tail( x, ndraws_tail = ndraws_tail, - tail = tail + tail = tail, + ... ) k <- smoothed$k x <- smoothed$x } - diags <- list(khat = k) + diags_list <- list(khat = k) if (extra_diags) { ext_diags <- .pareto_smooth_extra_diags(k, S) - diags <- c(diags, ext_diags) + diags_list <- c(diags_list, ext_diags) } if (verbose) { if (!extra_diags) { - diags <- .pareto_smooth_extra_diags(diags$khat, length(x)) + diags_list <- .pareto_smooth_extra_diags(diags_list$khat, length(x)) } pareto_k_diagmsg( - diags = diags + diags = diags_list ) } - out <- list(x = x, diagnostics = diags) + if (return_k) { + out <- list(x = x, diagnostics = diags_list) + } else { + out <- x + } return(out) } @@ -220,6 +321,7 @@ pareto_smooth.default <- function(x, #' @noRd .pareto_smooth_tail <- function(x, ndraws_tail, + smooth_draws = TRUE, tail = c("right", "left"), ... ) { @@ -255,9 +357,11 @@ pareto_smooth.default <- function(x, fit <- gpdfit(draws_tail - cutoff, sort_x = FALSE) k <- fit$k sigma <- fit$sigma - if (is.finite(k)) { + 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) + } else { + smoothed <- NULL } } } else { @@ -308,7 +412,7 @@ pareto_smooth.default <- function(x, #' #' Given Pareto-k computes the minimum sample size for reliable Pareto #' smoothed estimate (to have small probability of large error) -#' Equqtion (11) in PSIS paper +#' Equation (11) in PSIS paper #' @param k pareto k value #' @param ... unused #' @return minimum sample size diff --git a/man-roxygen/args-pareto.R b/man-roxygen/args-pareto.R index daff8977..51cae9c3 100644 --- a/man-roxygen/args-pareto.R +++ b/man-roxygen/args-pareto.R @@ -2,10 +2,10 @@ #' - A matrix of draws for a single variable (iterations x chains). See #' [extract_variable_matrix()]. #' - An [`rvar`]. -#' @param tail (string) The tail to smooth: -#' * `"right"`: smooth only the right (upper) tail -#' * `"left"`: smooth only the left (lower) tail -#' * `"both"`: smooth both tails and return the maximum k value +#' @param tail (string) The tail to diagnose/smooth: +#' * `"right"`: diagnose/smooth only the right (upper) tail +#' * `"left"`: diagnose/smooth only the left (lower) tail +#' * `"both"`: diagnose/smooth both tails and return the maximum k value #' #' The default is `"both"`. #' @param ndraws_tail (numeric) number of draws for the tail. If @@ -16,9 +16,5 @@ #' `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 smoothing diagnostics will be +#' `TRUE`, messages related to Pareto diagnostics will be #' printed. Default is `FALSE`. -#' @param extra_diags (logical) Should extra Pareto smoothing -#' diagnostics be included in output? If `TRUE`, `min_ss`, -#' `khat_threshold` and `convergence_rate` for the calculated k -#' value will be returned. Default is `FALSE`. diff --git a/man/pareto_khat.Rd b/man/pareto_khat.Rd index cc024b9f..2c3a148b 100644 --- a/man/pareto_khat.Rd +++ b/man/pareto_khat.Rd @@ -13,12 +13,11 @@ pareto_khat(x, ...) tail = c("both", "right", "left"), r_eff = NULL, ndraws_tail = NULL, - extra_diags = FALSE, verbose = FALSE, ... ) -\method{pareto_khat}{rvar}(x, extra_diags = FALSE, ...) +\method{pareto_khat}{rvar}(x, ...) } \arguments{ \item{x}{(multiple options) One of: @@ -30,11 +29,11 @@ pareto_khat(x, ...) \item{...}{Arguments passed to individual methods (if applicable).} -\item{tail}{(string) The tail to smooth: +\item{tail}{(string) The tail to diagnose/smooth: \itemize{ -\item \code{"right"}: smooth only the right (upper) tail -\item \code{"left"}: smooth only the left (lower) tail -\item \code{"both"}: smooth both tails and return the maximum k value +\item \code{"right"}: diagnose/smooth only the right (upper) tail +\item \code{"left"}: diagnose/smooth only the left (lower) tail +\item \code{"both"}: diagnose/smooth both tails and return the maximum k value } The default is \code{"both"}.} @@ -48,26 +47,15 @@ from MCMC.} ceiling(3 * sqrt(length(x) / r_eff)) if length(x) > 225 and length(x) / 5 otherwise (see Appendix H in Vehtari et al. (2022)).} -\item{extra_diags}{(logical) Should extra Pareto smoothing -diagnostics be included in output? If \code{TRUE}, \code{min_ss}, -\code{khat_threshold} and \code{convergence_rate} for the calculated k -value will be returned. Default is \code{FALSE}.} - \item{verbose}{(logical) Should diagnostic messages be printed? If -\code{TRUE}, messages related to Pareto smoothing diagnostics will be +\code{TRUE}, messages related to Pareto diagnostics will be printed. Default is \code{FALSE}.} } \value{ -List of Pareto smoothing diagnostics: -\itemize{ -\item \code{khat}: estimated Pareto k shape parameter, and optionally -\item \code{min_ss}: minimum sample size for reliable Pareto smoothed estimate -\item \code{khat_threshold}: khat-threshold for reliable Pareto smoothed estimate -\item \code{convergence_rate}: Pareto smoothed estimate RMSE convergence rate. -} +\code{khat} estimated Generalized Pareto Distribution shape parameter k } \description{ -Estimate Pareto k value from Pareto smoothing the tail(s) of x. For +Estimate Pareto k value by fitting a Generalized Pareto Distribution to one or two tails of x. For further details see Vehtari et al. (2022). } \examples{ diff --git a/man/pareto_smooth.Rd b/man/pareto_smooth.Rd index e5344d59..52d085d0 100644 --- a/man/pareto_smooth.Rd +++ b/man/pareto_smooth.Rd @@ -8,13 +8,14 @@ \usage{ pareto_smooth(x, ...) -\method{pareto_smooth}{rvar}(x, extra_diags = FALSE, ...) +\method{pareto_smooth}{rvar}(x, return_k = TRUE, extra_diags = FALSE, ...) \method{pareto_smooth}{default}( x, tail = c("both", "right", "left"), r_eff = NULL, ndraws_tail = NULL, + return_k = TRUE, extra_diags = FALSE, verbose = FALSE, ... @@ -30,16 +31,20 @@ pareto_smooth(x, ...) \item{...}{Arguments passed to individual methods (if applicable).} -\item{extra_diags}{(logical) Should extra Pareto smoothing -diagnostics be included in output? If \code{TRUE}, \code{min_ss}, -\code{khat_threshold} and \code{convergence_rate} for the calculated k -value will be returned. Default is \code{FALSE}.} +\item{return_k}{(logical) Should the Pareto khat be included in +output? If \code{TRUE}, output will be a list containing of smoothed +draws and diagnostics. Default is \code{TRUE}.} -\item{tail}{(string) The tail to smooth: +\item{extra_diags}{(logical) Should extra Pareto khat diagnostics +be included in output? If \code{TRUE}, \code{min_ss}, \code{khat_threshold} and +\code{convergence_rate} for the calculated k value will be +returned. Default is \code{FALSE}.} + +\item{tail}{(string) The tail to diagnose/smooth: \itemize{ -\item \code{"right"}: smooth only the right (upper) tail -\item \code{"left"}: smooth only the left (lower) tail -\item \code{"both"}: smooth both tails and return the maximum k value +\item \code{"right"}: diagnose/smooth only the right (upper) tail +\item \code{"left"}: diagnose/smooth only the left (lower) tail +\item \code{"both"}: diagnose/smooth both tails and return the maximum k value } The default is \code{"both"}.} @@ -54,22 +59,28 @@ ceiling(3 * sqrt(length(x) / r_eff)) if length(x) > 225 and 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 smoothing diagnostics will be +\code{TRUE}, messages related to Pareto diagnostics will be printed. Default is \code{FALSE}.} } \value{ -List containing vector \code{x} of smoothed values and list -\code{diagnostics} of Pareto smoothing diagnostics: +Either a vector \code{x} of smoothed values and or a list +containing the vector \code{x} and \code{diagnostics} of Pareto smoothing +diagnostics: \itemize{ -\item \code{khat}: estimated Pareto k shape parameter, and optionally -\item \code{min_ss}: minimum sample size for reliable Pareto smoothed estimate -\item \code{khat_threshold}: khat-threshold for reliable Pareto smoothed estimate -\item \code{convergence_rate}: Pareto smoothed estimate RMSE convergence rate. +\item \code{khat}: estimated Pareto k shape parameter, and +optionally +\item \code{min_ss}: minimum sample size for reliable Pareto +smoothed estimate +\item \code{khat_threshold}: khat-threshold for reliable +Pareto smoothed estimat +\item \code{convergence_rate}: Pareto smoothed +estimate RMSE convergence rate. } } \description{ -Smooth the tail(s) of x by fitting a generalized Pareto -distribution. For further details see Vehtari et al. (2022). +Smooth the tail draws of x by replacing tail draws by order +statistics of a generalized Pareto distribution fit to the +tail(s). For further details see Vehtari et al. (2022). } \examples{ mu <- extract_variable_matrix(example_draws(), "mu") diff --git a/tests/testthat/test-pareto_smooth.R b/tests/testthat/test-pareto_smooth.R index e2bddb60..28b16a4a 100644 --- a/tests/testthat/test-pareto_smooth.R +++ b/tests/testthat/test-pareto_smooth.R @@ -80,9 +80,10 @@ test_that("pareto_khat diagnostics messages are as expected", { }) test_that("pareto_khat extra diags are returned", { + tau <- extract_variable_matrix(example_draws(), "tau") - pk <- pareto_khat(tau, extra_diag = TRUE) + pk <- pareto_diags(tau) expect_true(all(names(pk) == c("khat", "min_ss", "khat_threshold", "convergence_rate"))) }) @@ -90,9 +91,9 @@ test_that("pareto_khat extra diags are returned", { test_that("pareto_khat diagnostics handles verbose argument", { tau <- extract_variable_matrix(example_draws(), "tau") - expect_message(pareto_khat(tau, extra_diag = TRUE, verbose = TRUE)) - expect_message(pareto_khat(tau, extra_diag = FALSE, verbose = TRUE)) - expect_no_message(pareto_khat(tau, extra_diag = TRUE, verbose = FALSE)) + expect_message(pareto_khat(tau, extra_diags = TRUE, verbose = TRUE)) + expect_message(pareto_khat(tau, extra_diags = FALSE, verbose = TRUE)) + expect_no_message(pareto_khat(tau, extra_diags = TRUE, verbose = FALSE)) }) @@ -108,13 +109,16 @@ test_that("pareto_khat diagnostics handle special cases correctly", { set.seed(1234) x <- c(rnorm(50), NA) - expect_true(is.na(pareto_khat(x))) + expect_warning(pareto_khat(x)) + expect_true(is.na(suppressWarnings(pareto_khat(x)))) x <- c(rnorm(50), Inf) - expect_true(is.na(pareto_khat(x))) + expect_warning(pareto_khat(x)) + expect_true(is.na(suppressWarnings(pareto_khat(x)))) x <- rep(1, 50) - expect_true(is.na(pareto_khat(x))) + expect_warning(pareto_khat(x)) + expect_true(is.na(suppressWarnings(pareto_khat(x)))) }) From c383f500111c5968596177502c29abbb5b1f19e1 Mon Sep 17 00:00:00 2001 From: n-kall Date: Wed, 24 May 2023 09:40:34 +0300 Subject: [PATCH 27/30] expand pareto diagnostic documentation --- R/pareto_smooth.R | 62 +++++++++++++++++++++++++++++++------------- man/pareto_khat.Rd | 6 +++-- man/pareto_smooth.Rd | 11 ++++---- 3 files changed, 53 insertions(+), 26 deletions(-) diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 3418cf88..735775b3 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -1,7 +1,9 @@ #' Pareto khat diagnostic #' -#' Estimate Pareto k value by fitting a Generalized Pareto Distribution to one or two tails of x. For -#' further details see Vehtari et al. (2022). +#' Estimate Pareto k value by fitting a Generalized Pareto +#' Distribution to one or two tails of x. This can be used to etimate +#' the number of fractional moments. For further details see Vehtari +#' et al. (2022). #' #' @template args-pareto #' @template args-methods-dots @@ -61,8 +63,7 @@ pareto_khat.rvar <- function(x, ...) { #' #' Compute diagnostics for Pareto smoothing the tail draws of x by #' replacing tail draws by order statistics of a generalized Pareto -#' distribution fit to the tail(s). For further details see Vehtari et -#' al. (2022). +#' distribution fit to the tail(s). #' #' @template args-pareto #' @template args-methods-dots @@ -73,6 +74,32 @@ pareto_khat.rvar <- function(x, ...) { #' * `khat_threshold`: khat-threshold for reliable Pareto smoothed estimate, #' * `convergence_rate`: Pareto smoothed estimate RMSE convergence rate. #' +#' @details When the fitted Generalized Pareto Distribution is used to +#' smooth the tail values and these smoothed values are used to +#' compute expectations, the following diagnostics can give further +#' information about the reliability of these estimates. +#' +#' * `min_ss`: Minimum sample size for reliable Pareto smoothed +#' estimate. If the actual sample size is greater than `min_ss`, then +#' Pareto smoothed estimates can be considered reliable. If the actual +#' sample size is lower than `min_ss`, increasing the sample size +#' might result in more reliable estimates. For further details, see +#' Section 3.2.3, Equation 11 in Vehtari et al. (2022). +#' +#' * `khat_threshold`: Threshold below which k-hat values result in +#' reliable Pareto smoothed estimates. The threshold is lower for +#' smaller effective sample sizes. If k-hat is larger than the +#' threshold, increasing the total sample size may improve reliability +#' of estimates. For further details, see Section 3.2.4, Equation 13 +#' in Vehtari et al. (2022). +#' +#' * `convergence_rate`: Relative convergence rate compared to the +#' central limit theorem. Applicable only if the actual sample size +#' is sufficiently large (greater than `min_ss`). The convergence +#' rate tells the rate at which the variance of an estimate reduces +#' when the sample size is increased, compared to the central limit +#' theorem convergence rate. See Appendix B in Vehtari et al. (2022). +#' #' @examples #' mu <- extract_variable_matrix(example_draws(), "mu") #' pareto_diags(mu) @@ -90,7 +117,7 @@ pareto_diags.default <- function(x, r_eff = NULL, ndraws_tail = NULL, verbose = FALSE, - ...) { + ...) { smoothed <- pareto_smooth.default( x, @@ -102,9 +129,9 @@ pareto_diags.default <- function(x, verbose = verbose, smooth_draws = FALSE, ...) - + return(smoothed$diagnostics) - + } @@ -119,7 +146,7 @@ pareto_diags.rvar <- function(x, ...) { extra_diags = TRUE, ... ) - + dim(draws_diags) <- dim(draws_diags) %||% length(draws_diags) margins <- seq_along(dim(draws_diags)) @@ -147,21 +174,20 @@ pareto_diags.rvar <- function(x, ...) { #' draws and diagnostics. Default is `TRUE`. #' @param extra_diags (logical) Should extra Pareto khat diagnostics #' be included in output? If `TRUE`, `min_ss`, `khat_threshold` and -#' `convergence_rate` for the calculated k value will be +#' `convergence_rate` for the estimated k value will be #' returned. Default is `FALSE`. #' @template args-methods-dots #' @template ref-vehtari-paretosmooth-2022 -#' @return Either a vector `x` of smoothed values and or a list -#' containing the vector `x` and `diagnostics` of Pareto smoothing +#' @return Either a vector `x` of smoothed values or a named list +#' containing the vector `x` and a named list `diagnostics` containing Pareto smoothing #' diagnostics: #' * `khat`: estimated Pareto k shape parameter, and #' optionally #' * `min_ss`: minimum sample size for reliable Pareto #' smoothed estimate #' * `khat_threshold`: khat-threshold for reliable -#' Pareto smoothed estimat -#' * `convergence_rate`: Pareto smoothed -#' estimate RMSE convergence rate. +#' Pareto smoothed estimates +#' * `convergence_rate`: Relative convergence rate for Pareto smoothed estimates #' #' @examples #' mu <- extract_variable_matrix(example_draws(), "mu") @@ -179,14 +205,14 @@ pareto_smooth.rvar <- function(x, return_k = TRUE, extra_diags = FALSE, ...) { if (extra_diags) { return_k <- TRUE } - + draws_diags <- summarise_rvar_by_element_with_chains(x, pareto_smooth.default, return_k = return_k, extra_diags = extra_diags, ...) dim(draws_diags) <- dim(draws_diags) %||% length(draws_diags) margins <- seq_along(dim(draws_diags)) if (return_k) { if (extra_diags) { - + diags <- list( khat = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat), min_ss = apply(draws_diags, margins, function(x) x[[1]]$diagnostics$min_ss), @@ -217,14 +243,14 @@ pareto_smooth.default <- function(x, return_k = TRUE, extra_diags = FALSE, verbose = 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) - + # check for infinite or na values if (should_return_NA(x)) { warning_no_call("Input contains infinite or NA values, Pareto smoothing not performed.") diff --git a/man/pareto_khat.Rd b/man/pareto_khat.Rd index 2c3a148b..7a8eb780 100644 --- a/man/pareto_khat.Rd +++ b/man/pareto_khat.Rd @@ -55,8 +55,10 @@ printed. Default is \code{FALSE}.} \code{khat} estimated Generalized Pareto Distribution shape parameter k } \description{ -Estimate Pareto k value by fitting a Generalized Pareto Distribution to one or two tails of x. For -further details see Vehtari et al. (2022). +Estimate Pareto k value by fitting a Generalized Pareto +Distribution to one or two tails of x. This can be used to etimate +the number of fractional moments. For further details see Vehtari +et al. (2022). } \examples{ mu <- extract_variable_matrix(example_draws(), "mu") diff --git a/man/pareto_smooth.Rd b/man/pareto_smooth.Rd index 52d085d0..e3d629c5 100644 --- a/man/pareto_smooth.Rd +++ b/man/pareto_smooth.Rd @@ -37,7 +37,7 @@ draws and diagnostics. Default is \code{TRUE}.} \item{extra_diags}{(logical) Should extra Pareto khat diagnostics be included in output? If \code{TRUE}, \code{min_ss}, \code{khat_threshold} and -\code{convergence_rate} for the calculated k value will be +\code{convergence_rate} for the estimated k value will be returned. Default is \code{FALSE}.} \item{tail}{(string) The tail to diagnose/smooth: @@ -63,8 +63,8 @@ length(x) / 5 otherwise (see Appendix H in Vehtari et al. (2022)).} printed. Default is \code{FALSE}.} } \value{ -Either a vector \code{x} of smoothed values and or a list -containing the vector \code{x} and \code{diagnostics} of Pareto smoothing +Either a vector \code{x} of smoothed values or a named list +containing the vector \code{x} and a named list \code{diagnostics} containing Pareto smoothing diagnostics: \itemize{ \item \code{khat}: estimated Pareto k shape parameter, and @@ -72,9 +72,8 @@ optionally \item \code{min_ss}: minimum sample size for reliable Pareto smoothed estimate \item \code{khat_threshold}: khat-threshold for reliable -Pareto smoothed estimat -\item \code{convergence_rate}: Pareto smoothed -estimate RMSE convergence rate. +Pareto smoothed estimates +\item \code{convergence_rate}: Relative convergence rate for Pareto smoothed estimates } } \description{ From 4af4ede9cdc246312f633e894d029d05057b007a Mon Sep 17 00:00:00 2001 From: n-kall Date: Thu, 25 May 2023 12:24:03 +0300 Subject: [PATCH 28/30] add further tests for pareto functions --- tests/testthat/test-pareto_smooth.R | 40 +++++++++++++++++++++-------- 1 file changed, 29 insertions(+), 11 deletions(-) diff --git a/tests/testthat/test-pareto_smooth.R b/tests/testthat/test-pareto_smooth.R index 28b16a4a..03e06f27 100644 --- a/tests/testthat/test-pareto_smooth.R +++ b/tests/testthat/test-pareto_smooth.R @@ -32,7 +32,7 @@ test_that("pareto_khat handles ndraws_tail argument", { expect_warning(pareto_khat(tau, tail = "both", ndraws_tail = 4), "Number of tail draws cannot be less than 5. Changing to 5.") - + }) @@ -42,7 +42,7 @@ test_that("pareto_khat handles r_eff argument", { pk1 <- pareto_khat(tau, r_eff = 1) pk0.6 <- pareto_khat(tau, r_eff = 0.6) expect_true(pk1$khat < pk0.6$khat) - + }) @@ -59,7 +59,7 @@ test_that("pareto_khat diagnostics messages are as expected", { paste0('To halve the RMSE, approximately 4.1 times bigger S is needed.')) diags$khat <- 0.6 - + expect_message(pareto_k_diagmsg(diags), paste0('S is too small, and sample size larger than 10 is neeeded for reliable results.\n')) @@ -71,7 +71,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', 'further draws may improve the estimates, but it is not possible to predict\n', @@ -79,7 +79,7 @@ test_that("pareto_khat diagnostics messages are as expected", { }) -test_that("pareto_khat extra diags are returned", { +test_that("pareto_diags returns expected diagnostics", { tau <- extract_variable_matrix(example_draws(), "tau") @@ -127,7 +127,7 @@ test_that("pareto_khat functions work with matrix with chains", { tau_chains <- extract_variable_matrix(example_draws(), "tau") tau_nochains <- extract_variable(example_draws(), "tau") - + expect_equal(pareto_khat(tau_chains, ndraws_tail = 20), pareto_khat(tau_nochains, ndraws_tail = 20)) @@ -142,7 +142,7 @@ test_that("pareto_khat functions work with rvars with and without chains", { tau_chains <- extract_variable_matrix(example_draws(), "tau") tau_nochains <- extract_variable(example_draws(), "tau") - + tau_rvar_chains <- rvar(tau_chains, with_chains = TRUE) tau_rvar_nochains <- rvar(tau_nochains) @@ -152,16 +152,34 @@ test_that("pareto_khat functions work with rvars with and without chains", { expect_equal(pareto_khat(tau_rvar_chains, ndraws_tail = 20), pareto_khat(tau_rvar_nochains, ndraws_tail = 20)) + + expect_equal(pareto_diags(tau_chains, ndraws_tail = 20), + pareto_diags(tau_rvar_chains, ndraws_tail = 20)) + + expect_equal(pareto_diags(tau_rvar_chains, ndraws_tail = 20), + pareto_diags(tau_rvar_nochains, ndraws_tail = 20)) + ps_chains <- pareto_smooth(tau_chains, ndraws_tail = 20) ps_rvar_chains <- pareto_smooth(tau_rvar_chains, ndraws_tail = 20) ps_nochains <- pareto_smooth(tau_nochains, ndraws_tail = 20) ps_rvar_nochains <- pareto_smooth(tau_rvar_nochains, ndraws_tail = 20) - + expect_equal(ps_rvar_chains$x, rvar(ps_chains$x, with_chains = TRUE)) - + + expect_equal(ps_rvar_nochains$x, rvar(ps_nochains$x)) + + + ps_chains <- pareto_smooth(tau_chains, ndraws_tail = 20, extra_diags = TRUE) + ps_rvar_chains <- pareto_smooth(tau_rvar_chains, ndraws_tail = 20, extra_diags = TRUE) + + ps_nochains <- pareto_smooth(tau_nochains, ndraws_tail = 20, extra_diags = TRUE) + ps_rvar_nochains <- pareto_smooth(tau_rvar_nochains, ndraws_tail = 20, extra_diags = TRUE) + + expect_equal(ps_rvar_chains$x, rvar(ps_chains$x, with_chains = TRUE)) + expect_equal(ps_rvar_nochains$x, rvar(ps_nochains$x)) - + }) test_that("pareto_smooth returns x with smoothed tail", { @@ -172,5 +190,5 @@ test_that("pareto_smooth returns x with smoothed tail", { expect_equal(sort(tau)[1:390], sort(tau_smoothed)[1:390]) expect_false(isTRUE(all.equal(sort(tau), sort(tau_smoothed)))) - + }) From 54011ec16ad503ff314801f6c86ad9855d35396d Mon Sep 17 00:00:00 2001 From: n-kall Date: Thu, 25 May 2023 13:22:53 +0300 Subject: [PATCH 29/30] add documentation file fo pareto_diags --- man/pareto_diags.Rd | 105 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 105 insertions(+) create mode 100644 man/pareto_diags.Rd diff --git a/man/pareto_diags.Rd b/man/pareto_diags.Rd new file mode 100644 index 00000000..6513c690 --- /dev/null +++ b/man/pareto_diags.Rd @@ -0,0 +1,105 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pareto_smooth.R +\name{pareto_diags} +\alias{pareto_diags} +\alias{pareto_diags.default} +\alias{pareto_diags.rvar} +\title{Pareto smoothing diagnostics} +\usage{ +pareto_diags(x, ...) + +\method{pareto_diags}{default}( + x, + tail = c("both", "right", "left"), + r_eff = NULL, + ndraws_tail = NULL, + verbose = FALSE, + ... +) + +\method{pareto_diags}{rvar}(x, ...) +} +\arguments{ +\item{x}{(multiple options) One of: +\itemize{ +\item A matrix of draws for a single variable (iterations x chains). See +\code{\link[=extract_variable_matrix]{extract_variable_matrix()}}. +\item An \code{\link{rvar}}. +}} + +\item{...}{Arguments passed to individual methods (if applicable).} + +\item{tail}{(string) The tail to diagnose/smooth: +\itemize{ +\item \code{"right"}: diagnose/smooth only the right (upper) tail +\item \code{"left"}: diagnose/smooth only the left (lower) tail +\item \code{"both"}: diagnose/smooth both tails and return the maximum k value +} + +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.} + +\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)).} + +\item{verbose}{(logical) Should diagnostic messages be printed? If +\code{TRUE}, messages related to Pareto diagnostics will be +printed. Default is \code{FALSE}.} +} +\value{ +List of Pareto smoothing diagnostics: +\itemize{ +\item \code{khat}: estimated Pareto k shape parameter, +\item \code{min_ss}: minimum sample size for reliable Pareto smoothed estimate, +\item \code{khat_threshold}: khat-threshold for reliable Pareto smoothed estimate, +\item \code{convergence_rate}: Pareto smoothed estimate RMSE convergence rate. +} +} +\description{ +Compute diagnostics for Pareto smoothing the tail draws of x by +replacing tail draws by order statistics of a generalized Pareto +distribution fit to the tail(s). +} +\details{ +When the fitted Generalized Pareto Distribution is used to +smooth the tail values and these smoothed values are used to +compute expectations, the following diagnostics can give further +information about the reliability of these estimates. +\itemize{ +\item \code{min_ss}: Minimum sample size for reliable Pareto smoothed +estimate. If the actual sample size is greater than \code{min_ss}, then +Pareto smoothed estimates can be considered reliable. If the actual +sample size is lower than \code{min_ss}, increasing the sample size +might result in more reliable estimates. For further details, see +Section 3.2.3, Equation 11 in Vehtari et al. (2022). +\item \code{khat_threshold}: Threshold below which k-hat values result in +reliable Pareto smoothed estimates. The threshold is lower for +smaller effective sample sizes. If k-hat is larger than the +threshold, increasing the total sample size may improve reliability +of estimates. For further details, see Section 3.2.4, Equation 13 +in Vehtari et al. (2022). +\item \code{convergence_rate}: Relative convergence rate compared to the +central limit theorem. Applicable only if the actual sample size +is sufficiently large (greater than \code{min_ss}). The convergence +rate tells the rate at which the variance of an estimate reduces +when the sample size is increased, compared to the central limit +theorem convergence rate. See Appendix B in Vehtari et al. (2022). +} +} +\examples{ +mu <- extract_variable_matrix(example_draws(), "mu") +pareto_diags(mu) + +d <- as_draws_rvars(example_draws("multi_normal")) +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 +} From a85b9b58b4dfe426eb679406f48a762208e36555 Mon Sep 17 00:00:00 2001 From: n-kall Date: Wed, 30 Aug 2023 15:57:43 +0300 Subject: [PATCH 30/30] use ess_tail for r_eff, documentation fixes, typofixes --- R/gpd.R | 6 +++--- R/pareto_smooth.R | 14 +++++++------- man-roxygen/args-pareto.R | 2 +- man/pareto_diags.Rd | 2 +- man/pareto_khat.Rd | 8 ++++---- man/pareto_smooth.Rd | 2 +- man/set_variables.Rd | 2 +- tests/testthat/test-pareto_smooth.R | 6 +++--- 8 files changed, 21 insertions(+), 21 deletions(-) diff --git a/R/gpd.R b/R/gpd.R index a893dd38..924f8c3f 100644 --- a/R/gpd.R +++ b/R/gpd.R @@ -23,10 +23,10 @@ qgeneralized_pareto <- function(p, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, #' #' Given a sample \eqn{x}, Estimate the parameters \eqn{k} and \eqn{\sigma} of #' the generalized Pareto distribution (GPD), assuming the location parameter is -#' 0. By default the fit uses a prior for \eqn{k}, which will stabilize +#' 0. By default the fit uses a prior for \eqn{k} (this is in addition to the prior described by Zhang and Stephens, 2009), which will stabilize #' estimates for very small sample sizes (and low effective sample sizes in the #' case of MCMC samples). The weakly informative prior is a Gaussian prior -#' centered at 0.5. +#' centered at 0.5 (see details in Vehtari et al., 2022). #' #' @param x A numeric vector. The sample from which to estimate the parameters. #' @param wip Logical indicating whether to adjust \eqn{k} based on a weakly @@ -67,7 +67,7 @@ gpdfit <- function(x, wip = TRUE, min_grid_pts = 30, sort_x = TRUE) { sigma_hat <- -k_hat / theta_hat # adjust k_hat based on weakly informative prior, Gaussian centered on 0.5. - # this stabilizes estimates for very small Monte Carlo sample sizes and low neff + # this stabilizes estimates for very small Monte Carlo sample sizes and low neff (see Vehtari et al., 2022 for details) if (wip) { k_hat <- (k_hat * N + 0.5 * 10) / (N + 10) } diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R index 735775b3..228aede9 100644 --- a/R/pareto_smooth.R +++ b/R/pareto_smooth.R @@ -1,9 +1,9 @@ #' Pareto khat diagnostic #' #' Estimate Pareto k value by fitting a Generalized Pareto -#' Distribution to one or two tails of x. This can be used to etimate -#' the number of fractional moments. For further details see Vehtari -#' et al. (2022). +#' Distribution to one or two tails of x. This can be used to estimate +#' the number of fractional moments that is useful for convergence +#' diagnostics. For further details see Vehtari et al. (2022). #' #' @template args-pareto #' @template args-methods-dots @@ -262,7 +262,7 @@ pareto_smooth.default <- function(x, # automatically calculate relative efficiency if (is.null(r_eff)) { - r_eff <- ess_basic(x) / S + r_eff <- ess_tail(x) / S } # automatically calculate tail length @@ -274,7 +274,7 @@ pareto_smooth.default <- function(x, if (ndraws_tail > S / 2) { warning("Number of tail draws cannot be more than half ", - "the length of the draws if both tails are fit, ", + "the total number of draws if both tails are fit, ", "changing to ", S / 2, ".") ndraws_tail <- S / 2 } @@ -416,7 +416,7 @@ pareto_smooth.default <- function(x, #' Extra pareto-k diagnostics #' -#' internal function to calculate the extra diagnostics or a given +#' internal function to calculate the extra diagnostics for a given #' pareto k and sample size S #' @noRd .pareto_smooth_extra_diags <- function(k, S, ...) { @@ -521,7 +521,7 @@ pareto_k_diagmsg <- function(diags, ...) { '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 neeeded 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') } diff --git a/man-roxygen/args-pareto.R b/man-roxygen/args-pareto.R index 51cae9c3..a406dbde 100644 --- a/man-roxygen/args-pareto.R +++ b/man-roxygen/args-pareto.R @@ -5,7 +5,7 @@ #' @param tail (string) The tail to diagnose/smooth: #' * `"right"`: diagnose/smooth only the right (upper) tail #' * `"left"`: diagnose/smooth only the left (lower) tail -#' * `"both"`: diagnose/smooth both tails and return the maximum k value +#' * `"both"`: diagnose/smooth both tails and return the maximum k-hat value #' #' The default is `"both"`. #' @param ndraws_tail (numeric) number of draws for the tail. If diff --git a/man/pareto_diags.Rd b/man/pareto_diags.Rd index 6513c690..c2558a9a 100644 --- a/man/pareto_diags.Rd +++ b/man/pareto_diags.Rd @@ -33,7 +33,7 @@ pareto_diags(x, ...) \itemize{ \item \code{"right"}: diagnose/smooth only the right (upper) tail \item \code{"left"}: diagnose/smooth only the left (lower) tail -\item \code{"both"}: diagnose/smooth both tails and return the maximum k value +\item \code{"both"}: diagnose/smooth both tails and return the maximum k-hat value } The default is \code{"both"}.} diff --git a/man/pareto_khat.Rd b/man/pareto_khat.Rd index 7a8eb780..c3383eb4 100644 --- a/man/pareto_khat.Rd +++ b/man/pareto_khat.Rd @@ -33,7 +33,7 @@ pareto_khat(x, ...) \itemize{ \item \code{"right"}: diagnose/smooth only the right (upper) tail \item \code{"left"}: diagnose/smooth only the left (lower) tail -\item \code{"both"}: diagnose/smooth both tails and return the maximum k value +\item \code{"both"}: diagnose/smooth both tails and return the maximum k-hat value } The default is \code{"both"}.} @@ -56,9 +56,9 @@ printed. Default is \code{FALSE}.} } \description{ Estimate Pareto k value by fitting a Generalized Pareto -Distribution to one or two tails of x. This can be used to etimate -the number of fractional moments. For further details see Vehtari -et al. (2022). +Distribution to one or two tails of x. This can be used to estimate +the number of fractional moments that is useful for convergence +diagnostics. For further details see Vehtari et al. (2022). } \examples{ mu <- extract_variable_matrix(example_draws(), "mu") diff --git a/man/pareto_smooth.Rd b/man/pareto_smooth.Rd index e3d629c5..259421de 100644 --- a/man/pareto_smooth.Rd +++ b/man/pareto_smooth.Rd @@ -44,7 +44,7 @@ returned. Default is \code{FALSE}.} \itemize{ \item \code{"right"}: diagnose/smooth only the right (upper) tail \item \code{"left"}: diagnose/smooth only the left (lower) tail -\item \code{"both"}: diagnose/smooth both tails and return the maximum k value +\item \code{"both"}: diagnose/smooth both tails and return the maximum k-hat value } The default is \code{"both"}.} diff --git a/man/set_variables.Rd b/man/set_variables.Rd index fbbccc78..ba67f440 100644 --- a/man/set_variables.Rd +++ b/man/set_variables.Rd @@ -25,7 +25,7 @@ when using pipe operators. x <- as_draws(matrix(rnorm(100), ncol = 2)) variables(x) -set_variables(x, c("theta[1]", "theta[2]")) +x <- set_variables(x, c("theta[1]", "theta[2]")) variables(x) # this is equivalent to diff --git a/tests/testthat/test-pareto_smooth.R b/tests/testthat/test-pareto_smooth.R index 03e06f27..dff22ee1 100644 --- a/tests/testthat/test-pareto_smooth.R +++ b/tests/testthat/test-pareto_smooth.R @@ -27,8 +27,8 @@ test_that("pareto_khat handles ndraws_tail argument", { expect_warning(pareto_khat(tau, tail = "both", ndraws_tail = 201), "Number of tail draws cannot be more than half ", - "the length of the draws if both tails are fit, ", - "changing to 200.") + "the total number of draws if both tails are fit, ", + "changing to 200.") expect_warning(pareto_khat(tau, tail = "both", ndraws_tail = 4), "Number of tail draws cannot be less than 5. Changing to 5.") @@ -61,7 +61,7 @@ test_that("pareto_khat diagnostics messages are as expected", { diags$khat <- 0.6 expect_message(pareto_k_diagmsg(diags), - paste0('S is too small, and sample size larger than 10 is neeeded for reliable results.\n')) + paste0('S is too small, and sample size larger than 10 is needed for reliable results.\n')) diags$khat <- 0.71 diags$khat_threshold <- 0.8