diff --git a/.github/workflows/covr.yml b/.github/workflows/covr.yml
index 72f46b83..c3724c93 100644
--- a/.github/workflows/covr.yml
+++ b/.github/workflows/covr.yml
@@ -35,7 +35,7 @@ jobs:
shell: Rscript {0}
- name: Cache R packages
- uses: actions/cache@v2
+ uses: actions/cache@v4
with:
path: ${{ env.R_LIBS_USER }}
key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }}
diff --git a/.github/workflows/rcmdcheck.yml b/.github/workflows/rcmdcheck.yml
index a380a7fa..bf665d2a 100644
--- a/.github/workflows/rcmdcheck.yml
+++ b/.github/workflows/rcmdcheck.yml
@@ -58,7 +58,7 @@ jobs:
- name: Cache R packages
if: runner.os != 'Windows'
- uses: actions/cache@v2
+ uses: actions/cache@v4
with:
path: ${{ env.R_LIBS_USER }}
key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }}
@@ -98,7 +98,7 @@ jobs:
- name: Upload check results
if: failure()
- uses: actions/upload-artifact@main
+ uses: actions/upload-artifact@v4
with:
name: ${{ runner.os }}-r${{ matrix.config.r }}-results
path: check
diff --git a/DESCRIPTION b/DESCRIPTION
index ce462430..b898582b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
Package: posterior
Title: Tools for Working with Posterior Distributions
-Version: 1.5.0.9000
-Date: 2023-10-31
+Version: 1.6.1
+Date: 2025-02-27
Authors@R: c(person("Paul-Christian", "Bürkner", email = "paul.buerkner@gmail.com", role = c("aut", "cre")),
person("Jonah", "Gabry", email = "jsg2201@columbia.edu", role = c("aut")),
person("Matthew", "Kay", email = "mjskay@northwestern.edu", role = c("aut")),
@@ -11,7 +11,8 @@ Authors@R: c(person("Paul-Christian", "Bürkner", email = "paul.buerkner@gmail.c
person("Ben", "Lambert", role = c("ctb")),
person("Ozan", "Adıgüzel", role = c("ctb")),
person("Jacob", "Socolar", role = c("ctb")),
- person("Noa", "Kallioinen", role = c("ctb")))
+ person("Noa", "Kallioinen", role = c("ctb")),
+ person("Teemu", "Säilynoja", role = c("ctb")))
Description: Provides useful tools for both users and developers of packages
for fitting Bayesian models or working with output from Bayesian models.
The primary goals of the package are to:
@@ -56,5 +57,5 @@ LazyData: false
URL: https://mc-stan.org/posterior/, https://discourse.mc-stan.org/
BugReports: https://github.com/stan-dev/posterior/issues
Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.0
+RoxygenNote: 7.3.2
VignetteBuilder: knitr
diff --git a/NAMESPACE b/NAMESPACE
index f8497c96..7c721a6b 100755
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -245,6 +245,9 @@ S3method(pareto_min_ss,rvar)
S3method(pareto_smooth,default)
S3method(pareto_smooth,rvar)
S3method(pillar_shaft,rvar)
+S3method(pit,default)
+S3method(pit,draws_matrix)
+S3method(pit,rvar)
S3method(print,draws_array)
S3method(print,draws_df)
S3method(print,draws_list)
@@ -480,6 +483,11 @@ export(pareto_khat)
export(pareto_khat_threshold)
export(pareto_min_ss)
export(pareto_smooth)
+export(pit)
+export(ps_convergence_rate)
+export(ps_khat_threshold)
+export(ps_min_ss)
+export(ps_tail_length)
export(quantile2)
export(r_scale)
export(rdo)
diff --git a/NEWS.md b/NEWS.md
index 1d0e079d..2fae4868 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,11 +1,33 @@
-# posterior (development version)
+# posterior 1.6.1
+
+### Bug Fixes
+
+* Fix a test issue that led to an R CMD check failure on R devel.
### Enhancements
-* Add `pareto_smooth` option to `weight_draws`, to Pareto smooth
- weights before adding to a draws object.
+* Convert lists of matrices to `draws_array` objects.
+* Improve the documentation in various places.
+* Implement `pit()` for draws and rvar objects. LOO-PIT can be computed using `weights`.
* Add support for applying weights to individual `rvar` objects.
* Add `log_weights()` function for easy access to raw internal weights.
+* Add support for weights in `ess` and `mcse` functions.
+
+# posterior 1.6.0
+
+### Enhancements
+
+* Add `exclude` option to `subset_draws()`, which can be used to exclude
+ the matched selection.
+* Add `are_log_weights` option to `pareto_smooth()`, which is necessary
+ for correct Pareto smoothing computation if the input vector
+ consists of log weights.
+* Add `pareto_smooth` option to `weight_draws()`, to Pareto smooth
+ weights before adding to a draws object.
+* Add individual Pareto diagnostic functions (`pareto_khat()`,
+ `pareto_khat_threshold()`, `pareto_min_ss()`, `pareto_convergence_rate()`)
+* `thin_draws()` now automatically thins draws based on ESS by default,
+ and non-integer thinning is possible.
* Matrix multiplication of `rvar`s can now be done with the base matrix
multiplication operator (`%*%`) instead of `%**%` in R >= 4.3.
* `variables()`, `variables<-()`, `set_variables()`, and `nvariables()` now
diff --git a/R/as_draws.R b/R/as_draws.R
index c52e75a5..84c4e669 100644
--- a/R/as_draws.R
+++ b/R/as_draws.R
@@ -76,10 +76,11 @@ closest_draws_format <- function(x) {
out <- "rvars"
} else if (is_draws_list_like(x)) {
out <- "list"
- }
- else {
- stop_no_call("Don't know how to transform an object of class ",
- "'", class(x)[1L], "' to any supported draws format.")
+ } else {
+ stop_no_call(
+ "Don't know how to transform an object of class '",
+ class(x)[1L], "' to any supported draws format."
+ )
}
paste0("draws_", out)
}
diff --git a/R/as_draws_array.R b/R/as_draws_array.R
index 84f4ec3a..11f77308 100644
--- a/R/as_draws_array.R
+++ b/R/as_draws_array.R
@@ -131,7 +131,11 @@ as_draws_array.mcmc.list <- function(x, ...) {
# try to convert any R object into a 'draws_array' object
.as_draws_array <- function(x) {
- x <- as.array(x)
+ if (is_matrix_list_like(x)) {
+ x <- as_array_matrix_list(x)
+ } else {
+ x <- as.array(x)
+ }
new_dimnames <- list(iteration = NULL, chain = NULL, variable = NULL)
if (!is.null(dimnames(x)[[3]])) {
new_dimnames[[3]] <- dimnames(x)[[3]]
@@ -178,7 +182,14 @@ is_draws_array <- function(x) {
# is an object looking like a 'draws_array' object?
is_draws_array_like <- function(x) {
- is.array(x) && length(dim(x)) == 3L
+ is.array(x) && length(dim(x)) == 3L ||
+ is_matrix_list_like(x)
+}
+
+# is an object likely a list of matrices?
+# such an object can be easily converted to a draws_array
+is_matrix_list_like <- function(x) {
+ is.list(x) && length(dim(x[[1]])) == 2L
}
#' Extract parts of a `draws_array` object
@@ -217,15 +228,8 @@ variance.draws_array <- function(x, ...) {
# convert a list of matrices to an array
as_array_matrix_list <- function(x) {
stopifnot(is.list(x))
- if (length(x) == 1) {
- tmp <- dimnames(x[[1]])
- x <- x[[1]]
- dim(x) <- c(dim(x), 1)
- dimnames(x) <- tmp
- } else {
- x <- abind::abind(x, along = 3L)
- }
- x <- aperm(x, c(1, 3, 2))
+ x <- abind::abind(x, along = 3L)
+ aperm(x, c(1, 3, 2))
}
# create an empty draws_array object
@@ -246,4 +250,3 @@ empty_draws_array <- function(variables = character(0), nchains = 0,
class(out) <- class_draws_array()
out
}
-
diff --git a/R/convergence.R b/R/convergence.R
index 0fb6e95f..c0ca11ed 100644
--- a/R/convergence.R
+++ b/R/convergence.R
@@ -16,11 +16,14 @@
#' | [ess_basic()] | Basic version of effective sample size |
#' | [ess_bulk()] | Bulk effective sample size |
#' | [ess_tail()] | Tail effective sample size |
+#' | [ess_mean()] | Effective sample sizes for the mean |
+#' | [ess_median()] | Effective sample sizes for the median |
#' | [ess_quantile()] | Effective sample sizes for quantiles |
-#' | [ess_sd()] | Effective sample sizes for standard deviations |
+#' | [ess_sd()] | Effective sample sizes for the standard deviation |
#' | [mcse_mean()] | Monte Carlo standard error for the mean |
+#' | [mcse_median()] | Monte Carlo standard error for the median |
#' | [mcse_quantile()] | Monte Carlo standard error for quantiles |
-#' | [mcse_sd()] | Monte Carlo standard error for standard deviations |
+#' | [mcse_sd()] | Monte Carlo standard error for the standard deviation |
#' | [pareto_khat()] | Pareto khat diagnostic for tail(s) |
#' | [pareto_diags()] | Additional diagnostics related to Pareto khat |
#' | [rhat_basic()] | Basic version of Rhat |
@@ -88,6 +91,7 @@ rhat_basic.rvar <- function(x, split = TRUE, ...) {
#' @family diagnostics
#' @template args-conv
#' @template args-conv-split
+#' @template args-conv-weights
#' @template args-methods-dots
#' @template return-conv
#' @template ref-gelman-bda-2013
@@ -178,6 +182,7 @@ rhat.rvar <- function(x, ...) {
#'
#' @family diagnostics
#' @template args-conv
+#' @template args-conv-weights
#' @template args-methods-dots
#' @template return-conv
#' @template ref-vehtari-rhat-2021
@@ -231,6 +236,7 @@ ess_bulk.rvar <- function(x, ...) {
#'
#' @family diagnostics
#' @template args-conv
+#' @template args-conv-weights
#' @template args-methods-dots
#' @template return-conv
#' @template ref-vehtari-rhat-2021
@@ -270,6 +276,7 @@ ess_tail.rvar <- function(x, ...) {
#' @family diagnostics
#' @template args-conv
#' @template args-conv-quantile
+#' @template args-conv-weights
#' @template args-methods-dots
#' @template return-conv-quantile
#' @template ref-vehtari-rhat-2021
@@ -360,6 +367,7 @@ ess_median <- function(x, ...) {
#' estimate of a single variable.
#'
#' @template args-conv
+#' @template args-conv-weights
#' @template args-methods-dots
#' @template return-conv
#' @template ref-gelman-bda-2013
@@ -374,7 +382,7 @@ ess_median <- function(x, ...) {
#' @export
ess_mean <- function(x, ...) UseMethod("ess_mean")
-#' @rdname ess_quantile
+#' @rdname ess_mean
#' @export
ess_mean.default <- function(x, weights = NULL, ...) {
@@ -410,6 +418,7 @@ ess_mean.rvar <- function(x, ...) {
#'
#' @family diagnostics
#' @template args-conv
+#' @template args-conv-weights
#' @template args-methods-dots
#' @template return-conv
#' @template ref-vehtari-rhat-2021
@@ -428,6 +437,9 @@ ess_sd <- function(x, ...) UseMethod("ess_sd")
#' @export
ess_sd.default <- function(x, weights = NULL, ...) {
if (is.null(weights)) {
+ # var/sd are not a simple expectation of g(X), e.g. variance
+ # has (X-E[X])^2. The following ESS is based on a relevant quantity
+ # in the computation and is empirically a good choice.
.ess(.split_chains(abs(x - mean(x))))
} else {
@@ -449,7 +461,7 @@ ess_sd.rvar <- function(x, ...) {
summarise_rvar_by_element_with_chains(x, ess_sd, weights = weights, ...)
}
-# TODO: ess_weights
+# TODO: ess_weights function
#' Monte Carlo standard error for quantiles
#'
@@ -460,6 +472,7 @@ ess_sd.rvar <- function(x, ...) {
#' @family diagnostics
#' @template args-conv
#' @template args-conv-quantile
+#' @template args-conv-weights
#' @template args-methods-dots
#' @template return-conv-quantile
#' @template ref-vehtari-rhat-2021
@@ -555,6 +568,7 @@ mcse_median <- function(x, ...) {
#'
#' @family diagnostics
#' @template args-conv
+#' @template args-conv-weights
#' @template args-methods-dots
#' @template return-conv
#' @template ref-gelman-bda-2013
@@ -583,7 +597,7 @@ mcse_mean.default <- function(x, weights = NULL, ...) {
x <- as.matrix(x)
r_eff <- .ess(.split_chains(x)) / (nrow(x) * ncol(x))
- .mcse_weighted(x, weights, r_eff, ...)
+ .mcse_weighted(x, weights = weights, r_eff = r_eff, ...)
}
}
@@ -602,6 +616,7 @@ mcse_mean.rvar <- function(x, ...) {
#'
#' @family diagnostics
#' @template args-conv
+#' @template args-conv-weights
#' @template args-methods-dots
#' @template return-conv
#' @template ref-vehtari-rhat-2021
@@ -632,7 +647,6 @@ mcse_sd.default <- function(x, weights = NULL, ...) {
# which doesn't assume normality of sims.
Evar <- mean(sims_c^2)
varvar <- (mean(sims_c^4) - Evar^2) / ess # (Equation 6.20)
-
# The first order Taylor series approximation of variance of sd.
# Kenney and Keeping (1951, p. 141) write "...since fluctuations of
# any moment are of order N^{-1/2}, squares and higher powers of
@@ -656,7 +670,7 @@ mcse_sd.default <- function(x, weights = NULL, ...) {
second_moment_weighted <- weighted.mean(x_centered^2, w = weights)
fourth_moment_weighted <- weighted.mean(x_centered^4, w = weights)
- r_eff <- .ess(x_centered^2) / (nrow(x) * ncol(x))
+ r_eff <- .ess(.split_chains(x_centered^2)) / (nrow(x) * ncol(x))
weighted_ess <- .ess_weighted(x_centered^2, weights = weights, r_eff = r_eff)
# Kenney and Keeping (1951, eq 6.20)
@@ -681,6 +695,7 @@ mcse_sd.rvar <- function(x, ...) {
#' with other summary functions in the \pkg{posterior} package.
#'
#' @template args-conv
+#' @template args-conv-weights
#' @template args-conv-quantile
#' @param na.rm (logical) Should `NA` and `NaN` values be removed from `x` prior
#' to computing quantiles? The default is `FALSE`.
@@ -697,19 +712,19 @@ mcse_sd.rvar <- function(x, ...) {
#' quantile2(mu)
#'
#' @export
-quantile2 <- function(x, probs = c(0.05, 0.95), na.rm = FALSE, ...) {
+quantile2 <- function(x, ...) {
UseMethod("quantile2")
}
#' @rdname quantile2
#' @export
quantile2.default <- function(
- x, probs = c(0.05, 0.95), na.rm = FALSE, names = TRUE, ...
+ x, probs = c(0.05, 0.95), na.rm = FALSE, names = TRUE, weights = NULL, ...
) {
names <- as_one_logical(names)
na.rm <- as_one_logical(na.rm)
- out <- weighted_quantile(x, probs = probs, na.rm = na.rm, ...)
+ out <- weighted_quantile(x, probs = probs, na.rm = na.rm, weights = weights, ...)
if (names) {
names(out) <- paste0("q", probs * 100)
@@ -726,8 +741,8 @@ quantile2.rvar <- function(
) {
weights <- weights(x)
summarise_rvar_by_element(x, function(draws) {
- quantile2(
- draws, probs = probs, weights = weights, na.rm = na.rm, names = names, ...
+ quantile2.default(
+ draws, probs = probs, na.rm = na.rm, names = names, weights = weights, ...
)
})
}
@@ -777,6 +792,7 @@ autocovariance <- function(x) {
autocorrelation <- function(x) {
ac <- autocovariance(x)
ac <- ac / ac[1]
+ ac
}
#' Rank normalization
@@ -807,8 +823,10 @@ z_scale <- function(x, c = 3/8) {
#'
#' Compute rank uniformization for a numeric array. First replace each value by
#' its rank. Average rank for ties are used to conserve the number of unique
-#' values of discrete quantities. Second, uniformize ranks to the scale
-#' `[1/(2S), 1-1/(2S)]`, where `S` is the number of values.
+#' values of discrete quantities. Second, uniformize ranks using formula
+#' `(r - c) / (S - 2 * c + 1)`, where `r` is a rank, `S` is the number of
+#' values, and `c` is a fractional offset which defaults to c = 3/8 as
+#' recommend by Blom (1958).
#'
#' @template args-scale
#' @template args-frac-offset
diff --git a/R/draws-index.R b/R/draws-index.R
index fec535e6..6ed92c39 100644
--- a/R/draws-index.R
+++ b/R/draws-index.R
@@ -68,6 +68,10 @@ iteration_ids.draws_df <- function(x) {
#' @export
iteration_ids.draws_list <- function(x) {
+ if (!length(x) || !length(x[[1]])) {
+ # no chains or no variables within chains
+ return(integer(0))
+ }
seq_along(x[[1]][[1]])
}
diff --git a/R/for_each_draw.R b/R/for_each_draw.R
index 6aabbecb..d7cb91b7 100644
--- a/R/for_each_draw.R
+++ b/R/for_each_draw.R
@@ -7,7 +7,8 @@
#' @param expr (expression) A bare expression that can contain references to
#' variables in `x` by name. This expression will be executed once per draw
#' of `x`, where references to variables in `x` resolve to the value of that
-#' variable in that draw. The expression supports [quasiquotation].
+#' variable in that draw. The expression supports
+#' [quasiquotation][rlang::quasiquotation].
#'
#' @details
#' If `x` is not in the [`draws_rvars`] format, it is first converted to that
diff --git a/R/gpd.R b/R/gpd.R
index 924f8c3f..3cffd77c 100644
--- a/R/gpd.R
+++ b/R/gpd.R
@@ -21,22 +21,28 @@ qgeneralized_pareto <- function(p, mu = 0, sigma = 1, k = 0, lower.tail = TRUE,
#' 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} (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 (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
-#' 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.
+#' 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} (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 (see details in Vehtari et al., 2024). This is used
+#' internally but is exported for use by other packages.
+#' @family helper-functions
+#' @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 &
@@ -67,7 +73,8 @@ 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 (see Vehtari et al., 2022 for details)
+ # this stabilizes estimates for very small Monte Carlo sample sizes and low ess
+ # (see Vehtari et al., 2024 for details)
if (wip) {
k_hat <- (k_hat * N + 0.5 * 10) / (N + 10)
}
diff --git a/R/misc.R b/R/misc.R
index c1bce76c..da26fd0f 100644
--- a/R/misc.R
+++ b/R/misc.R
@@ -1,4 +1,4 @@
-# initialize a named list
+ # initialize a named list
# @param names names of the elements
# @param values optional values of the elements
named_list <- function(names, values = NULL) {
@@ -209,9 +209,16 @@ escape_all <- function(x) {
# numerically stable version of log(sum(exp(x)))
log_sum_exp <- function(x) {
- max <- max(x)
- sum <- sum(exp(x - max))
- max + log(sum)
+ max <- max(as.numeric(x), warnings = FALSE)
+ if (max == -Inf) {
+ res <- 0
+ } else if (max == Inf) {
+ res <- Inf
+ } else {
+ sum <- sum(exp(x - max))
+ res <- max + log(sum)
+ }
+ res
}
# simple version of destructuring assignment
@@ -238,3 +245,4 @@ has_s3_method <- function(f, signature) {
}
FALSE
}
+
diff --git a/R/nested_rhat.R b/R/nested_rhat.R
index f1b81901..c65f0c95 100644
--- a/R/nested_rhat.R
+++ b/R/nested_rhat.R
@@ -1,7 +1,7 @@
#' Nested Rhat convergence diagnostic
#'
#' Compute the nested Rhat convergence diagnostic for a single
-#' variable as proposed in Margossian et al. (2023).
+#' variable as proposed in Margossian et al. (2024).
#'
#' @family diagnostics
#' @template args-conv
@@ -20,11 +20,11 @@
#' Note that there is a slight difference in the calculation of Rhat
#' and nested Rhat, as nested Rhat is lower bounded by 1. This means
#' that nested Rhat with one chain per superchain will not be
-#' exactly equal to basic Rhat (see Footnote 1 in Margossian et
-#' al. (2023)).
+#' exactly equal to basic Rhat (see Footnote 3 in Margossian et
+#' al. (2024)).
#'
#' @template return-conv
-#' @template ref-margossian-nestedrhat-2023
+#' @template ref-margossian-nestedrhat-2024
#'
#' @examples
#' mu <- extract_variable_matrix(example_draws(), "mu")
@@ -78,17 +78,17 @@ rhat_nested.rvar <- function(x, superchain_ids, ...) {
superchains <- unique(superchain_ids)
- # mean and variance of chains calculated as in rhat
+ # mean and variance of chains calculated as in Rhat
chain_mean <- matrixStats::colMeans2(x)
chain_var <- matrixStats::colVars(x, center = chain_mean)
# mean of superchains calculated by only including specified chains
- # (equation 15 in Margossian et al. 2023)
+ # (equation 4 in Margossian et al. 2024)
superchain_mean <- sapply(
superchains, function(k) mean(x[, which(superchain_ids == k)])
)
- # between-chain variance estimate (B_k in equation 18 in Margossian et al. 2023)
+ # between-chain variance estimate (Bhat_k in equation 7 in Margossian et al. 2024)
if (nchains_per_superchain == 1) {
var_between_chain <- 0
} else {
@@ -98,7 +98,7 @@ rhat_nested.rvar <- function(x, superchain_ids, ...) {
)
}
- # within-chain variance estimate (W_k in equation 18 in Margossian et al. 2023)
+ # within-chain variance estimate (What_k in equation 7 in Margossian et al. 2024)
if (niterations == 1) {
var_within_chain <- 0
} else {
@@ -108,12 +108,12 @@ rhat_nested.rvar <- function(x, superchain_ids, ...) {
)
}
- # between-superchain variance (nB in equation 17 in Margossian et al. 2023)
+ # between-superchain variance (Bhat_nu in equation 6 in Margossian et al. 2024)
var_between_superchain <- var(superchain_mean)
- # within-superchain variance (nW in equation 18 in Margossian et al. 2023)
+ # within-superchain variance (What_nu in equation 7 in Margossian et al. 2024)
var_within_superchain <- mean(var_within_chain + var_between_chain)
- # nested Rhat (nRhat in equation 19 in Margossian et al. 2023)
+ # nested Rhat (Rhat_nu in equation 8 in Margossian et al. 2024)
sqrt(1 + var_between_superchain / var_within_superchain)
}
diff --git a/R/pareto_smooth.R b/R/pareto_smooth.R
index 794c7188..ed01c4ed 100644
--- a/R/pareto_smooth.R
+++ b/R/pareto_smooth.R
@@ -3,13 +3,13 @@
#' Estimate Pareto k value by fitting a Generalized Pareto
#' 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).
+#' diagnostics. For further details see Vehtari et al. (2024).
#'
#' @family diagnostics
#' @template args-pareto
#' @template args-methods-dots
#' @template ref-vehtari-paretosmooth-2022
-#' @return `khat` estimated Generalized Pareto Distribution shape parameter k
+#' @template return-conv
#'
#' @seealso [`pareto_diags`] for additional related diagnostics, and
#' [`pareto_smooth`] for Pareto smoothed draws.
@@ -39,35 +39,35 @@ pareto_khat.default <- function(x,
verbose = verbose,
return_k = TRUE,
smooth_draws = FALSE,
- are_log_weights = are_log_weights,
- ...)
- return(smoothed$diagnostics)
+ are_log_weights = are_log_weights
+ )
+ return(smoothed$diagnostics$khat)
}
#' @rdname pareto_khat
#' @export
-pareto_khat.rvar <- function(x, verbose = FALSE, ...) {
+pareto_khat.rvar <- function(x, ...) {
if (is.null(weights(x))) {
draws_diags <- summarise_rvar_by_element_with_chains(
x,
pareto_smooth.default,
smooth_draws = FALSE,
return_k = TRUE,
- verbose = verbose,
...
)
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)
- )
+ diags <- apply(draws_diags, margins, function(x) x[[1]]$diagnostics$khat)
+
} else {
# take the max of khat for x * weights and khat for weights
+ lw <- weights(x, log = TRUE)
+
weights_diags <- pareto_khat(
- weights(x, log = TRUE),
+ lw,
are_log_weights = TRUE,
...
)
@@ -80,20 +80,17 @@ pareto_khat.rvar <- function(x, verbose = FALSE, ...) {
product_diags <- summarise_rvar_by_element_with_chains(
xu,
pareto_khat.default,
- verbose = verbose,
...
)
dim(product_diags) <- dim(product_diags) %||% length(product_diags)
margins <- seq_along(dim(product_diags))
- diags <- list(
- khat = apply(product_diags, margins,
+
+ diags <- apply(product_diags, margins,
function(x) {
- max(x[[1]]$khat,
- weights_diags$khat)
+ max(x, weights_diags)
})
- )
}
diags
}
@@ -125,24 +122,26 @@ pareto_khat.rvar <- function(x, verbose = FALSE, ...) {
#' 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).
+#' Section 3.2.3, Equation 11 in Vehtari et al. (2024).
#'
#' * `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).
+#' in Vehtari et al. (2024).
#'
#' * `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).
+#' theorem convergence rate. See Appendix B in Vehtari et al. (2024).
#'
-#' @seealso [`pareto_khat`] for only calculating khat, and
-#' [`pareto_smooth`] for Pareto smoothed draws.
+#' @seealso [`pareto_khat`], [`pareto_min_ss`],
+#' [`pareto_khat_threshold`], and [`pareto_convergence_rate`] for
+#' individual diagnostics; and [`pareto_smooth`] for Pareto smoothing
+#' draws.
#' @examples
#' mu <- extract_variable_matrix(example_draws(), "mu")
#' pareto_diags(mu)
@@ -187,10 +186,7 @@ pareto_diags.rvar <- function(x, ...) {
if (is.null(weights(x))) {
draws_diags <- summarise_rvar_by_element_with_chains(
x,
- pareto_smooth.default,
- return_k = TRUE,
- smooth_draws = FALSE,
- extra_diags = TRUE,
+ pareto_diags.default,
...
)
@@ -198,10 +194,10 @@ pareto_diags.rvar <- function(x, ...) {
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)
+ khat = apply(draws_diags, margins, function(x) x[[1]]$khat),
+ min_ss = apply(draws_diags, margins, function(x) x[[1]]$min_ss),
+ khat_threshold = apply(draws_diags, margins, function(x) x[[1]]$khat_threshold),
+ convergence_rate = apply(draws_diags, margins, function(x) x[[1]]$convergence_rate)
)
} else {
@@ -242,7 +238,7 @@ pareto_diags.rvar <- function(x, ...) {
#'
#' 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).
+#' tail(s). For further details see Vehtari et al. (2024).
#'
#' @template args-pareto
#' @param return_k (logical) Should the Pareto khat be included in
@@ -257,13 +253,20 @@ pareto_diags.rvar <- function(x, ...) {
#' @template ref-vehtari-paretosmooth-2022
#' @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
-#' estimates * `convergence_rate`: Relative convergence rate for
+#' containing numeric values:
+#'
+#' * `khat`: estimated Pareto k shape parameter, and optionally
+#' * `min_ss`: minimum sample size for reliable Pareto smoothed
+#' estimate
+#' * `khat_threshold`: sample size specific khat threshold for
+#' reliable Pareto smoothed estimates
+#' * `convergence_rate`: Relative convergence rate for
#' Pareto smoothed estimates
#'
+#' If any of the draws is non-finite, that is, `NA`, `NaN`, `Inf`, or
+#' `-Inf`, Pareto smoothing will not be performed, and the original
+#' draws will be returned and and diagnostics will be `NA` (numeric).
+#'
#' @seealso [`pareto_khat`] for only calculating khat, and
#' [`pareto_diags`] for additional diagnostics.
#' @examples
@@ -323,22 +326,25 @@ pareto_smooth.default <- function(x,
are_log_weights = FALSE,
...) {
- checkmate::assert_numeric(ndraws_tail, null.ok = TRUE)
- checkmate::assert_numeric(r_eff, null.ok = TRUE)
+ if (!is.null(r_eff)) {
+ r_eff <- as_one_numeric(r_eff)
+ }
+ if (!is.null(ndraws_tail)) {
+ ndraws_tail <- as_one_integer(ndraws_tail)
+ }
extra_diags <- as_one_logical(extra_diags)
return_k <- as_one_logical(return_k)
verbose <- as_one_logical(verbose)
are_log_weights <- as_one_logical(are_log_weights)
+ if (extra_diags) {
+ return_k <- TRUE
+ }
+
# check for infinite or na values
if (should_return_NA(x)) {
- warning_no_call("Input contains infinite or NA values, or is constant. Fitting of generalized Pareto distribution not performed.")
- if (!return_k) {
- out <- x
- } else {
- out <- list(x = x, diagnostics = NA_real_)
- }
- return(out)
+ warning_no_call("Input contains infinite or NA values, is constant or has constant tail. Fitting of generalized Pareto distribution not performed.")
+ return(pareto_diags_na(x, return_k, extra_diags))
}
if (are_log_weights) {
@@ -346,25 +352,30 @@ pareto_smooth.default <- function(x,
}
tail <- match.arg(tail)
- S <- length(x)
+ ndraws <- length(x)
# automatically calculate relative efficiency
if (is.null(r_eff)) {
- r_eff <- ess_tail(x) / S
+ r_eff <- ess_tail(x) / ndraws
}
# automatically calculate tail length
if (is.null(ndraws_tail)) {
- ndraws_tail <- ps_tail_length(S, r_eff)
+ ndraws_tail <- ps_tail_length(ndraws, r_eff)
+ }
+
+ if (is.na(ndraws_tail)) {
+ warning_no_call("Input contains infinite or NA values, is constant, or has constant tail. Fitting of generalized Pareto distribution not performed.")
+ return(pareto_diags_na(x, return_k, extra_diags))
}
if (tail == "both") {
- if (ndraws_tail > S / 2) {
+ if (ndraws_tail > ndraws / 2) {
warning("Number of tail draws cannot be more than half ",
"the total number of draws if both tails are fit, ",
- "changing to ", S / 2, ".")
- ndraws_tail <- S / 2
+ "changing to ", ndraws / 2, ".")
+ ndraws_tail <- ndraws / 2
}
if (ndraws_tail < 5) {
@@ -393,7 +404,7 @@ pareto_smooth.default <- function(x,
)
right_k <- smoothed$k
- k <- max(left_k, right_k)
+ k <- max(left_k, right_k, na.rm = TRUE)
x <- smoothed$x
} else {
@@ -412,7 +423,7 @@ pareto_smooth.default <- function(x,
diags_list <- list(khat = k)
if (extra_diags) {
- ext_diags <- .pareto_smooth_extra_diags(k, S)
+ ext_diags <- .pareto_smooth_extra_diags(k, ndraws)
diags_list <- c(diags_list, ext_diags)
}
@@ -444,13 +455,13 @@ pareto_khat_threshold <- function(x, ...) {
#' @rdname pareto_diags
#' @export
pareto_khat_threshold.default <- function(x, ...) {
- c(khat_threshold = ps_khat_threshold(length(x)))
+ ps_khat_threshold(length(x))
}
#' @rdname pareto_diags
#' @export
pareto_khat_threshold.rvar <- function(x, ...) {
- c(khat_threshold = ps_khat_threshold(ndraws(x)))
+ ps_khat_threshold(ndraws(x))
}
#' @rdname pareto_diags
@@ -462,15 +473,15 @@ pareto_min_ss <- function(x, ...) {
#' @rdname pareto_diags
#' @export
pareto_min_ss.default <- function(x, ...) {
- k <- pareto_khat(x)$k
- c(min_ss = ps_min_ss(k))
+ k <- pareto_khat(x)
+ ps_min_ss(k)
}
#' @rdname pareto_diags
#' @export
pareto_min_ss.rvar <- function(x, ...) {
- k <- pareto_khat(x)$k
- c(min_ss = ps_min_ss(k))
+ k <- pareto_khat(x)
+ ps_min_ss(k)
}
#' @rdname pareto_diags
@@ -482,15 +493,15 @@ pareto_convergence_rate <- function(x, ...) {
#' @rdname pareto_diags
#' @export
pareto_convergence_rate.default <- function(x, ...) {
- k <- pareto_khat(x)$khat
- c(convergence_rate = ps_convergence_rate(k, length(x)))
+ k <- pareto_khat(x)
+ ps_convergence_rate(k, length(x))
}
#' @rdname pareto_diags
#' @export
pareto_convergence_rate.rvar <- function(x, ...) {
k <- pareto_khat(x)
- c(convergence_rate = ps_convergence_rate(k, ndraws(x)))
+ ps_convergence_rate(k, ndraws(x))
}
@@ -512,8 +523,8 @@ pareto_convergence_rate.rvar <- function(x, ...) {
tail <- match.arg(tail)
- S <- length(x)
- tail_ids <- seq(S - ndraws_tail + 1, S)
+ ndraws <- length(x)
+ tail_ids <- seq(ndraws - ndraws_tail + 1, ndraws)
if (tail == "left") {
x <- -x
@@ -522,6 +533,16 @@ pareto_convergence_rate.rvar <- function(x, ...) {
ord <- sort.int(x, index.return = TRUE)
draws_tail <- ord$x[tail_ids]
+ if (is_constant(draws_tail)) {
+
+ if (tail == "left") {
+ x <- -x
+ }
+
+ out <- list(x = x, k = NA)
+ return(out)
+ }
+
cutoff <- ord$x[min(tail_ids) - 1] # largest value smaller than tail values
max_tail <- max(draws_tail)
@@ -582,15 +603,15 @@ pareto_convergence_rate.rvar <- function(x, ...) {
#' Extra pareto-k diagnostics
#'
#' internal function to calculate the extra diagnostics for a given
-#' pareto k and sample size S
+#' pareto k and number of draws ndraws
#' @noRd
-.pareto_smooth_extra_diags <- function(k, S, ...) {
+.pareto_smooth_extra_diags <- function(k, ndraws, ...) {
min_ss <- ps_min_ss(k)
- khat_threshold <- ps_khat_threshold(S)
+ khat_threshold <- ps_khat_threshold(ndraws)
- convergence_rate <- ps_convergence_rate(k, S)
+ convergence_rate <- ps_convergence_rate(k, ndraws)
other_diags <- list(
min_ss = min_ss,
@@ -602,12 +623,15 @@ pareto_convergence_rate.rvar <- function(x, ...) {
#' 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)
-#' Equation (11) in PSIS paper
+#' smoothed estimate (to have small probability of large error). See
+#' section 3.2.3 in Vehtari et al. (2024). This function is exported
+#' to be usable by other packages. For user-facing diagnostic functions, see
+#' [`pareto_min_ss`] and [`pareto_diags`].
+#' @family helper-functions
#' @param k pareto k value
#' @param ... unused
#' @return minimum sample size
-#' @noRd
+#' @export
ps_min_ss <- function(k, ...) {
if (k < 1) {
out <- 10^(1 / (1 - max(0, k)))
@@ -617,29 +641,37 @@ ps_min_ss <- function(k, ...) {
out
}
-#' Pareto-smoothing k-hat threshold
+#' Pareto k-hat threshold
#'
-#' Given sample size S computes khat threshold for reliable Pareto
+#' Given number of draws, 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
+#' section 3.2.4, equation (13) of Vehtari et al. (2024). This
+#' function is exported to be usable by other packages. For
+#' user-facing diagnostic functions, see [`pareto_khat_threshold`] and
+#' [`pareto_diags`].
+#' @family helper-functions
+#' @param ndraws number of draws
#' @param ... unused
#' @return threshold
-#' @noRd
-ps_khat_threshold <- function(S, ...) {
- 1 - 1 / log10(S)
+#' @export
+ps_khat_threshold <- function(ndraws, ...) {
+ 1 - 1 / log10(ndraws)
}
-#' Pareto-smoothing convergence rate
+#' Pareto convergence rate
#'
-#' Given S and scalar or array of k's, compute the relative
-#' convergence rate of PSIS estimate RMSE
+#' Given number of draws and scalar or array of k's, compute the
+#' relative convergence rate of PSIS estimate RMSE. See Appendix B of
+#' Vehtari et al. (2024). This function is exported to be usable by
+#' other packages. For user-facing diagnostic functions, see
+#' [`pareto_convergence_rate`] and [`pareto_diags`].
+#' @family helper-functions
#' @param k pareto-k values
-#' @param S sample size
+#' @param ndraws number of draws
#' @param ... unused
#' @return convergence rate
-#' @noRd
-ps_convergence_rate <- function(k, S, ...) {
+#' @export
+ps_convergence_rate <- function(k, ndraws, ...) {
# allow array of k's
rate <- numeric(length(k))
# k<0 bounded distribution
@@ -647,28 +679,41 @@ ps_convergence_rate <- function(k, S, ...) {
# 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)
+ rate[k == 0.5] <- 1 - 1 / log(ndraws)
+ # smooth approximation for the rest (see Appendix B 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)))
+ (2 * (kk - 1) * ndraws^(2 * kk + 1) + (1 - 2 * kk) * ndraws^(2 * kk) + ndraws^2) /
+ ((ndraws - 1) * (ndraws - ndraws^(2 * kk)))
)
rate
}
-#' Calculate the tail length from S and r_eff
-#' Appendix H in PSIS paper
-#' @noRd
-ps_tail_length <- function(S, r_eff, ...) {
- ifelse(S > 225, ceiling(3 * sqrt(S / r_eff)), S / 5)
+#' Pareto tail length
+#'
+#' Calculate the tail length from number of draws and relative efficiency
+#' r_eff. See Appendix H in Vehtari et al. (2024). This function is
+#' used internally and is exported to be available for other packages.
+#' @family helper-functions
+#' @param ndraws number of draws
+#' @param r_eff relative efficiency
+#' @param ... unused
+#' @return tail length
+#' @export
+ps_tail_length <- function(ndraws, r_eff, ...) {
+ floor(
+ ifelse(
+ ndraws * r_eff > 225,
+ 3 * sqrt(ndraws / r_eff),
+ ndraws / 5)
+ )
}
#' Pareto-k diagnostic message
#'
-#' Given S and scalar and k, form a diagnostic message string
+#' Given number of draws and k, form a diagnostic message string.
#' @param diags (numeric) named vector of diagnostic values
#' @param are_weights (logical) are the diagnostics for weights
#' @param ... unused
@@ -701,3 +746,23 @@ pareto_k_diagmsg <- function(diags, are_weights = FALSE, ...) {
message("Pareto k-hat = ", round(khat, 2), ".", msg)
invisible(diags)
}
+
+
+pareto_diags_na <- function(x, return_k, extra_diags) {
+ if (!return_k) {
+ out <- x
+ } else if (!extra_diags) {
+ out <- list(x = x, diagnostics = list(khat = NA_real_))
+ } else {
+ out <- list(
+ x = x,
+ diagnostics = list(
+ khat = NA_real_,
+ min_ss = NA_real_,
+ khat_threshold = NA_real_,
+ convergence_rate = NA_real_
+ )
+ )
+ }
+ return(out)
+}
diff --git a/R/pit.R b/R/pit.R
new file mode 100644
index 00000000..4f4006f9
--- /dev/null
+++ b/R/pit.R
@@ -0,0 +1,179 @@
+#' Probability integral transform
+#'
+#' Probability integral transform (PIT). LOO-PIT is given by a weighted sample.
+#'
+#' @name pit
+#'
+#' @param x (draws) A [`draws_matrix`] object or one coercible to a
+#' `draws_matrix` object, or an [`rvar`] object.
+#'
+#' @param y (observations) A 1D vector, or an array of dim(x), if x is `rvar`.
+#' Each element of `y` corresponds to a variable in `x`.
+#'
+#' @param weights A matrix of weights for each draw and variable. `weights`
+#' should have one column per variable in `x`, and `ndraws(x)` rows.
+#'
+#' @param log (logical) Are the weights passed already on the log scale? The
+#' default is `FALSE`, that is, expecting `weights` to be on the standard
+#' (non-log) scale.
+#'
+#' @template args-methods-dots
+#'
+#' @details The `pit()` function computes the probability integral transform of
+#' `y` using the empirical cumulative distribution computed from the samples
+#' in `x`. For continuous valued `y` and `x`, the PIT for the elements of `y`
+#' is computed as the empirical cumulative distribution value:
+#'
+#' PIT(y_i) = Pr(x_i < y_i),
+#'
+#' where x_i, is the corresponding set of draws in `x`. For `draws` objects,
+#' this corresponds to the draws of the *i*th variable, and for `rvar`
+#' the elements of `y` and `x` are matched.
+#'
+#' The draws in `x` can further be provided (log-)weights in
+# `weights`, which enables for example the computation of LOO-PITs.
+#'
+#' If `y` and `x` are discrete, randomisation is used to obtain continuous PIT
+#' values. (see, e.g., Czado, C., Gneiting, T., Held, L.: Predictive model
+#' assessment for count data. Biometrics 65(4), 1254–1261 (2009).)
+#'
+#' @return A numeric vector of length `length(y)` containing the PIT values, or
+#' an array of shape `dim(y)`, if `x` is an `rvar`.
+
+#' @examples
+#' # PIT for a draws object
+#' x <- example_draws()
+#' # Create a vector of observations
+#' y <- rnorm(nvariables(x), 5, 5)
+#' pit(x, y)
+#'
+#' # Compute weighted PIT (for example LOO-PIT)
+#' weights <- matrix(runif(length(x)), ncol = nvariables(x))
+#'
+#' pit(x, y, weights)
+#'
+#' # PIT for an rvar
+#' x <- rvar(example_draws())
+#' # Create an array of observations with the same dimensions as x.
+#' y_arr <- array(rnorm(length(x), 5, 5), dim = dim(x))
+#' pit(x, y_arr)
+#'
+NULL
+
+#' @rdname pit
+#' @export
+pit <- function(x, y, ...) UseMethod("pit")
+
+#' @rdname pit
+#' @export
+pit.default <- function(x, y, weights = NULL, log = FALSE, ...) {
+ x <- as_draws_matrix(x)
+ if (!is.null(weights)) {
+ weights <- as_draws_matrix(weights)
+ }
+ pit(x, y, weights, log)
+}
+
+#' @rdname pit
+#' @export
+pit.draws_matrix <- function(x, y, weights = NULL, log = FALSE, ...) {
+ y <- validate_y(y, x)
+ if (!is.null(weights)) {
+ weights <- sapply(seq_len(nvariables(x)), \(var_idx) {
+ validate_weights(weights[, var_idx], ndraws(x), log)
+ })
+ weights <- normalize_log_weights(weights)
+ }
+ pit <- vapply(seq_len(ncol(x)), function(j) {
+ sel_min <- x[, j] < y[j]
+ if (!any(sel_min)) {
+ pit <- 0
+ } else {
+ if (is.null(weights)) {
+ pit <- mean(sel_min)
+ } else {
+ pit <- exp(log_sum_exp(weights[sel_min, j]))
+ }
+ }
+
+ sel_sup <- x[, j] == y[j]
+ if (any(sel_sup)) {
+ # randomized PIT for discrete y (see, e.g., Czado, C., Gneiting, T.,
+ # Held, L.: Predictive model assessment for count data.
+ # Biometrics 65(4), 1254–1261 (2009).)
+ if (is.null(weights)) {
+ pit_sup <- pit + mean(sel_sup)
+ } else {
+ pit_sup <- pit + exp(log_sum_exp(weights[sel_sup, j]))
+ }
+
+ pit <- runif(1, pit, pit_sup)
+ }
+ pit
+ }, FUN.VALUE = 1.0)
+
+ if (any(pit > 1 + 1e-10)) {
+ warning_no_call(
+ paste(
+ "Some PIT values larger than 1. ",
+ "This is usually due to numerical inaccuracies. ",
+ "Largest value: ",
+ max(pit),
+ "\nRounding PIT > 1 to 1.",
+ sep = ""
+ )
+ )
+ }
+
+ setNames(pmin(1, pit), variables(x))
+}
+
+#' @rdname pit
+#' @export
+pit.rvar <- function(x, y, weights = NULL, log = FALSE, ...) {
+ y <- validate_y(y, x)
+ if (is.null(weights)) {
+ out <- array(
+ runif(length(y), Pr(x < y), Pr(x <= y)),
+ dim(x),
+ dimnames(x)
+ )
+ } else {
+ out <- array(
+ data = pit(
+ x = as_draws_matrix(c(x)),
+ y = c(y),
+ weights = weights,
+ log = log
+ ),
+ dim = dim(x),
+ dimnames = dimnames(x)
+ )
+ }
+ out
+}
+
+# internal ----------------------------------------------------------------
+
+validate_y <- function(y, x = NULL) {
+ if (!is.numeric(y)) {
+ stop_no_call("`y` must be numeric.")
+ }
+ if (anyNA(y)) {
+ stop_no_call("NAs not allowed in `y`.")
+ }
+ if (is_rvar(x)) {
+ if (length(x) != length(y) || any(dim(y) != dim(x))) {
+ stop_no_call("`dim(y)` must match `dim(x)`.")
+ }
+ } else if (is_draws(x)) {
+ if (!is.vector(y, mode = "numeric") || length(y) != nvariables(x)) {
+ stop_no_call("`y` must be a vector of length `nvariables(x)`.")
+ }
+ }
+ y
+}
+
+normalize_log_weights <- function(log_weights) {
+ apply(log_weights, 2, function(col) col - log_sum_exp(col))
+}
diff --git a/R/rvar-.R b/R/rvar-.R
index 32ecdd05..3f2ba2e4 100755
--- a/R/rvar-.R
+++ b/R/rvar-.R
@@ -920,6 +920,13 @@ drop_chain_dim <- function(x) {
#' (first one is draws), etc.
#' @noRd
cleanup_rvar_draws <- function(x) {
+ assert(
+ check_numeric(x),
+ check_logical(x),
+ check_factor(x),
+ check_character(x)
+ )
+
if (length(x) == 0) {
# canonical NULL rvar is at least 1 draw of nothing
# this ensures that (e.g.) extending a null rvar
diff --git a/R/rvar-rfun.R b/R/rvar-rfun.R
index 79210e6a..47bdeb0b 100755
--- a/R/rvar-rfun.R
+++ b/R/rvar-rfun.R
@@ -116,7 +116,7 @@ rfun <- function (.f, rvar_args = NULL, rvar_dots = TRUE, ndraws = NULL) {
#' producing a new [`rvar`].
#'
#' @param expr (expression) A bare expression that can (optionally) contain
-#' [`rvar`]s. The expression supports [quasiquotation].
+#' [`rvar`]s. The expression supports [quasiquotation][rlang::quasiquotation].
#' @param ndraws (positive integer) The number of draws used to construct new
#' random variables if no [`rvar`]s are supplied in `expr`. If `NULL`,
#' `getOption("posterior.rvar_ndraws")` is used (default 4000). If `expr`
diff --git a/R/subset_draws.R b/R/subset_draws.R
index a765f435..b2b7b9a9 100644
--- a/R/subset_draws.R
+++ b/R/subset_draws.R
@@ -22,6 +22,9 @@
#' If `FALSE` (the default) only the selected subset will be
#' returned. If `TRUE` everything but the selected subset will be
#' returned.
+#' @param scalar (logical) Should only scalar variables be selected?
+#' If `FALSE` (the default), all variables with matching names and
+#' *arbitrary* indices will be selected (see examples).
#'
#' @template args-methods-dots
#' @template return-draws
@@ -43,6 +46,9 @@
#' # extract all elements of 'theta'
#' subset_draws(x, variable = "theta")
#'
+#' # trying to extract only a scalar 'theta' will fail
+#' # subset_draws(x, variable = "theta", scalar = TRUE)
+#'
#' @export
subset_draws <- function(x, ...) {
UseMethod("subset_draws")
@@ -52,12 +58,16 @@ subset_draws <- function(x, ...) {
#' @export
subset_draws.draws_matrix <- function(x, variable = NULL, iteration = NULL,
chain = NULL, draw = NULL, regex = FALSE,
- unique = TRUE, exclude = FALSE, ...) {
+ unique = TRUE, exclude = FALSE,
+ scalar = FALSE, ...) {
if (all_null(variable, iteration, chain, draw)) {
return(x)
}
x <- repair_draws(x)
- variable <- check_existing_variables(variable, x, regex = regex, exclude = exclude)
+ variable <- check_existing_variables(
+ variable, x, regex = regex, exclude = exclude,
+ scalar = scalar
+ )
iteration <- check_iteration_ids(iteration, x, unique = unique, exclude = exclude)
chain <- check_chain_ids(chain, x, unique = unique, exclude = exclude)
draw <- check_draw_ids(draw, x, unique = unique, exclude = exclude)
@@ -74,13 +84,17 @@ subset_draws.draws_matrix <- function(x, variable = NULL, iteration = NULL,
#' @export
subset_draws.draws_array <- function(x, variable = NULL, iteration = NULL,
chain = NULL, draw = NULL, regex = FALSE,
- unique = TRUE, exclude = FALSE, ...) {
+ unique = TRUE, exclude = FALSE,
+ scalar = FALSE, ...) {
if (all_null(variable, iteration, chain, draw)) {
return(x)
}
x <- repair_draws(x)
- variable <- check_existing_variables(variable, x, regex = regex, exclude = exclude)
+ variable <- check_existing_variables(
+ variable, x, regex = regex, exclude = exclude,
+ scalar = scalar
+ )
iteration <- check_iteration_ids(iteration, x, unique = unique, exclude = exclude)
chain <- check_chain_ids(chain, x, unique = unique, exclude = exclude)
draw <- check_draw_ids(draw, x, unique = unique, exclude = exclude)
@@ -100,17 +114,21 @@ subset_draws.draws_array <- function(x, variable = NULL, iteration = NULL,
#' @export
subset_draws.draws_df <- function(x, variable = NULL, iteration = NULL,
chain = NULL, draw = NULL, regex = FALSE,
- unique = TRUE, exclude = FALSE, ...) {
+ unique = TRUE, exclude = FALSE,
+ scalar = FALSE, ...) {
if (all_null(variable, iteration, chain, draw)) {
return(x)
}
x <- repair_draws(x)
unique <- as_one_logical(unique)
- variable <- check_existing_variables(variable, x, regex = regex, exclude= exclude)
- iteration <- check_iteration_ids(iteration, x, unique = unique, exclude= exclude)
- chain <- check_chain_ids(chain, x, unique = unique, exclude= exclude)
- draw <- check_draw_ids(draw, x, unique = unique, exclude= exclude)
+ variable <- check_existing_variables(
+ variable, x, regex = regex, exclude = exclude,
+ scalar = scalar
+ )
+ iteration <- check_iteration_ids(iteration, x, unique = unique, exclude = exclude)
+ chain <- check_chain_ids(chain, x, unique = unique, exclude = exclude)
+ draw <- check_draw_ids(draw, x, unique = unique, exclude = exclude)
x <- prepare_subsetting(x, iteration, chain, draw)
x <- .subset_draws(
@@ -124,13 +142,17 @@ subset_draws.draws_df <- function(x, variable = NULL, iteration = NULL,
#' @export
subset_draws.draws_list <- function(x, variable = NULL, iteration = NULL,
chain = NULL, draw = NULL, regex = FALSE,
- unique = TRUE, exclude = FALSE, ...) {
+ unique = TRUE, exclude = FALSE,
+ scalar = FALSE, ...) {
if (all_null(variable, iteration, chain, draw)) {
return(x)
}
x <- repair_draws(x)
- variable <- check_existing_variables(variable, x, regex = regex, exclude = exclude)
+ variable <- check_existing_variables(
+ variable, x, regex = regex, exclude = exclude,
+ scalar = scalar
+ )
iteration <- check_iteration_ids(iteration, x, unique = unique, exclude = exclude)
chain <- check_chain_ids(chain, x, unique = unique, exclude = exclude)
draw <- check_draw_ids(draw, x, unique = unique, exclude = exclude)
@@ -150,16 +172,20 @@ subset_draws.draws_list <- function(x, variable = NULL, iteration = NULL,
#' @export
subset_draws.draws_rvars <- function(x, variable = NULL, iteration = NULL,
chain = NULL, draw = NULL, regex = FALSE,
- unique = TRUE, exclude = FALSE, ...) {
+ unique = TRUE, exclude = FALSE,
+ scalar = FALSE, ...) {
if (all_null(variable, iteration, chain, draw)) {
return(x)
}
x <- repair_draws(x)
- variable <- check_existing_variables(variable, x, regex = regex, exclude = exclude)
- iteration <- check_iteration_ids(iteration, x, unique = unique, exclude= exclude)
- chain <- check_chain_ids(chain, x, unique = unique, exclude= exclude)
- draw <- check_draw_ids(draw, x, unique = unique, exclude= exclude)
+ variable <- check_existing_variables(
+ variable, x, regex = regex, exclude = exclude,
+ scalar = scalar
+ )
+ iteration <- check_iteration_ids(iteration, x, unique = unique, exclude = exclude)
+ chain <- check_chain_ids(chain, x, unique = unique, exclude = exclude)
+ draw <- check_draw_ids(draw, x, unique = unique, exclude = exclude)
x <- prepare_subsetting(x, iteration, chain, draw)
if (!is.null(draw)) {
diff --git a/R/summarise_draws.R b/R/summarise_draws.R
index 05acd83b..6f13755a 100644
--- a/R/summarise_draws.R
+++ b/R/summarise_draws.R
@@ -329,9 +329,9 @@ empty_draws_summary <- function(dimensions = "variable") {
create_summary_list <- function(x, v, funs, .args) {
draws <- drop_dims_or_classes(x[, , v], dims = 3, reset_class = FALSE)
+ args <- c(list(draws), .args)
v_summary <- named_list(names(funs))
for (m in names(funs)) {
- args <- c(list(draws), .args[[m]])
v_summary[[m]] <- do.call(funs[[m]], args)
}
v_summary
diff --git a/R/variables.R b/R/variables.R
index 49911d15..e0c732a0 100644
--- a/R/variables.R
+++ b/R/variables.R
@@ -224,19 +224,22 @@ nvariables.draws <- function(x, ...) {
# check validity of existing variable names: e.g., that
# all `variables` exist in `x` and that no `variables`are reserved words
-# Additionally, this returns the cannonical name, so e.g. "theta" will get
+# Additionally, this returns the canonical name, so e.g. "theta" will get
# converted to c("theta[1]", "theta[2]", ...) if those variables exist.
# @param regex should 'variables' be treated as regular expressions?
-# @param scalar_only should only scalar variables be matched?
+# @param scalar should only scalar variables be matched?
check_existing_variables <- function(variables, x, regex = FALSE,
- scalar_only = FALSE, exclude = FALSE) {
+ scalar = FALSE, exclude = FALSE) {
check_draws_object(x)
if (is.null(variables)) {
return(NULL)
}
+ # do not allow NA values in variables
+ checkmate::assertCharacter(variables, any.missing = FALSE)
+
regex <- as_one_logical(regex)
- scalar_only <- as_one_logical(scalar_only)
+ scalar <- as_one_logical(scalar)
exclude <- as_one_logical(exclude)
variables <- unique(as.character(variables))
all_variables <- variables(x, reserved = TRUE)
@@ -249,7 +252,7 @@ check_existing_variables <- function(variables, x, regex = FALSE,
# regular expressions are not required to match anything
missing_variables <- NULL
variables <- as.character(all_variables[unique(unlist(tmp))])
- } else if (!scalar_only) {
+ } else if (!scalar) {
# need to find variables that are matched by either a scalar or vector
# variable in x and what the matching variable is, while keeping original
# order of input `variables`
diff --git a/R/weight_draws.R b/R/weight_draws.R
index a336e177..2469dd06 100644
--- a/R/weight_draws.R
+++ b/R/weight_draws.R
@@ -213,14 +213,17 @@ log_weights.rvar <- function(object, ...) {
# validate weights and return log weights
validate_weights <- function(weights, ndraws, log = FALSE, pareto_smooth = FALSE) {
if (is.null(weights)) return(NULL)
- assert_numeric(weights)
- assert_atomic_vector(weights)
- assert_flag(log)
- assert_flag(pareto_smooth)
+ checkmate::assert_numeric(weights)
+# checkmate::assert_atomic_vector(weights)
+ checkmate::assert_flag(log)
+ checkmate::assert_flag(pareto_smooth)
if (length(weights) != ndraws) {
stop_no_call("Number of weights must match the number of draws.")
}
+ if (any(weights == Inf)) {
+ stop_no_call("Weights must not be positive infinite.")
+ }
if (!log) {
if (any(weights < 0)) {
stop_no_call("Weights must be non-negative.")
diff --git a/R/weighted.R b/R/weighted.R
index f444ef74..98f8d2bc 100644
--- a/R/weighted.R
+++ b/R/weighted.R
@@ -42,8 +42,7 @@ weighted_quantile = function(x,
)(probs)
}
-#' @rdname weighted_quantile
-#' @export
+#' @noRd
weighted_quantile_fun = function(x, weights = NULL, na.rm = FALSE, type = 7, ...) {
na.rm <- as_one_logical(na.rm)
assert_number(type, lower = 1, upper = 9)
diff --git a/README.Rmd b/README.Rmd
index e6716c17..6f771fa7 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -245,10 +245,12 @@ print(line)
### Contributing to posterior
-We welcome contributions! The **posterior** package is under active development.
+We welcome contributions!
+The **posterior** package is under active development.
If you find bugs or have ideas for new features (for us or yourself to
-implement) please open an issue on GitHub
-(https://github.com/stan-dev/posterior/issues).
+implement) please [open an issue](https://github.com/stan-dev/posterior/issues) on GitHub.
+See [CONTRIBUTING.md](https://github.com/stan-dev/posterior/blob/master/.github/CONTRIBUTING.md)
+for more details.
### Citing posterior
@@ -263,7 +265,8 @@ When using **posterior**, please cite it as follows:
Working with Posterior Distributions.” R package version XXX, This outlines how to propose a change to posterior and is based on similar instructions for tidyverse packages, including the contributing guidelines generated by You can fix typos, spelling mistakes, or grammatical errors in the documentation directly using the GitHub web interface, as long as the changes are made in the source file. This generally means you’ll need to edit roxygen2 comments in an If you want to make a bigger change, it’s a good idea to first file an issue and make sure someone from the team agrees that it’s needed. If you’ve found a bug, please file an issue that illustrates the bug with a minimal reproducible example (see e.g. the tidyverse reprex instructions). The tidyverse guide on how to create a great issue has more advice. If you are new to creating pull requests here are some tips. Using the functions from the Fork the package and clone onto your computer. If you haven’t done this before, we recommend using Install all development dependencies with Create a Git branch for your pull request (PR). We recommend using Make your changes, commit to git, and then create a PR by running For user-facing changes, add a bullet to the top of New code should attempt to follow the style used in the package. When in doubt follow the tidyverse style guide. We use roxygen2, with Markdown syntax, for documentation. We use testthat for unit tests. Contributions with test cases included are easier to accept. Please note that the posterior project follows the Stan project’s Code of Conduct. By contributing to this project you agree to abide by its terms. An overview of diagnostics related to Pareto-smoothed importance sampling The paper presents Pareto smoothed importance sampling, but also Pareto-\(\hat{k}\) diagnostic that can be used when
+estimating any expectation based on finite sample. This vignette
+illustrates the use of these diagnostics. The individual diagnostic
+functions are Additionally, the Generate Transform We examine the draws with the default
+ All the usual convergence diagnostics \(\widehat{R}\), Bulk-ESS, and Tail-ESS look
+good, which is fine as they have been designed to work also with
+infinite variance (Vehtari et al., 2020). If these variables would present variables of interest for which we
+would like to estimate means, we would be also interested in Monte Carlo
+standard error (MCSE, see case study How
+many iterations to run and how many digits to report). Here MCSE for mean is based on standard deviation and Basic-ESS, but
+these assume finite variance. We did sample also from distributions with
+infinite variance, but given a finite sample size, the empirical
+variance estimates are always finite, and thus we get overoptimistic
+MCSE. To diagnose whether our variables of interest may have infinite
+variance and even infinite mean, we can use Pareto-\(\hat{k}\) diagnostic. \(\hat{k} \leq 0\) indicates that
+all moments exist, and the inverse of positive \(\hat{k}\) tells estimate for the number of
+finite (fractional) moments. Thus, \(\hat{k}\geq 1/2\) indicates infinite
+variance, and \(\hat{k}\geq 1\)
+indicates infinite mean. Sometimes very thick distribution tails may
+affect also sampling, but assuming sampling did go well, and we would be
+interested only in quantiles, infinite variance and mean are not a
+problem. But if we are interested in mean, then we need to care about
+the number of (fractional) moments. Here we see \(\hat{k} \geq 1/2\) for \(t_2\), \(t_{1.5}\), and \(t_{1}\), and we should not trust their
+ If we really do need those mean estimates, we can improve
+trustworthiness by Pareto smoothing, which replaces extreme tail draws
+with expected ordered statistics of Pareto distribution fitted to the
+tails of the distribution. Pareto smoothed mean estimate (computed using
+Pareto smoothed draws) has finite variance with a cost of some bias
+which we know when it is negligible. As a thumb rule when \(\hat{k}<0.7\), the bias is
+negligible. We do Pareto smoothing for all the variables. Now the The bias and variance depend on the sample size, and we can use
+additional diagnostic Here required Given finite variance, the central limit theorem states that to halve
+MCSE we need four times bigger sample size. With Pareto smoothing, we
+can go further, but the convergence rate decreases when \(\hat{k}\) increases. We see that with \(t_2\), \(t_{1.5}\), and \(t_1\) we need \(4^{1/0.86}\approx 5\), \(4^{1/0.60}\approx 10\), and \(4^{1/0}\approx \infty\) times bigger sample
+sizes to halve MCSE for mean. The final Pareto diagnostic, \(\hat{k}\)-threshold, is useful for smaller
+sample sizes. Here we select only 100 iterations per chain to get total
+of 400 draws. With only 400 draws, we can trust the Pareto smoothed result only
+when \(\hat{k}<0.62\). For \(t_{1.5}\) \(\hat{k}\approx 0.64\), and
+ We can get all these diagnostics with All these diagnostics are presented in Section 3 and summarized in
+Table 1 in PSIS paper (Vehtari et al., 2024). If you don’t need to estimate means of thick tailed distributions,
+and there are no sampling issues due to thick tails, then you don’t need
+to check existence of finite variance, and thus there is no need to
+check Pareto-\(\hat{k}\) for all the
+parameters and derived quantities. It is possible that the distribution has finite variance, but
+pre-asymptotically given a finite sample size the behavior can be
+similar to infinite variance. Thus the diagnostic is useful even in
+cases where theory guarantees finite variance. Vehtari, A., Simpson, D., Gelman, A., Yao, Y., & Gabry, J.
+(2024). Pareto smoothed importance sampling. Journal of Machine
+Learning Research, 25(72):1-58. Vehtari A., Gelman A., Simpson D., Carpenter B., & Bürkner P. C.
+(2020). Rank-normalization, folding, and localization: An improved Rhat
+for assessing convergence of MCMC. Bayesian Analysis,
+16(2):667-718. Because the matrix was converted to a Analogous functions exist for the other draws formats and are used
similarly. In the call to Contributing to posterior
+ usethis::use_tidy_contributing().Fixing typos
+.R, not a .Rd file. You can find the .R file that generates the .Rd by reading the comment in the first line.Bigger changes
+Pull request process
+usethis package is not required but can be helpful if this process is new to you.usethis::create_from_github("stan-dev/posterior", fork = TRUE).devtools::install_dev_deps(), and then make sure the package passes R CMD check by running devtools::check(). If R CMD check doesn’t pass cleanly, it’s a good idea to ask for help before continuing.usethis::pr_init("brief-description-of-change").usethis::pr_push(), and following the prompts in your browser. The title of your PR should briefly describe the change. The body of your PR should contain Fixes #issue-number.NEWS.md (i.e. just below the first header). Follow the style already used in NEWS.md.Code style
+Code of Conduct
+rvar: The Random Variable Datatype
+ Pareto-khat diagnostics
+ Pareto-khat diagnostics
+ Aki
+Vehtari
+
+
+ Source: vignettes/pareto_diagnostics.Rmd
+ pareto_diagnostics.RmdIntroduction
+
+
+
+pareto_khat(), pareto_min_ss(),
+pareto_convergence_rate() and
+pareto_khat_threshold(). The function
+pareto_diags() will return all of these.pareto_smooth() function can be used
+to transform draws by smoothing the tail(s). ## Example
+## This is posterior version 1.6.0
+##
+## Attaching package: 'posterior'
+## The following objects are masked from 'package:stats':
+##
+## mad, sd, var
+
+## The following objects are masked from 'package:base':
+##
+## %in%, match
+##
+## Attaching package: 'dplyr'
+## The following objects are masked from 'package:stats':
+##
+## filter, lag
+## The following objects are masked from 'package:base':
+##
+## intersect, setdiff, setequal, union
+
options(pillar.neg = FALSE, pillar.subtle=FALSE, pillar.sigfig=2)Simulated data
+
+xn a simulated MCMC sample with 4 chains each
+with 1000 iterations using AR process with marginal normal(0,1)
+
N <- 1000
+phi <- 0.3
+set.seed(6534)
+dr <- array(data=replicate(4,as.numeric(arima.sim(n = N,
+ list(ar = c(phi)),
+ sd = sqrt((1-phi^2))))),
+ dim=c(N,4,1)) %>%
+ as_draws_df() %>%
+ set_variables('xn')xn via cdf-inverse-cdf so that we have
+variables that have marginally distributions \(t_3\), \(t_{2.5}\), \(t_2\), \(t_{1.5}\), and \(t_1\). These all have thick tails. In
+addition \(t_2\), \(t_{1.5}\), and \(t_1\) have infinite variance, and \(t_1\) (aka Cauchy) has infinite mean.MCMC convergence diagnostics
+
+summarise_draws().
+
drt %>%
+ summarise_draws()
+## # A tibble: 6 × 10
+## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
+## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
+## 1 xn 0.010 0.00094 0.99 0.99 -1.6 1.6 1.0 2284. 3189.
+## 2 xt3 0.025 0.0010 1.6 1.1 -2.3 2.3 1.0 2284. 3189.
+## 3 xt2_5 0.031 0.0010 2.0 1.1 -2.5 2.5 1.0 2284. 3189.
+## 4 xt2 0.046 0.0011 2.9 1.2 -2.8 2.9 1.0 2284. 3189.
+## 5 xt1_5 0.092 0.0011 7.6 1.3 -3.5 3.6 1.0 2284. 3189.
+## 6 xt1 0.33 0.0012 93. 1.5 -5.8 6.1 1.0 2284. 3189.
+
drt %>%
+ summarise_draws(mean, sd, mcse_mean, ess_bulk, ess_basic)
+## # A tibble: 6 × 6
+## variable mean sd mcse_mean ess_bulk ess_basic
+## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
+## 1 xn 0.010 0.99 0.021 2284. 2280.
+## 2 xt3 0.025 1.6 0.033 2284. 2452.
+## 3 xt2_5 0.031 2.0 0.039 2284. 2584.
+## 4 xt2 0.046 2.9 0.054 2284. 2903.
+## 5 xt1_5 0.092 7.6 0.13 2284. 3553.
+## 6 xt1 0.33 93. 1.5 2284. 3976.Pareto-\(\hat{k}\)
+
+
+
drt %>%
+ summarise_draws(mean, sd, mcse_mean, ess_basic, pareto_khat)
+## # A tibble: 6 × 6
+## variable mean sd mcse_mean ess_basic pareto_khat
+## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
+## 1 xn 0.010 0.99 0.021 2280. -0.072
+## 2 xt3 0.025 1.6 0.033 2452. 0.33
+## 3 xt2_5 0.031 2.0 0.039 2584. 0.41
+## 4 xt2 0.046 2.9 0.054 2903. 0.52
+## 5 xt1_5 0.092 7.6 0.13 3553. 0.69
+## 6 xt1 0.33 93. 1.5 3976. 1.0mcse_mean values. Without trustworthy MCSE estimate we
+don’t have good estimate of how accurate the mean estimate is.
+Furthermore, as \(\hat{k} \geq 1\) for
+\(t_{1}\), the mean is not finite and
+the mean estimate is not valid.Pareto smoothing
+
+
+
drts <- drt %>%
+ mutate_variables(xt3_s=pareto_smooth(xt3),
+ xt2_5_s=pareto_smooth(xt2_5),
+ xt2_s=pareto_smooth(xt2),
+ xt1_5_s=pareto_smooth(xt1_5),
+ xt1_s=pareto_smooth(xt1)) %>%
+ subset_draws(variable="_s", regex=TRUE)
+## Pareto k-hat = 0.32.
+## Pareto k-hat = 0.4.
+## Pareto k-hat = 0.51.
+## Pareto k-hat = 0.68.
+## Pareto k-hat = 1.02. Mean does not exist, making empirical mean estimate of the draws not applicable.mcse_mean values are more trustworthy when \(\hat{k} < 0.7\). When \(\hat{k}>0.7\) both bias and variance
+grow so fast that Pareto smoothing rarely helps (see more details in the
+paper).
+
drts %>%
+ summarise_draws(mean, mcse_mean, ess_basic, pareto_khat)
+## # A tibble: 5 × 5
+## variable mean mcse_mean ess_basic pareto_khat
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 xt3_s 0.026 0.033 2438. 0.33
+## 2 xt2_5_s 0.033 0.038 2536. 0.40
+## 3 xt2_s 0.052 0.051 2763. 0.50
+## 4 xt1_5_s 0.12 0.10 3293. 0.67
+## 5 xt1_s 0.97 0.80 3903. 0.98Minimum sample size required
+
+min_ss which tells the minimum sample
+size needed so that mcse_mean can be trusted.
+
drt %>%
+ summarise_draws(mean, mcse_mean, ess_basic,
+ pareto_khat, min_ss=pareto_min_ss)
+## # A tibble: 6 × 6
+## variable mean mcse_mean ess_basic pareto_khat min_ss
+## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
+## 1 xn 0.010 0.021 2280. -0.072 10
+## 2 xt3 0.025 0.033 2452. 0.33 31.
+## 3 xt2_5 0.031 0.039 2584. 0.41 48.
+## 4 xt2 0.046 0.054 2903. 0.52 116.
+## 5 xt1_5 0.092 0.13 3553. 0.69 1735.
+## 6 xt1 0.33 1.5 3976. 1.0 Infmin_ss is smaller than
+ess_basic for all except \(t_1\), for which there is no hope.Convergence rate
+
+
+
drt %>%
+ summarise_draws(mean, mcse_mean, ess_basic,
+ pareto_khat, min_ss=pareto_min_ss,
+ conv_rate=pareto_convergence_rate)
+## # A tibble: 6 × 7
+## variable mean mcse_mean ess_basic pareto_khat min_ss conv_rate
+## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
+## 1 xn 0.010 0.021 2280. -0.072 10 1
+## 2 xt3 0.025 0.033 2452. 0.33 31. 0.98
+## 3 xt2_5 0.031 0.039 2584. 0.41 48. 0.95
+## 4 xt2 0.046 0.054 2903. 0.52 116. 0.86
+## 5 xt1_5 0.092 0.13 3553. 0.69 1735. 0.60
+## 6 xt1 0.33 1.5 3976. 1.0 Inf 0Pareto-\(\hat{k}\)-threshold
+
+
+
drt %>%
+ subset_draws(iteration=1:100) %>%
+ summarise_draws(mean, mcse_mean, ess_basic,
+ pareto_khat, min_ss=pareto_min_ss,
+ khat_thres=pareto_khat_threshold,
+ conv_rate=pareto_convergence_rate)
+## # A tibble: 6 × 8
+## variable mean mcse_mean ess_basic pareto_khat min_ss khat_thres conv_rate
+## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
+## 1 xn 0.038 0.066 244. -0.012 1 e 1 0.62 1
+## 2 xt3 0.054 0.11 237. 0.32 3.0e 1 0.62 0.96
+## 3 xt2_5 0.057 0.13 237. 0.38 4.2e 1 0.62 0.93
+## 4 xt2 0.063 0.17 243. 0.48 8.3e 1 0.62 0.86
+## 5 xt1_5 0.061 0.31 271. 0.64 5.6e 2 0.62 0.66
+## 6 xt1 -0.26 1.6 344. 0.95 2.2e18 0.62 0.11min_ss reveals we would probably need more than 560 draws
+to be on the safe side.Pareto diagnostics
+
+pareto_diags(),
+and it’s easy to use it also for derived quantities.
+
drt %>%
+ mutate_variables(xt2_5_sq=xt2_5^2) %>%
+ subset_draws(variable="xt2_5_sq") %>%
+ summarise_draws(mean, mcse_mean,
+ pareto_diags)
+## # A tibble: 1 × 7
+## variable mean mcse_mean khat min_ss khat_threshold convergence_rate
+## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
+## 1 xt2_5_sq 3.9 0.56 0.67 1124. 0.72 0.63Discussion
+
+Reference
+
+Example
-
+## This is posterior version 1.5.0## This is posterior version 1.6.0##
## Attaching package: 'posterior'## The following objects are masked from 'package:stats':
@@ -339,17 +339,17 @@ Example: create draws_matrix
print(x)
+## draw V1 V2 V3 V4 V5
+## 1 -1.14 -1.386 0.54 -0.83 0.700
+## 2 -1.52 2.124 -0.35 0.64 1.078
+## 3 -2.01 -1.062 0.40 0.94 -0.071
+## 4 0.31 -0.032 1.69 0.55 -0.071
+## 5 -0.63 -1.580 0.83 0.69 -0.870
+## 6 0.12 0.687 1.76 -0.25 -0.448
+## 7 -0.25 0.278 0.79 0.68 2.251
+## 8 0.98 0.932 -1.02 0.45 2.170
+## 9 -0.30 -0.384 -0.12 1.34 0.434
+## 10 -0.38 -0.872 -2.81 -0.36 -0.379
## # A draws_matrix: 10 iterations, 1 chains, and 5 variables
## variable
-## draw V1 V2 V3 V4 V5
-## 1 0.123 -0.27 -0.59 0.60 -1.789
-## 2 0.353 0.24 -0.12 0.92 0.013
-## 3 -0.316 -0.60 -0.21 0.81 -1.121
-## 4 0.755 -1.38 -0.63 -0.17 -1.363
-## 5 0.042 0.88 0.66 0.28 -0.469
-## 6 -0.210 0.31 0.93 -0.80 -1.809
-## 7 0.135 0.63 0.81 -0.55 -2.086
-## 8 -0.615 0.18 -0.10 -0.34 0.790
-## 9 -0.372 0.13 -2.56 -1.22 1.520
-## 10 -0.354 0.33 0.60 0.79 0.153draws_matrix, all
of the methods for working with draws objects described in
subsequent sections of this vignette will now be available.Example: create draws
print(x)
## # A draws_matrix: 50 iterations, 1 chains, and 2 variables
## variable
-## draw alpha beta
-## 1 -0.7025 -0.72
-## 2 0.0350 0.44
-## 3 -0.5237 -2.36
-## 4 0.0415 -0.23
-## 5 -0.0051 1.01
-## 6 -0.2882 -1.28
-## 7 -1.5631 -0.11
-## 8 -0.4261 0.37
-## 9 -1.5636 1.00
-## 10 -0.2385 -1.28
+## draw alpha beta
+## 1 -0.044 -0.885
+## 2 -1.048 1.247
+## 3 0.983 0.725
+## 4 1.390 1.201
+## 5 -1.002 -0.693
+## 6 0.247 -1.470
+## 7 1.265 -0.074
+## 8 0.113 -0.132
+## 9 -0.094 0.310
+## 10 -0.794 -0.014
## # ... with 40 more drawsRenaming
+variables(x)
# mu is a scalar, theta is a vector
x <- rename_variables(eight_schools_df, mean = mu, alpha = theta)
-variables(x)## [1] "mean" "tau" "alpha[1]" "alpha[2]" "alpha[3]" "alpha[4]"
## [7] "alpha[5]" "alpha[6]" "alpha[7]" "alpha[8]"rename_variables() above, mu
@@ -501,7 +501,7 @@ Renamingalpha:
x <- rename_variables(x, a1 = `alpha[1]`)
-variables(x)## [1] "mean" "tau" "a1" "alpha[2]" "alpha[3]" "alpha[4]"
## [7] "alpha[5]" "alpha[6]" "alpha[7]" "alpha[8]"
@@ -524,12 +524,12 @@ ## # A draws_matrix: 5 iterations, 1 chains, and 3 variables
## variable
-## draw alpha beta theta
-## 1 -0.83 -2.78 0.073
-## 2 -0.37 3.17 0.053
-## 3 0.21 -0.48 0.110
-## 4 0.31 -0.65 0.534
-## 5 -0.24 -1.51 0.694
+## draw alpha beta theta
+## 1 0.793 -0.073 0.22
+## 2 0.483 0.267 0.19
+## 3 0.331 0.238 6.60
+## 4 -0.837 1.301 0.47
+## 5 0.066 0.722 1.33
Because x1 and x2 have the same variables,
we can bind them along the 'draw' dimension to create a
single draws_matrix with more draws:
## # A draws_matrix: 10 iterations, 1 chains, and 2 variables
## variable
-## draw alpha beta
-## 1 -0.83 -2.78
-## 2 -0.37 3.17
-## 3 0.21 -0.48
-## 4 0.31 -0.65
-## 5 -0.24 -1.51
-## 6 -0.92 0.29
-## 7 1.09 0.51
-## 8 -0.15 -1.99
-## 9 -2.01 -0.63
-## 10 1.80 0.27
+## draw alpha beta
+## 1 0.793 -0.073
+## 2 0.483 0.267
+## 3 0.331 0.238
+## 4 -0.837 1.301
+## 5 0.066 0.722
+## 6 1.374 1.284
+## 7 -0.061 -0.212
+## 8 -0.386 -0.052
+## 9 -1.356 0.061
+## 10 0.627 -0.815
As with all posterior methods, bind_draws() can be used
with all draws formats and depending on the format different dimensions
are available to bind on. For example, we can bind
@@ -661,16 +661,16 @@
## # A tibble: 10 × 2
## variable weighted_mean
## <chr> <dbl>
-## 1 mu 4.00
-## 2 tau 4.09
-## 3 theta[1] 6.56
-## 4 theta[2] 4.95
-## 5 theta[3] 2.98
-## 6 theta[4] 4.72
-## 7 theta[5] 3.15
-## 8 theta[6] 3.59
-## 9 theta[7] 6.06
-## 10 theta[8] 4.45
+## 1 mu 4.13
+## 2 tau 4.07
+## 3 theta[1] 6.93
+## 4 theta[2] 5.54
+## 5 theta[3] 3.35
+## 6 theta[4] 5.03
+## 7 theta[5] 3.23
+## 8 theta[6] 3.52
+## 9 theta[7] 6.59
+## 10 theta[8] 4.58
Gelman A., Carlin J. B., Stern H. S., David B. Dunson D. B., Aki
-Vehtari A., & Rubin D. B. (2013). Bayesian Data Analysis, Third
+ Gelman A., Carlin J. B., Stern H. S., David B. Dunson D. B., Vehtari,
+A., & Rubin D. B. (2013). Bayesian Data Analysis, Third
Edition. Chapman and Hall/CRC. Vehtari A., Gelman A., Simpson D., Carpenter B., & Bürkner P. C.
(2020). Rank-normalization, folding, and localization: An improved Rhat
-for assessing convergence of MCMC. Bayesian Analysis.
vignettes/rvar.Rmd
rvar.Rmddata.frame()s and
-tibble()s, and to be used with distribution visualizations
+tibble()s, and to be used with distribution visualizations
in the ggdist
package.
@@ -342,12 +342,12 @@ draws_rvars datatypemu and
Sigma:
-variables(post)variables(post)
## [1] "mu" "Sigma"
But converted to a draws_list(), it contains one
variable for each combination of the dimensions of its variables:
-variables(as_draws_list(post))variables(as_draws_list(post))
## [1] "mu[1]" "mu[2]" "mu[3]" "Sigma[1,1]" "Sigma[2,1]"
## [6] "Sigma[3,1]" "Sigma[1,2]" "Sigma[2,2]" "Sigma[3,2]" "Sigma[1,3]"
## [11] "Sigma[2,3]" "Sigma[3,3]"
@@ -370,9 +370,12 @@ rvars## rvar<100,4>[3] mean ± sd:
## [1] 1.1 ± 0.11 1.1 ± 0.20 1.2 ± 0.31
Matrix multiplication is also implemented (using a tensor product
-under the hood). Because the normal matrix multiplication operator in R
-(%*%) cannot be properly implemented for S3 datatypes,
-rvar uses %**% instead. A trivial example:
%*%) cannot be properly implemented for S3
+datatypes, so rvar uses %**% instead. In R ≥
+4.3, which does support matrix multiplication for S3 datatypes, you can
+use %*% to matrix-multiply rvars.
+A trivial example:
## rvar<100,4>[3,3] mean ± sd:
@@ -420,7 +423,8 @@ Math with rvars
Matrix multiplication
-%**%
+
+%**%, %*% (R ≥ 4.3 only)
Basic functions
@@ -935,9 +939,10 @@ Subsetting rvars by
## [1] 4 ± 2
The resulting mixture looks like this:
+library(ggplot2)
+## Warning: package 'ggplot2' was built under R version 4.2.3
+

See vignette("slabinterval", package = "ggdist") for
more examples of visualizing distribution-type objects, including
@@ -953,7 +958,7 @@
Conditionals using rvar_ifelse()
yes where test == TRUE and draws from
no where test == FALSE.
Thus, we can create the mixture as follows:
-
+
x = rvar_ifelse(i == 1, component[[1]], component[[2]])
x
## rvar<4000>[1] mean ± sd:
@@ -972,7 +977,7 @@ Selecting diffe
rvar whose values are either 1 or
2 within each draw, you can use it as an index directly on
component to create the mixture:
-
+
x = component[[i]]
x
## rvar<4000>[1] mean ± sd:
@@ -993,7 +998,7 @@ Applying functions over rvar
its first dimension, which may be necessary for compatibility with some
functions (like purrr:map()).
For example, given this multidimensional rvar…
-
+
set.seed(3456)
x <- rvar_rng(rnorm, 24, mean = 1:24)
dim(x) <- c(2,3,4)
@@ -1024,7 +1029,7 @@ Applying functions over rvar
## [2,] 20 ± 1.00 22 ± 1.00 24 ± 1.00
… you can apply functions along the margins using
apply() (here, a silly example):
-
+
## [,1] [,2] [,3]
## [1,] 4 4 4
@@ -1040,12 +1045,12 @@ Applying functions over rvar
For example, you can use rvar_apply() with
rvar_mean() to compute the distributions of means along one
margin of an array:
-
+
rvar_apply(x, 1, rvar_mean)
## rvar<4000>[2] mean ± sd:
## [1] 12 ± 0.29 13 ± 0.29
Or along multiple dimensions:
-
+
rvar_apply(x, c(2,3), rvar_mean)
## rvar<4000>[3,4] mean ± sd:
## [,1] [,2] [,3] [,4]
@@ -1068,7 +1073,7 @@ Looping over draws and
base-R plots of individual draws (for ggplot2-based
plotting of rvars, see the next section and the ggdist package). For
example, it can be used to construct a parallel coordinates plot:
-
+
eight_schools <- as_draws_rvars(example_draws())
plot(1, type = "n",
@@ -1097,8 +1102,8 @@ Looping over draws and
Using rvars in data frames and in ggplot2
rvars can be used as columns in
-data.frame() or tibble() objects:
-
+data.frame() or tibble() objects:
+
df <- data.frame(group = c("a","b","c","d"), mu)
df
## group mu
@@ -1113,11 +1118,12 @@ Using rvars in d
stat_... family of geometries in the ggdist package, such as
stat_halfeye(), stat_lineribbon(), and
stat_dotsinterval(). For example:
-
+
+## Warning: package 'ggdist' was built under R version 4.2.3
+
+ggplot(df) +
stat_halfeye(aes(y = group, xdist = mu))

See vignette("slabinterval", package = "ggdist") or
diff --git a/docs/articles/rvar_files/figure-html/data_frame_plot-1.png b/docs/articles/rvar_files/figure-html/data_frame_plot-1.png
index 4cf02c0b..30e22dfe 100644
Binary files a/docs/articles/rvar_files/figure-html/data_frame_plot-1.png and b/docs/articles/rvar_files/figure-html/data_frame_plot-1.png differ
diff --git a/docs/articles/rvar_files/figure-html/mixture-1.png b/docs/articles/rvar_files/figure-html/mixture-1.png
index a63f63ad..add78e0e 100644
Binary files a/docs/articles/rvar_files/figure-html/mixture-1.png and b/docs/articles/rvar_files/figure-html/mixture-1.png differ
diff --git a/docs/authors.html b/docs/authors.html
index f3e2af9b..acb26f56 100644
--- a/docs/authors.html
+++ b/docs/authors.html
@@ -17,7 +17,7 @@
@@ -139,6 +139,10 @@ Authors
Jacob Socolar. Contributor.
+
+ Noa Kallioinen. Contributor.
+
+
@@ -148,15 +152,15 @@ Citation
- Bürkner P, Gabry J, Kay M, Vehtari A (2023).
+
Bürkner P, Gabry J, Kay M, Vehtari A (2024).
“posterior: Tools for Working with Posterior Distributions.”
-R package version 1.5.0, https://mc-stan.org/posterior/.
+R package version 1.6.0, https://mc-stan.org/posterior/.
@Misc{,
title = {posterior: Tools for Working with Posterior Distributions},
author = {Paul-Christian Bürkner and Jonah Gabry and Matthew Kay and Aki Vehtari},
- year = {2023},
- note = {R package version 1.5.0},
+ year = {2024},
+ note = {R package version 1.6.0},
url = {https://mc-stan.org/posterior/},
}
Vehtari A, Gelman A, Simpson D, Carpenter B, Bürkner P (2021).
diff --git a/docs/index.html b/docs/index.html
index ae75cc26..714c5c83 100644
--- a/docs/index.html
+++ b/docs/index.html
@@ -50,7 +50,7 @@
@@ -328,7 +328,7 @@ Mutating and renaming draws
x <- rename_variables(eight_schools_df, mean = mu, alpha = theta)
-variables(x)
+variables(x)
#> [1] "mean" "tau" "alpha[1]" "alpha[2]" "alpha[3]" "alpha[4]" "alpha[5]"
#> [8] "alpha[6]" "alpha[7]" "alpha[8]"
As with all posterior methods, mutate_variables and rename_variables can be used with all draws formats.
@@ -409,7 +409,7 @@ Converting from regu
x <- array(data=rnorm(200), dim=c(10, 2, 5))
x <- as_draws_matrix(x)
-variables(x) <- paste0("V", 1:5)
+variables(x) <- paste0("V", 1:5)
print(x)
#> # A draws_matrix: 10 iterations, 2 chains, and 5 variables
#> variable
@@ -480,7 +480,7 @@ Licensinghttps://opensource.org/license/bsd-3-clause/)
+
Code: BSD 3-clause (https://opensource.org/license/bsd-3-clause)
Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)
@@ -506,6 +506,12 @@ License
+
+Community
+
+
Citation
diff --git a/docs/news/index.html b/docs/news/index.html
index 9d738f3f..d94aa8a1 100644
--- a/docs/news/index.html
+++ b/docs/news/index.html
@@ -17,7 +17,7 @@
@@ -102,6 +102,23 @@ Changelog
Source: NEWS.md
+
+posterior 1.6.02024-07-03
+
+Enhancements
+- Add
exclude option to subset_draws(), which can be used to exclude the matched selection.
+- Add
are_log_weights option to pareto_smooth(), which is necessary for correct Pareto smoothing computation if the input vector consists of log weights.
+- Add
pareto_smooth option to weight_draws(), to Pareto smooth weights before adding to a draws object.
+- Add individual Pareto diagnostic functions (
pareto_khat(), pareto_khat_threshold(), pareto_min_ss(), pareto_convergence_rate())
+-
+
thin_draws() now automatically thins draws based on ESS by default, and non-integer thinning is possible.
+- Matrix multiplication of
rvars can now be done with the base matrix multiplication operator (%*%) instead of %**% in R >= 4.3.
+-
+
variables(), variables<-(), set_variables(), and nvariables() now support a with_indices argument, which determines whether variable names are retrieved/set with ("x[1]", "x[2]" …) or without ("x") indices (#208).
+- Add
extract_variable_array() function to extract variables with indices into arrays of iterations x chains x any remaining dimensions (#340).
+- For types that support
factor variables (draws_df, draws_list, and draws_rvars), extract_variable() and extract_variable_matrix() can now return factors.
+
+
posterior 1.5.02023-10-31
diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml
index 8ad1e3a7..ef76ad77 100644
--- a/docs/pkgdown.yml
+++ b/docs/pkgdown.yml
@@ -1,10 +1,11 @@
-pandoc: 3.1.1
+pandoc: 3.1.11
pkgdown: 2.0.7
pkgdown_sha: ~
articles:
+ pareto_diagnostics: pareto_diagnostics.html
posterior: posterior.html
rvar: rvar.html
-last_built: 2023-11-10T09:53Z
+last_built: 2024-07-05T08:44Z
urls:
reference: https://mc-stan.org/posterior/reference
article: https://mc-stan.org/posterior/articles
diff --git a/docs/pull_request_template.html b/docs/pull_request_template.html
index 4bc6c54d..d7bc28be 100644
--- a/docs/pull_request_template.html
+++ b/docs/pull_request_template.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/as_rvar.html b/docs/reference/as_rvar.html
index e0cce8e9..b5fe7754 100644
--- a/docs/reference/as_rvar.html
+++ b/docs/reference/as_rvar.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/as_rvar_factor.html b/docs/reference/as_rvar_factor.html
index 31187484..417130e0 100644
--- a/docs/reference/as_rvar_factor.html
+++ b/docs/reference/as_rvar_factor.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/autocorrelation.html b/docs/reference/autocorrelation.html
index 46947ec1..dac6b794 100644
--- a/docs/reference/autocorrelation.html
+++ b/docs/reference/autocorrelation.html
@@ -19,7 +19,7 @@
diff --git a/docs/reference/autocovariance.html b/docs/reference/autocovariance.html
index 3faaf667..76e996eb 100644
--- a/docs/reference/autocovariance.html
+++ b/docs/reference/autocovariance.html
@@ -19,7 +19,7 @@
diff --git a/docs/reference/bind_draws.html b/docs/reference/bind_draws.html
index 0d7a9cc1..937411cc 100644
--- a/docs/reference/bind_draws.html
+++ b/docs/reference/bind_draws.html
@@ -17,7 +17,7 @@
@@ -165,7 +165,7 @@ Examples
x4 <- draws_matrix(theta = rexp(5))
x5 <- bind_draws(x1, x4, along = "variable")
-variables(x5)
+variables(x5)
#> [1] "alpha" "beta" "theta"
diff --git a/docs/reference/chol.rvar.html b/docs/reference/chol.rvar.html
index 0ec7c5fd..7c17b693 100644
--- a/docs/reference/chol.rvar.html
+++ b/docs/reference/chol.rvar.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/diag-rvar-method.html b/docs/reference/diag-rvar-method.html
index a7d4690a..467cd5b3 100644
--- a/docs/reference/diag-rvar-method.html
+++ b/docs/reference/diag-rvar-method.html
@@ -18,7 +18,7 @@
diff --git a/docs/reference/diagnostics.html b/docs/reference/diagnostics.html
index 186cec0c..6c5404e6 100644
--- a/docs/reference/diagnostics.html
+++ b/docs/reference/diagnostics.html
@@ -17,7 +17,7 @@
@@ -117,7 +117,7 @@ Value
Details
-Function Description ess_basic()Basic version of effective sample size ess_bulk()Bulk effective sample size ess_tail()Tail effective sample size ess_quantile()Effective sample sizes for quantiles ess_sd()Effective sample sizes for standard deviations mcse_mean()Monte Carlo standard error for the mean mcse_quantile()Monte Carlo standard error for quantiles mcse_sd()Monte Carlo standard error for standard deviations rhat_basic()Basic version of Rhat rhat()Improved, rank-based version of Rhat rhat_nested()Rhat for use with many short chains rstar()R* diagnostic
+Function Description ess_basic()Basic version of effective sample size ess_bulk()Bulk effective sample size ess_tail()Tail effective sample size ess_quantile()Effective sample sizes for quantiles ess_sd()Effective sample sizes for standard deviations mcse_mean()Monte Carlo standard error for the mean mcse_quantile()Monte Carlo standard error for quantiles mcse_sd()Monte Carlo standard error for standard deviations pareto_khat()Pareto khat diagnostic for tail(s) pareto_diags()Additional diagnostics related to Pareto khat rhat_basic()Basic version of Rhat rhat()Improved, rank-based version of Rhat rhat_nested()Rhat for use with many short chains rstar()R* diagnostic
diff --git a/docs/reference/draws-index.html b/docs/reference/draws-index.html
index 004fadd1..3658f418 100644
--- a/docs/reference/draws-index.html
+++ b/docs/reference/draws-index.html
@@ -1,5 +1,5 @@
-Index draws objects — draws-index • posterior Index draws objects — draws-index • posterior
@@ -17,7 +17,7 @@
@@ -104,22 +104,16 @@ Index draws objects
- Index variables, iterations, chains, and draws.
+ Index iterations, chains, and draws of draws objects.
- variables(x, ...)
-
-variables(x) <- value
-
-iteration_ids(x)
+ iteration_ids(x)
chain_ids(x)
draw_ids(x)
-nvariables(x, ...)
-
niterations(x)
nchains(x)
@@ -133,23 +127,11 @@ Arguments
(draws) A draws object or another R object for which the method
is defined.
-
-...
-Arguments passed to individual methods (if applicable).
-
-
-value
-(character vector) For variables(x) <- value, the new variable
-names to use.
-
Value
-For variables(), a character vector.
-
-
For iteration_ids(), chain_ids(), and draw_ids(), an integer vector.
@@ -157,27 +139,21 @@ Value
Details
- The methods variables(), iteration_ids(), chain_ids(), and draw_ids() return
-vectors of all variables, iterations, chains, and draws, respectively. In
-contrast, the methods nvariables(), niterations(), nchains(), and
+
The methods iteration_ids(), chain_ids(), and draw_ids() return
+vectors of all iterations, chains, and draws, respectively. In
+contrast, the methods niterations(), nchains(), and
ndraws() return the number of variables, iterations, chains, and draws,
respectively.
-variables(x) <- value allows you to modify the vector of variable names,
-similar to how names(x) <- value works for vectors and lists. For renaming
-specific variables, set_variables() works equivalently, but is more intuitive when using the pipe operator. rename_variables() may offer a more convenient approach.
+
+
+ See also
+
Examples
x <- example_draws()
-variables(x)
-#> [1] "mu" "tau" "theta[1]" "theta[2]" "theta[3]" "theta[4]"
-#> [7] "theta[5]" "theta[6]" "theta[7]" "theta[8]"
-nvariables(x)
-#> [1] 10
-variables(x) <- letters[1:nvariables(x)]
-
iteration_ids(x)
#> [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
#> [19] 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
diff --git a/docs/reference/draws.html b/docs/reference/draws.html
index 74a4b3b7..9bd04ab8 100644
--- a/docs/reference/draws.html
+++ b/docs/reference/draws.html
@@ -18,7 +18,7 @@
diff --git a/docs/reference/draws_array.html b/docs/reference/draws_array.html
index 29dc038a..44e66e6d 100644
--- a/docs/reference/draws_array.html
+++ b/docs/reference/draws_array.html
@@ -21,7 +21,7 @@
@@ -179,11 +179,11 @@ Details
See also
Other formats:
+draws,
draws_df(),
draws_list(),
draws_matrix(),
-draws_rvars(),
-draws
+draws_rvars()
diff --git a/docs/reference/draws_df.html b/docs/reference/draws_df.html
index 19f52749..cee18fa3 100644
--- a/docs/reference/draws_df.html
+++ b/docs/reference/draws_df.html
@@ -21,7 +21,7 @@
@@ -189,11 +189,11 @@ Details
See also
Other formats:
+draws,
draws_array(),
draws_list(),
draws_matrix(),
-draws_rvars(),
-draws
+draws_rvars()
diff --git a/docs/reference/draws_list.html b/docs/reference/draws_list.html
index 3e8c53b4..b94afc4e 100644
--- a/docs/reference/draws_list.html
+++ b/docs/reference/draws_list.html
@@ -21,7 +21,7 @@
@@ -181,11 +181,11 @@ Details
See also
Other formats:
+draws,
draws_array(),
draws_df(),
draws_matrix(),
-draws_rvars(),
-draws
+draws_rvars()
diff --git a/docs/reference/draws_matrix.html b/docs/reference/draws_matrix.html
index 8b79b02f..8ebdabfe 100644
--- a/docs/reference/draws_matrix.html
+++ b/docs/reference/draws_matrix.html
@@ -21,7 +21,7 @@
@@ -179,11 +179,11 @@ Details
See also
Other formats:
+draws,
draws_array(),
draws_df(),
draws_list(),
-draws_rvars(),
-draws
+draws_rvars()
diff --git a/docs/reference/draws_of.html b/docs/reference/draws_of.html
index bb599298..81420a41 100644
--- a/docs/reference/draws_of.html
+++ b/docs/reference/draws_of.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/draws_rvars.html b/docs/reference/draws_rvars.html
index bb8ede22..188a1371 100644
--- a/docs/reference/draws_rvars.html
+++ b/docs/reference/draws_rvars.html
@@ -21,7 +21,7 @@
@@ -182,11 +182,11 @@ Details
See also
Other formats:
+draws,
draws_array(),
draws_df(),
draws_list(),
-draws_matrix(),
-draws
+draws_matrix()
diff --git a/docs/reference/draws_summary.html b/docs/reference/draws_summary.html
index e2416de6..9b1ffd29 100644
--- a/docs/reference/draws_summary.html
+++ b/docs/reference/draws_summary.html
@@ -22,7 +22,7 @@
diff --git a/docs/reference/drop-rvar-method.html b/docs/reference/drop-rvar-method.html
index 7370fa52..8cc9f12b 100644
--- a/docs/reference/drop-rvar-method.html
+++ b/docs/reference/drop-rvar-method.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/entropy.html b/docs/reference/entropy.html
index 424579ff..92d2abe9 100644
--- a/docs/reference/entropy.html
+++ b/docs/reference/entropy.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/ess_basic.html b/docs/reference/ess_basic.html
index 52b65cdf..7b11c10f 100644
--- a/docs/reference/ess_basic.html
+++ b/docs/reference/ess_basic.html
@@ -22,7 +22,7 @@
@@ -189,9 +189,11 @@ See also
mcse_mean(),
mcse_quantile(),
mcse_sd(),
+pareto_diags(),
+pareto_khat(),
+rhat(),
rhat_basic(),
rhat_nested(),
-rhat(),
rstar()
diff --git a/docs/reference/ess_bulk.html b/docs/reference/ess_bulk.html
index aad47857..b3ef7cd1 100644
--- a/docs/reference/ess_bulk.html
+++ b/docs/reference/ess_bulk.html
@@ -22,7 +22,7 @@
@@ -181,9 +181,11 @@ See also
mcse_mean(),
mcse_quantile(),
mcse_sd(),
+pareto_diags(),
+pareto_khat(),
+rhat(),
rhat_basic(),
rhat_nested(),
-rhat(),
rstar()
diff --git a/docs/reference/ess_mean.html b/docs/reference/ess_mean.html
index 93a70275..0b5ff9e1 100644
--- a/docs/reference/ess_mean.html
+++ b/docs/reference/ess_mean.html
@@ -18,7 +18,7 @@
diff --git a/docs/reference/ess_quantile.html b/docs/reference/ess_quantile.html
index 98ca6c63..287678ab 100644
--- a/docs/reference/ess_quantile.html
+++ b/docs/reference/ess_quantile.html
@@ -18,7 +18,7 @@
@@ -190,9 +190,11 @@ See also
mcse_mean(),
mcse_quantile(),
mcse_sd(),
+pareto_diags(),
+pareto_khat(),
+rhat(),
rhat_basic(),
rhat_nested(),
-rhat(),
rstar()
diff --git a/docs/reference/ess_sd.html b/docs/reference/ess_sd.html
index 0ffc275b..4b186826 100644
--- a/docs/reference/ess_sd.html
+++ b/docs/reference/ess_sd.html
@@ -19,7 +19,7 @@
@@ -173,9 +173,11 @@ See also
mcse_mean(),
mcse_quantile(),
mcse_sd(),
+pareto_diags(),
+pareto_khat(),
+rhat(),
rhat_basic(),
rhat_nested(),
-rhat(),
rstar()
diff --git a/docs/reference/ess_tail.html b/docs/reference/ess_tail.html
index 4ee43339..c48cbfa7 100644
--- a/docs/reference/ess_tail.html
+++ b/docs/reference/ess_tail.html
@@ -22,7 +22,7 @@
@@ -181,9 +181,11 @@ See also
mcse_mean(),
mcse_quantile(),
mcse_sd(),
+pareto_diags(),
+pareto_khat(),
+rhat(),
rhat_basic(),
rhat_nested(),
-rhat(),
rstar()
diff --git a/docs/reference/example_draws.html b/docs/reference/example_draws.html
index 31881809..086226df 100644
--- a/docs/reference/example_draws.html
+++ b/docs/reference/example_draws.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/extract_variable.html b/docs/reference/extract_variable.html
index 6b278ec6..9a75d9ae 100644
--- a/docs/reference/extract_variable.html
+++ b/docs/reference/extract_variable.html
@@ -17,7 +17,7 @@
@@ -116,6 +116,12 @@ Extract draws of a single variable
# S3 method for draws
extract_variable(x, variable, ...)
+# S3 method for draws_df
+extract_variable(x, variable, ...)
+
+# S3 method for draws_list
+extract_variable(x, variable, ...)
+
# S3 method for draws_rvars
extract_variable(x, variable, ...)
@@ -128,7 +134,9 @@ Arguments
variable
-(string) The name of the variable to extract.
+(string) The name of the variable to extract. Must include
+indices for array variables (e.g. "x[1]", "y[1,2]"). To extract all
+dimensions from variables with indices, use extract_variable_array().
...
@@ -139,7 +147,13 @@ Arguments
Value
-
A numeric vector of length equal to the number of draws.
+A vector of length equal to the number of draws.
+
+
+ See also
+ Other variable extraction methods:
+extract_variable_array(),
+extract_variable_matrix()
diff --git a/docs/reference/extract_variable_array.html b/docs/reference/extract_variable_array.html
new file mode 100644
index 00000000..a0b3bd7d
--- /dev/null
+++ b/docs/reference/extract_variable_array.html
@@ -0,0 +1,207 @@
+
+Extract array of a single (possibly indexed) variable — extract_variable_array • posterior
+
+
+
+
+
+
+
+
+
+
+ Extract array of a single (possibly indexed) variable
+ Source: R/extract_variable_array.R
+ extract_variable_array.Rd
+
+
+
+ Extract an array of draws of a single variable, including any dimensions of
+variables with indices.
+
+
+
+ extract_variable_array(x, variable, ...)
+
+# S3 method for default
+extract_variable_array(x, variable, ...)
+
+# S3 method for draws
+extract_variable_array(x, variable, ...)
+
+
+
+ Arguments
+ - x
+(draws) A draws object or another R object for which the method
+is defined.
+
+
+- variable
+(string) The name of the variable to extract. To extract all
+dimensions from variables with indices (e.g. "x[1]"), provide the base
+variable name (e.g. "x").
+
+
+- ...
+Arguments passed to individual methods (if applicable).
+
+
+
+ Value
+
+
+An array with dimension niterations(x) x nchains(x) x any remaining
+dimensions determined by the indices of the variable x.
+
+
+ See also
+ Other variable extraction methods:
+extract_variable(),
+extract_variable_matrix()
+
+
+
+ Examples
+ x <- example_draws(example = "multi_normal")
+
+mu <- extract_variable_array(x, variable = "mu")
+str(mu)
+#> num [1:100, 1:4, 1:3] 0.18119 -0.03419 -0.05875 -0.1536 0.00989 ...
+#> - attr(*, "dimnames")=List of 3
+#> ..$ : NULL
+#> ..$ : NULL
+#> ..$ : NULL
+
+mu1 <- extract_variable_array(x, variable = "mu[1]")
+str(mu1)
+#> num [1:100, 1:4, 1] 0.18119 -0.03419 -0.05875 -0.1536 0.00989 ...
+#> - attr(*, "dimnames")=List of 3
+#> ..$ : NULL
+#> ..$ : NULL
+#> ..$ : NULL
+
+Sigma <- extract_variable_array(x, variable = "Sigma")
+str(Sigma)
+#> num [1:100, 1:4, 1:3, 1:3] 1.2 1.14 1.12 1.14 1.19 ...
+#> - attr(*, "dimnames")=List of 4
+#> ..$ : NULL
+#> ..$ : NULL
+#> ..$ : NULL
+#> ..$ : NULL
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/docs/reference/extract_variable_matrix.html b/docs/reference/extract_variable_matrix.html
index 186d9b8d..2be61fc4 100644
--- a/docs/reference/extract_variable_matrix.html
+++ b/docs/reference/extract_variable_matrix.html
@@ -18,7 +18,7 @@
@@ -118,6 +118,12 @@ Extract matrix of a single variable
# S3 method for draws
extract_variable_matrix(x, variable, ...)
+# S3 method for draws_df
+extract_variable_matrix(x, variable, ...)
+
+# S3 method for draws_list
+extract_variable_matrix(x, variable, ...)
+
# S3 method for draws_rvars
extract_variable_matrix(x, variable, ...)
@@ -130,7 +136,9 @@ Arguments
variable
-(string) The name of the variable to extract.
+(string) The name of the variable to extract. Must include
+indices for array variables (e.g. "x[1]", "y[1,2]"). To extract all
+dimensions from variables with indices, use extract_variable_array().
...
@@ -143,6 +151,12 @@ Value
A matrix with dimension iterations x chains.
+
+ See also
+ Other variable extraction methods:
+extract_variable(),
+extract_variable_array()
+
Examples
diff --git a/docs/reference/for_each_draw.html b/docs/reference/for_each_draw.html
index 9e3acfdc..4d62ba8d 100644
--- a/docs/reference/for_each_draw.html
+++ b/docs/reference/for_each_draw.html
@@ -18,7 +18,7 @@
@@ -124,7 +124,8 @@ Arguments
(expression) A bare expression that can contain references to
variables in x by name. This expression will be executed once per draw
of x, where references to variables in x resolve to the value of that
-variable in that draw. The expression supports quasiquotation.
+variable in that draw. The expression supports
+quasiquotation.
diff --git a/docs/reference/index.html b/docs/reference/index.html
index baab84ab..f5dc5818 100644
--- a/docs/reference/index.html
+++ b/docs/reference/index.html
@@ -17,7 +17,7 @@
@@ -166,9 +166,9 @@ Draws objects and formats variables() `variables<-`() iteration_ids() chain_ids() draw_ids() nvariables() niterations() nchains() ndraws()
+
- Index draws objects
+ Get variable names from draws objects
@@ -189,6 +189,10 @@ Working with draws objects extract_variable()
Extract draws of a single variable
+
+
+
+ Extract array of a single (possibly indexed) variable
@@ -202,7 +206,7 @@ Working with draws objects set_variables()
+
`variables<-`() set_variables()
Set variable names in draws objects
@@ -217,6 +221,10 @@ Working with draws objects subset_draws() subset(<draws>)
Subset draws objects
+
+ iteration_ids() chain_ids() draw_ids() niterations() nchains() ndraws()
+
+ Index draws objects
@@ -301,6 +309,34 @@ Summarizing and diagnosing dra
Monte Carlo standard error for the standard deviation
+
+ pareto_diags() pareto_khat_threshold() pareto_min_ss() pareto_convergence_rate()
+
+ Pareto smoothing diagnostics
+
+
+
+ Pareto khat diagnostic
+
+
+
+ Pareto smoothing
+
+
+
+ Pareto convergence rate
+
+
+
+ Pareto k-hat threshold
+
+
+
+ Pareto-smoothing minimum sample-size
+
+
+
+ Pareto tail length
@@ -321,18 +357,6 @@ Summarizing and diagnosing dra
Modal category
-
-
-
- Pareto smoothing diagnostics
-
-
-
- Pareto khat diagnostic
-
-
-
- Pareto smoothing
Functionality specific to the rvar datatype
The draws_rvar format (a structured list of rvar objects) has the same methods (e.g. bind_draws()) as the other draws formats. For individual rvar objects themselves, however, posterior provides additional functionality.
@@ -346,7 +370,7 @@ Functionality specific to t
Density, CDF, and quantile functions of random variables
-
+
Matrix multiplication of random variables
diff --git a/docs/reference/is_rvar.html b/docs/reference/is_rvar.html
index 4724eba3..f19f07f4 100644
--- a/docs/reference/is_rvar.html
+++ b/docs/reference/is_rvar.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/is_rvar_factor.html b/docs/reference/is_rvar_factor.html
index a771634d..38de9e33 100644
--- a/docs/reference/is_rvar_factor.html
+++ b/docs/reference/is_rvar_factor.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/match.html b/docs/reference/match.html
index c18dd181..42adcb70 100644
--- a/docs/reference/match.html
+++ b/docs/reference/match.html
@@ -19,7 +19,7 @@
diff --git a/docs/reference/mcse_mean.html b/docs/reference/mcse_mean.html
index c51f4260..d7965430 100644
--- a/docs/reference/mcse_mean.html
+++ b/docs/reference/mcse_mean.html
@@ -18,7 +18,7 @@
@@ -169,9 +169,11 @@ See also
ess_tail(),
mcse_quantile(),
mcse_sd(),
+pareto_diags(),
+pareto_khat(),
+rhat(),
rhat_basic(),
rhat_nested(),
-rhat(),
rstar()
diff --git a/docs/reference/mcse_quantile.html b/docs/reference/mcse_quantile.html
index 836681f7..c6760abe 100644
--- a/docs/reference/mcse_quantile.html
+++ b/docs/reference/mcse_quantile.html
@@ -18,7 +18,7 @@
@@ -187,9 +187,11 @@ See also
ess_tail(),
mcse_mean(),
mcse_sd(),
+pareto_diags(),
+pareto_khat(),
+rhat(),
rhat_basic(),
rhat_nested(),
-rhat(),
rstar()
diff --git a/docs/reference/mcse_sd.html b/docs/reference/mcse_sd.html
index 7907ad21..e1d0f975 100644
--- a/docs/reference/mcse_sd.html
+++ b/docs/reference/mcse_sd.html
@@ -19,7 +19,7 @@
@@ -174,9 +174,11 @@ See also
ess_tail(),
mcse_mean(),
mcse_quantile(),
+pareto_diags(),
+pareto_khat(),
+rhat(),
rhat_basic(),
rhat_nested(),
-rhat(),
rstar()
diff --git a/docs/reference/merge_chains.html b/docs/reference/merge_chains.html
index 1d65b227..041c85e0 100644
--- a/docs/reference/merge_chains.html
+++ b/docs/reference/merge_chains.html
@@ -21,7 +21,7 @@
diff --git a/docs/reference/modal_category.html b/docs/reference/modal_category.html
index 5e4a0421..4ba38bb7 100644
--- a/docs/reference/modal_category.html
+++ b/docs/reference/modal_category.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/mutate_variables.html b/docs/reference/mutate_variables.html
index 5a0abd21..7eeb05a7 100644
--- a/docs/reference/mutate_variables.html
+++ b/docs/reference/mutate_variables.html
@@ -17,7 +17,7 @@
@@ -160,7 +160,7 @@ Details
diff --git a/docs/reference/order_draws.html b/docs/reference/order_draws.html
index 8b024d90..d02aee2b 100644
--- a/docs/reference/order_draws.html
+++ b/docs/reference/order_draws.html
@@ -19,7 +19,7 @@
diff --git a/docs/reference/pareto_diags.html b/docs/reference/pareto_diags.html
index 7348b046..fb1240b6 100644
--- a/docs/reference/pareto_diags.html
+++ b/docs/reference/pareto_diags.html
@@ -19,7 +19,7 @@
@@ -121,11 +121,36 @@ Pareto smoothing diagnostics
r_eff = NULL,
ndraws_tail = NULL,
verbose = FALSE,
+ are_log_weights = FALSE,
...
)
# S3 method for rvar
-pareto_diags(x, ...)
+pareto_diags(x, ...)
+
+pareto_khat_threshold(x, ...)
+
+# S3 method for default
+pareto_khat_threshold(x, ...)
+
+# S3 method for rvar
+pareto_khat_threshold(x, ...)
+
+pareto_min_ss(x, ...)
+
+# S3 method for default
+pareto_min_ss(x, ...)
+
+# S3 method for rvar
+pareto_min_ss(x, ...)
+
+pareto_convergence_rate(x, ...)
+
+# S3 method for default
+pareto_convergence_rate(x, ...)
+
+# S3 method for rvar
+pareto_convergence_rate(x, ...)
@@ -150,15 +175,16 @@ Arguments
r_eff
(numeric) relative effective sample size estimate. If
-r_eff is omitted, it will be calculated assuming the draws are
-from MCMC.
+r_eff is NULL, it will be calculated assuming the draws are
+from MCMC. Default is NULL.
ndraws_tail
(numeric) number of draws for the tail. If
ndraws_tail is not specified, it will be calculated as
ceiling(3 * sqrt(length(x) / r_eff)) if length(x) > 225 and
-length(x) / 5 otherwise (see Appendix H in Vehtari et al. (2022)).
+length(x) / 5 otherwise (see Appendix H in Vehtari et
+al. (2024)).
verbose
@@ -166,6 +192,12 @@ Arguments
TRUE, messages related to Pareto diagnostics will be
printed. Default is FALSE.
+
+are_log_weights
+(logical) Are the draws log weights? Default is
+FALSE. If TRUE computation will take into account that the
+draws are log weights, and only right tail will be smoothed.
+
Value
@@ -186,25 +218,47 @@ Details
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).
+Section 3.2.3, Equation 11 in Vehtari et al. (2024).
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).
+in Vehtari et al. (2024).
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).
+theorem convergence rate. See Appendix B in Vehtari et al. (2024).
References
Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao and
-Jonah Gabry (2022). Pareto Smoothed Importance Sampling.
-arxiv:arXiv:1507.02646
+Jonah Gabry (2024). Pareto Smoothed Importance Sampling.
+Journal of Machine Learning Research, 25(72):1-58.
+PDF
+
+
+ See also
+ pareto_khat, pareto_min_ss,
+pareto_khat_threshold, and pareto_convergence_rate for
+individual diagnostics; and pareto_smooth for Pareto smoothing
+draws.
+Other diagnostics:
+ess_basic(),
+ess_bulk(),
+ess_quantile(),
+ess_sd(),
+ess_tail(),
+mcse_mean(),
+mcse_quantile(),
+mcse_sd(),
+pareto_khat(),
+rhat(),
+rhat_basic(),
+rhat_nested(),
+rstar()
diff --git a/docs/reference/pareto_khat.html b/docs/reference/pareto_khat.html
index 286f4195..f40686c0 100644
--- a/docs/reference/pareto_khat.html
+++ b/docs/reference/pareto_khat.html
@@ -2,7 +2,7 @@
Pareto khat diagnostic — pareto_khat • posterior
@@ -20,7 +20,7 @@
@@ -110,7 +110,7 @@ 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 estimate
the number of fractional moments that is useful for convergence
-diagnostics. For further details see Vehtari et al. (2022).
+diagnostics. For further details see Vehtari et al. (2024).
@@ -123,6 +123,7 @@ Pareto khat diagnostic
r_eff = NULL,
ndraws_tail = NULL,
verbose = FALSE,
+ are_log_weights = FALSE,
...
)
@@ -152,15 +153,16 @@ Arguments
r_eff
(numeric) relative effective sample size estimate. If
-r_eff is omitted, it will be calculated assuming the draws are
-from MCMC.
+r_eff is NULL, it will be calculated assuming the draws are
+from MCMC. Default is NULL.
ndraws_tail
(numeric) number of draws for the tail. If
ndraws_tail is not specified, it will be calculated as
ceiling(3 * sqrt(length(x) / r_eff)) if length(x) > 225 and
-length(x) / 5 otherwise (see Appendix H in Vehtari et al. (2022)).
+length(x) / 5 otherwise (see Appendix H in Vehtari et
+al. (2024)).
verbose
@@ -168,36 +170,73 @@ Arguments
TRUE, messages related to Pareto diagnostics will be
printed. Default is FALSE.
+
+are_log_weights
+(logical) Are the draws log weights? Default is
+FALSE. If TRUE computation will take into account that the
+draws are log weights, and only right tail will be smoothed.
+
Value
-khat estimated Generalized Pareto Distribution shape parameter k
+If the input is an array, returns a single numeric value. If any of the draws
+is non-finite, that is, NA, NaN, Inf, or -Inf, the returned output
+will be (numeric) NA. Also, if all draws within any of the chains of a
+variable are the same (constant), the returned output will be (numeric) NA
+
+
+as well. The reason for the latter is that, for constant draws, we cannot
+distinguish between variables that are supposed to be constant (e.g., a
+diagonal element of a correlation matrix is always 1) or variables that just
+happened to be constant because of a failure of convergence or other problems
+in the sampling process.
+
+
+If the input is an rvar, returns an array of the same dimensions as the
+rvar, where each element is equal to the value that would be returned by
+passing the draws array for that element of the rvar to this function.
References
Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao and
-Jonah Gabry (2022). Pareto Smoothed Importance Sampling.
-arxiv:arXiv:1507.02646
+Jonah Gabry (2024). Pareto Smoothed Importance Sampling.
+Journal of Machine Learning Research, 25(72):1-58.
+PDF
+
+
+ See also
+ pareto_diags for additional related diagnostics, and
+pareto_smooth for Pareto smoothed draws.
+Other diagnostics:
+ess_basic(),
+ess_bulk(),
+ess_quantile(),
+ess_sd(),
+ess_tail(),
+mcse_mean(),
+mcse_quantile(),
+mcse_sd(),
+pareto_diags(),
+rhat(),
+rhat_basic(),
+rhat_nested(),
+rstar()
Examples
mu <- extract_variable_matrix(example_draws(), "mu")
pareto_khat(mu)
-#> $khat
#> [1] 0.1979001
-#>
d <- as_draws_rvars(example_draws("multi_normal"))
pareto_khat(d$Sigma)
-#> $khat
#> [,1] [,2] [,3]
#> [1,] 0.05601935 0.04156719 0.05091481
#> [2,] 0.04156719 0.10157218 0.06191862
#> [3,] 0.05091481 0.06191862 -0.08123058
-#>
diff --git a/docs/reference/pareto_smooth.html b/docs/reference/pareto_smooth.html
index cbf7c19f..639bc4a8 100644
--- a/docs/reference/pareto_smooth.html
+++ b/docs/reference/pareto_smooth.html
@@ -1,7 +1,7 @@
Pareto smoothing — pareto_smooth • posterior
@@ -19,7 +19,7 @@
@@ -108,14 +108,14 @@ Pareto smoothing
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).
+tail(s). For further details see Vehtari et al. (2024).
pareto_smooth(x, ...)
# S3 method for rvar
-pareto_smooth(x, return_k = TRUE, extra_diags = FALSE, ...)
+pareto_smooth(x, return_k = FALSE, extra_diags = FALSE, ...)
# S3 method for default
pareto_smooth(
@@ -123,9 +123,10 @@ Pareto smoothing
tail = c("both", "right", "left"),
r_eff = NULL,
ndraws_tail = NULL,
- return_k = TRUE,
+ return_k = FALSE,
extra_diags = FALSE,
- verbose = FALSE,
+ verbose = TRUE,
+ are_log_weights = FALSE,
...
)
@@ -145,8 +146,9 @@ Arguments
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.
+output? If TRUE, output will be a list containing smoothed
+draws and diagnostics, otherwise it will be a numeric of the
+smoothed draws. Default is FALSE.
extra_diags
@@ -165,15 +167,16 @@ Arguments
r_eff
(numeric) relative effective sample size estimate. If
-r_eff is omitted, it will be calculated assuming the draws are
-from MCMC.
+r_eff is NULL, it will be calculated assuming the draws are
+from MCMC. Default is NULL.
ndraws_tail
(numeric) number of draws for the tail. If
ndraws_tail is not specified, it will be calculated as
ceiling(3 * sqrt(length(x) / r_eff)) if length(x) > 225 and
-length(x) / 5 otherwise (see Appendix H in Vehtari et al. (2022)).
+length(x) / 5 otherwise (see Appendix H in Vehtari et
+al. (2024)).
verbose
@@ -181,33 +184,50 @@ Arguments
TRUE, messages related to Pareto diagnostics will be
printed. Default is FALSE.
+
+are_log_weights
+(logical) Are the draws log weights? Default is
+FALSE. If TRUE computation will take into account that the
+draws are log weights, and only right tail will be smoothed.
+
Value
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
+containing the vector x and a named list diagnostics
+
+
+containing numeric values:
khat: estimated Pareto k shape parameter, and optionally
+min_ss: minimum sample size for reliable Pareto smoothed
+estimate
+khat_threshold: sample size specific khat threshold for
+reliable Pareto smoothed estimates
+convergence_rate: Relative convergence rate for
Pareto smoothed estimates
-convergence_rate: Relative convergence rate for Pareto smoothed estimates
-
+If any of the draws is non-finite, that is, NA, NaN, Inf, or
+-Inf, Pareto smoothing will not be performed, and the original
+draws will be returned and and diagnostics will be NA (numeric).
+
References
Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao and
-Jonah Gabry (2022). Pareto Smoothed Importance Sampling.
-arxiv:arXiv:1507.02646
+Jonah Gabry (2024). Pareto Smoothed Importance Sampling.
+Journal of Machine Learning Research, 25(72):1-58.
+PDF
+
+
+ See also
+ pareto_khat for only calculating khat, and
+pareto_diags for additional diagnostics.
Examples
mu <- extract_variable_matrix(example_draws(), "mu")
pareto_smooth(mu)
-#> $x
+#> Pareto k-hat = 0.2.
#> chain
#> iteration 1 2 3 4
#> 1 2.0058311 2.99038071 1.79436801 6.45897880
@@ -310,30 +330,23 @@ Examples
#> 98 6.1728283 1.51990114 0.15400719 2.72636415
#> 99 1.5485347 8.53166943 3.17470254 0.61142881
#> 100 7.5338964 -1.44854601 3.57555760 7.04723309
-#>
-#> $diagnostics
-#> $diagnostics$khat
-#> [1] 0.1979001
-#>
-#>
d <- as_draws_rvars(example_draws("multi_normal"))
pareto_smooth(d$Sigma)
-#> $x
+#> Pareto k-hat = 0.06.
+#> Pareto k-hat = 0.04.
+#> Pareto k-hat = 0.05.
+#> Pareto k-hat = 0.04.
+#> Pareto k-hat = 0.1.
+#> Pareto k-hat = 0.06.
+#> Pareto k-hat = 0.05.
+#> Pareto k-hat = 0.06.
+#> Pareto k-hat = -0.08.
#> rvar<100,4>[3,3] mean ± sd:
#> [,1] [,2] [,3]
#> [1,] 1.28 ± 0.17 0.53 ± 0.21 -0.40 ± 0.29
#> [2,] 0.53 ± 0.21 3.66 ± 0.45 -2.10 ± 0.49
#> [3,] -0.40 ± 0.29 -2.10 ± 0.49 8.12 ± 0.96
-#>
-#> $diagnostics
-#> $diagnostics$khat
-#> [,1] [,2] [,3]
-#> [1,] 0.05601935 0.04156719 0.05091481
-#> [2,] 0.04156719 0.10157218 0.06191862
-#> [3,] 0.05091481 0.06191862 -0.08123058
-#>
-#>
diff --git a/docs/reference/posterior-package.html b/docs/reference/posterior-package.html
index bf384314..fda7fc4b 100644
--- a/docs/reference/posterior-package.html
+++ b/docs/reference/posterior-package.html
@@ -31,7 +31,7 @@
@@ -161,6 +161,26 @@ Package options
trigger an automatic merging of chains, for example, because chains do not
match between two objects involved in a binary operation. Whether this
causes a warning can be controlled by this option.
+
+
+ See also
+ Useful links:
+
+
diff --git a/docs/reference/print.draws_array.html b/docs/reference/print.draws_array.html
index 1c65ef7d..c26582ad 100644
--- a/docs/reference/print.draws_array.html
+++ b/docs/reference/print.draws_array.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/print.draws_df.html b/docs/reference/print.draws_df.html
index 0472145a..9b95be58 100644
--- a/docs/reference/print.draws_df.html
+++ b/docs/reference/print.draws_df.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/print.draws_list.html b/docs/reference/print.draws_list.html
index 81574f4d..4ee47f81 100644
--- a/docs/reference/print.draws_list.html
+++ b/docs/reference/print.draws_list.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/print.draws_matrix.html b/docs/reference/print.draws_matrix.html
index 42bf6225..85db03e0 100644
--- a/docs/reference/print.draws_matrix.html
+++ b/docs/reference/print.draws_matrix.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/print.draws_rvars.html b/docs/reference/print.draws_rvars.html
index 2e2d42d9..3620971d 100644
--- a/docs/reference/print.draws_rvars.html
+++ b/docs/reference/print.draws_rvars.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/print.draws_summary.html b/docs/reference/print.draws_summary.html
index 9c830356..839f1c91 100644
--- a/docs/reference/print.draws_summary.html
+++ b/docs/reference/print.draws_summary.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/print.rvar.html b/docs/reference/print.rvar.html
index 67ece360..0f3029d0 100644
--- a/docs/reference/print.rvar.html
+++ b/docs/reference/print.rvar.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/ps_convergence_rate.html b/docs/reference/ps_convergence_rate.html
new file mode 100644
index 00000000..d31c252c
--- /dev/null
+++ b/docs/reference/ps_convergence_rate.html
@@ -0,0 +1,173 @@
+
+Pareto convergence rate — ps_convergence_rate • posterior
+
+
+
+
+
+
+
+
+
+
+
+
+ Given number of draws and scalar or array of k's, compute the
+relative convergence rate of PSIS estimate RMSE. See Appendix B of
+Vehtari et al. (2024). This function is exported to be usable by
+other packages. For user-facing diagnostic functions, see
+pareto_convergence_rate and pareto_diags.
+
+
+
+ ps_convergence_rate(k, ndraws, ...)
+
+
+
+ Arguments
+ - k
+pareto-k values
+
+
+- ndraws
+number of draws
+
+
+- ...
+unused
+
+
+
+ Value
+
+
+convergence rate
+
+
+ See also
+ Other helper-functions:
+ps_khat_threshold(),
+ps_min_ss(),
+ps_tail_length()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/docs/reference/ps_khat_threshold.html b/docs/reference/ps_khat_threshold.html
new file mode 100644
index 00000000..5e824367
--- /dev/null
+++ b/docs/reference/ps_khat_threshold.html
@@ -0,0 +1,171 @@
+
+Pareto k-hat threshold — ps_khat_threshold • posterior
+
+
+
+
+
+
+
+
+
+
+
+
+ Given number of draws, computes khat threshold for reliable Pareto
+smoothed estimate (to have small probability of large error). See
+section 3.2.4, equation (13) of Vehtari et al. (2024). This
+function is exported to be usable by other packages. For
+user-facing diagnostic functions, see pareto_khat_threshold and
+pareto_diags.
+
+
+
+ ps_khat_threshold(ndraws, ...)
+
+
+
+ Arguments
+ - ndraws
+number of draws
+
+
+- ...
+unused
+
+
+
+ Value
+
+
+threshold
+
+
+ See also
+ Other helper-functions:
+ps_convergence_rate(),
+ps_min_ss(),
+ps_tail_length()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/docs/reference/ps_min_ss.html b/docs/reference/ps_min_ss.html
new file mode 100644
index 00000000..bf028f3c
--- /dev/null
+++ b/docs/reference/ps_min_ss.html
@@ -0,0 +1,169 @@
+
+Pareto-smoothing minimum sample-size — ps_min_ss • posterior
+
+
+
+
+
+
+
+
+
+
+
+
+ Given Pareto-k computes the minimum sample size for reliable Pareto
+smoothed estimate (to have small probability of large error). See
+section 3.2.3 in Vehtari et al. (2024). This function is exported
+to be usable by other packages. For user-facing diagnostic functions, see
+pareto_min_ss and pareto_diags.
+
+
+
+ ps_min_ss(k, ...)
+
+
+
+ Arguments
+ - k
+pareto k value
+
+
+- ...
+unused
+
+
+
+ Value
+
+
+minimum sample size
+
+
+ See also
+ Other helper-functions:
+ps_convergence_rate(),
+ps_khat_threshold(),
+ps_tail_length()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/docs/reference/ps_tail_length.html b/docs/reference/ps_tail_length.html
new file mode 100644
index 00000000..514d8ad6
--- /dev/null
+++ b/docs/reference/ps_tail_length.html
@@ -0,0 +1,169 @@
+
+Pareto tail length — ps_tail_length • posterior
+
+
+
+
+
+
+
+
+
+
+
+
+ Calculate the tail length from number of draws and relative efficiency
+r_eff. See Appendix H in Vehtari et al. (2024). This function is
+used internally and is exported to be available for other packages.
+
+
+
+ ps_tail_length(ndraws, r_eff, ...)
+
+
+
+ Arguments
+ - ndraws
+number of draws
+
+
+- r_eff
+relative efficiency
+
+
+- ...
+unused
+
+
+
+ Value
+
+
+tail length
+
+
+ See also
+ Other helper-functions:
+ps_convergence_rate(),
+ps_khat_threshold(),
+ps_min_ss()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/docs/reference/quantile2.html b/docs/reference/quantile2.html
index 24a0721e..e2ce3067 100644
--- a/docs/reference/quantile2.html
+++ b/docs/reference/quantile2.html
@@ -18,7 +18,7 @@
diff --git a/docs/reference/r_scale.html b/docs/reference/r_scale.html
index 6aa29382..65b1e65f 100644
--- a/docs/reference/r_scale.html
+++ b/docs/reference/r_scale.html
@@ -19,7 +19,7 @@
diff --git a/docs/reference/rdo.html b/docs/reference/rdo.html
index 91851905..1ba8c754 100644
--- a/docs/reference/rdo.html
+++ b/docs/reference/rdo.html
@@ -18,7 +18,7 @@
@@ -117,7 +117,7 @@ Execute expressions of random variables
Arguments
- expr
(expression) A bare expression that can (optionally) contain
-rvars. The expression supports quasiquotation.
+rvars. The expression supports quasiquotation.
- dim
diff --git a/docs/reference/reexports.html b/docs/reference/reexports.html
index cd209a64..eb2e2209 100644
--- a/docs/reference/reexports.html
+++ b/docs/reference/reexports.html
@@ -24,7 +24,7 @@
diff --git a/docs/reference/rename_variables.html b/docs/reference/rename_variables.html
index 14fca7b3..fef58b4e 100644
--- a/docs/reference/rename_variables.html
+++ b/docs/reference/rename_variables.html
@@ -17,7 +17,7 @@
@@ -137,29 +137,29 @@ Value
Examples
x <- as_draws_df(example_draws())
-variables(x)
+variables(x)
#> [1] "mu" "tau" "theta[1]" "theta[2]" "theta[3]" "theta[4]"
#> [7] "theta[5]" "theta[6]" "theta[7]" "theta[8]"
x <- rename_variables(x, mean = mu, sigma = tau)
-variables(x)
+variables(x)
#> [1] "mean" "sigma" "theta[1]" "theta[2]" "theta[3]" "theta[4]"
#> [7] "theta[5]" "theta[6]" "theta[7]" "theta[8]"
x <- rename_variables(x, b = `theta[1]`) # or b = "theta[1]"
-variables(x)
+variables(x)
#> [1] "mean" "sigma" "b" "theta[2]" "theta[3]" "theta[4]"
#> [7] "theta[5]" "theta[6]" "theta[7]" "theta[8]"
# rename all elements of 'theta' at once
x <- rename_variables(x, alpha = theta)
-variables(x)
+variables(x)
#> [1] "mean" "sigma" "b" "alpha[2]" "alpha[3]" "alpha[4]"
#> [7] "alpha[5]" "alpha[6]" "alpha[7]" "alpha[8]"
diff --git a/docs/reference/repair_draws.html b/docs/reference/repair_draws.html
index 1ee2335a..2228cf9b 100644
--- a/docs/reference/repair_draws.html
+++ b/docs/reference/repair_draws.html
@@ -18,7 +18,7 @@
diff --git a/docs/reference/resample_draws.html b/docs/reference/resample_draws.html
index 79bc645f..c509774e 100644
--- a/docs/reference/resample_draws.html
+++ b/docs/reference/resample_draws.html
@@ -18,7 +18,7 @@
diff --git a/docs/reference/reserved_variables.html b/docs/reference/reserved_variables.html
index f0fb854a..5c5ec954 100644
--- a/docs/reference/reserved_variables.html
+++ b/docs/reference/reserved_variables.html
@@ -17,7 +17,7 @@
@@ -168,7 +168,6 @@ Examples
# if we add weights, the `.log_weight` reserved variable is used
x <- weight_draws(x, rexp(ndraws(x)))
-#> Loading required namespace: testthat
reserved_variables(x)
#> [1] ".log_weight"
diff --git a/docs/reference/rfun.html b/docs/reference/rfun.html
index dbe4c554..5489821f 100644
--- a/docs/reference/rfun.html
+++ b/docs/reference/rfun.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/rhat.html b/docs/reference/rhat.html
index e2c377c6..9ab5cd2a 100644
--- a/docs/reference/rhat.html
+++ b/docs/reference/rhat.html
@@ -19,7 +19,7 @@
@@ -174,6 +174,8 @@ See also
mcse_mean(),
mcse_quantile(),
mcse_sd(),
+pareto_diags(),
+pareto_khat(),
rhat_basic(),
rhat_nested(),
rstar()
diff --git a/docs/reference/rhat_basic.html b/docs/reference/rhat_basic.html
index 5b3e5330..42f61c06 100644
--- a/docs/reference/rhat_basic.html
+++ b/docs/reference/rhat_basic.html
@@ -20,7 +20,7 @@
@@ -184,8 +184,10 @@ See also
mcse_mean(),
mcse_quantile(),
mcse_sd(),
-rhat_nested(),
+pareto_diags(),
+pareto_khat(),
rhat(),
+rhat_nested(),
rstar()
diff --git a/docs/reference/rhat_nested.html b/docs/reference/rhat_nested.html
index 6561f222..327c5ff0 100644
--- a/docs/reference/rhat_nested.html
+++ b/docs/reference/rhat_nested.html
@@ -18,7 +18,7 @@
@@ -191,8 +191,10 @@ See also
mcse_mean(),
mcse_quantile(),
mcse_sd(),
-rhat_basic(),
+pareto_diags(),
+pareto_khat(),
rhat(),
+rhat_basic(),
rstar()
diff --git a/docs/reference/rstar.html b/docs/reference/rstar.html
index f64dab4d..a6a1957c 100644
--- a/docs/reference/rstar.html
+++ b/docs/reference/rstar.html
@@ -22,7 +22,7 @@
@@ -226,9 +226,11 @@ See also
mcse_mean(),
mcse_quantile(),
mcse_sd(),
+pareto_diags(),
+pareto_khat(),
+rhat(),
rhat_basic(),
-rhat_nested(),
-rhat()
+rhat_nested()
@@ -249,6 +251,7 @@ Examples
# can use other classification methods in caret library
print(rstar(x, method = "knn"))
}
+#> Warning: package ‘ggplot2’ was built under R version 4.2.3
#> randomForest 4.7-1.1
#> Type rfNews() to see new features/changes/bug fixes.
#>
diff --git a/docs/reference/rvar-dist.html b/docs/reference/rvar-dist.html
index 041060ea..bea313a8 100644
--- a/docs/reference/rvar-dist.html
+++ b/docs/reference/rvar-dist.html
@@ -19,7 +19,7 @@
diff --git a/docs/reference/rvar-matmult.html b/docs/reference/rvar-matmult.html
index eaed0eef..f14f8934 100644
--- a/docs/reference/rvar-matmult.html
+++ b/docs/reference/rvar-matmult.html
@@ -17,7 +17,7 @@
@@ -108,7 +108,10 @@ Matrix multiplication of random variables
- x %**% y
+ x %**% y
+
+# S3 method for rvar
+matrixOps(x, y)
@@ -140,8 +143,9 @@ Details
by rvars and are broadcasted across all draws of the rvar argument. Tensor multiplication
is used to efficiently multiply matrices across draws, so if either x or y is an rvar,
x %**% y will be much faster than rdo(x %*% y).
-Because rvar is an S3 class and S3 classes cannot properly override %*%, rvars use
-%**% for matrix multiplication.
+In R >= 4.3, you can also use %*% in place of %**% for matrix multiplication
+of rvars. In R < 4.3, S3 classes cannot properly override %*%, so
+you must use %**% for matrix multiplication of rvars.
diff --git a/docs/reference/rvar-slice.html b/docs/reference/rvar-slice.html
index c6d628fa..54bd7b12 100644
--- a/docs/reference/rvar-slice.html
+++ b/docs/reference/rvar-slice.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/rvar-summaries-over-draws.html b/docs/reference/rvar-summaries-over-draws.html
index 9e6802a1..da3cbcaa 100644
--- a/docs/reference/rvar-summaries-over-draws.html
+++ b/docs/reference/rvar-summaries-over-draws.html
@@ -19,7 +19,7 @@
diff --git a/docs/reference/rvar-summaries-within-draws.html b/docs/reference/rvar-summaries-within-draws.html
index cbc05372..0fc61b17 100644
--- a/docs/reference/rvar-summaries-within-draws.html
+++ b/docs/reference/rvar-summaries-within-draws.html
@@ -19,7 +19,7 @@
diff --git a/docs/reference/rvar.html b/docs/reference/rvar.html
index 4f350500..93f77141 100644
--- a/docs/reference/rvar.html
+++ b/docs/reference/rvar.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/rvar_apply.html b/docs/reference/rvar_apply.html
index bebd22b3..9e35dd06 100644
--- a/docs/reference/rvar_apply.html
+++ b/docs/reference/rvar_apply.html
@@ -19,7 +19,7 @@
diff --git a/docs/reference/rvar_factor.html b/docs/reference/rvar_factor.html
index f0226646..b76b6620 100644
--- a/docs/reference/rvar_factor.html
+++ b/docs/reference/rvar_factor.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/rvar_ifelse.html b/docs/reference/rvar_ifelse.html
index 90afcfa7..5d833dd2 100644
--- a/docs/reference/rvar_ifelse.html
+++ b/docs/reference/rvar_ifelse.html
@@ -17,7 +17,7 @@
diff --git a/docs/reference/rvar_is_finite.html b/docs/reference/rvar_is_finite.html
index 69560d36..296c67ea 100644
--- a/docs/reference/rvar_is_finite.html
+++ b/docs/reference/rvar_is_finite.html
@@ -18,7 +18,7 @@
diff --git a/docs/reference/rvar_rng.html b/docs/reference/rvar_rng.html
index 8b40e0fe..044e9f01 100644
--- a/docs/reference/rvar_rng.html
+++ b/docs/reference/rvar_rng.html
@@ -18,7 +18,7 @@
diff --git a/docs/reference/split_chains.html b/docs/reference/split_chains.html
index 2b1a475a..150f7759 100644
--- a/docs/reference/split_chains.html
+++ b/docs/reference/split_chains.html
@@ -18,7 +18,7 @@
diff --git a/docs/reference/sub-.draws_array.html b/docs/reference/sub-.draws_array.html
index 197198f7..c9f4c838 100644
--- a/docs/reference/sub-.draws_array.html
+++ b/docs/reference/sub-.draws_array.html
@@ -23,7 +23,7 @@
diff --git a/docs/reference/sub-.draws_matrix.html b/docs/reference/sub-.draws_matrix.html
index cd37d42a..314d6b44 100644
--- a/docs/reference/sub-.draws_matrix.html
+++ b/docs/reference/sub-.draws_matrix.html
@@ -23,7 +23,7 @@
diff --git a/docs/reference/subset_draws.html b/docs/reference/subset_draws.html
index 9c6594af..eb1fd8fd 100644
--- a/docs/reference/subset_draws.html
+++ b/docs/reference/subset_draws.html
@@ -17,7 +17,7 @@
@@ -119,6 +119,8 @@ Subset draws objects
draw = NULL,
regex = FALSE,
unique = TRUE,
+ exclude = FALSE,
+ scalar = FALSE,
...
)
@@ -131,6 +133,8 @@ Subset draws objects
draw = NULL,
regex = FALSE,
unique = TRUE,
+ exclude = FALSE,
+ scalar = FALSE,
...
)
@@ -143,6 +147,8 @@ Subset draws objects
draw = NULL,
regex = FALSE,
unique = TRUE,
+ exclude = FALSE,
+ scalar = FALSE,
...
)
@@ -155,6 +161,8 @@ Subset draws objects
draw = NULL,
regex = FALSE,
unique = TRUE,
+ exclude = FALSE,
+ scalar = FALSE,
...
)
@@ -167,6 +175,8 @@ Subset draws objects
draw = NULL,
regex = FALSE,
unique = TRUE,
+ exclude = FALSE,
+ scalar = FALSE,
...
)
@@ -189,8 +199,8 @@ Arguments
- variable
-(character vector) The variables to select. All elements of
-non-scalar variables can be selected at once.
+(character vector) The variables to select. All
+elements of non-scalar variables can be selected at once.
- iteration
@@ -202,21 +212,36 @@ Arguments
- draw
-(integer vector) The draw indices to be select. Subsetting draw
-indices will lead to an automatic merging of chains via merge_chains.
+(integer vector) The draw indices to be
+select. Subsetting draw indices will lead to an automatic merging
+of chains via merge_chains.
- regex
(logical) Should variable should be treated as a
-(vector of) regular expressions? Any variable in x matching at least one
-of the regular expressions will be selected. Defaults to FALSE.
+(vector of) regular expressions? Any variable in x matching at
+least one of the regular expressions will be selected. Defaults
+to FALSE.
- unique
-(logical) Should duplicated selection of chains, iterations, or
-draws be allowed? If TRUE (the default) only unique chains, iterations,
-and draws are selected regardless of how often they appear in the
-respective selecting arguments.
+(logical) Should duplicated selection of chains,
+iterations, or draws be allowed? If TRUE (the default) only
+unique chains, iterations, and draws are selected regardless of
+how often they appear in the respective selecting arguments.
+
+
+- exclude
+(logical) Should the selected subset be excluded?
+If FALSE (the default) only the selected subset will be
+returned. If TRUE everything but the selected subset will be
+returned.
+
+
+- scalar
+(logical) Should only scalar variables be selected?
+If FALSE (the default), all variables with matching names and
+arbitrary indices will be selected (see examples).
@@ -435,6 +460,9 @@ Examples
#>
#> # ... with 95 more iterations, and 4 more variables
+# trying to extract only a scalar 'theta' will fail
+# subset_draws(x, variable = "theta", scalar = TRUE)
+
diff --git a/docs/reference/thin_draws.html b/docs/reference/thin_draws.html
index afbc4697..16d871d3 100644
--- a/docs/reference/thin_draws.html
+++ b/docs/reference/thin_draws.html
@@ -1,5 +1,6 @@
-Thin draws objects — thin_draws • posterior Thin draws objects — thin_draws • posterior
@@ -17,7 +18,7 @@
@@ -104,17 +105,18 @@ Thin draws objects
- Thin draws objects to reduce their size and autocorrelation in the chains.
+ Thin draws objects to reduce their size and autocorrelation in
+the chains.
- thin_draws(x, thin, ...)
+ thin_draws(x, thin = NULL, ...)
# S3 method for draws
-thin_draws(x, thin, ...)
+thin_draws(x, thin = NULL, ...)
# S3 method for rvar
-thin_draws(x, thin, ...)
+thin_draws(x, thin = NULL, ...)
@@ -125,7 +127,15 @@ Arguments
thin
-(positive integer) The period for selecting draws.
+(positive numeric) The period for selecting draws. Must
+be between 1 and the number of iterations. If the value is not an
+integer, the draws will be selected such that the number of draws
+returned is equal to round(ndraws(x) / thin). Intervals between
+selected draws will be either ceiling(thin) or floor(thin), such
+that the average interval will be close to the thin value. If
+NULL, it will be automatically calculated based on bulk and
+tail effective sample size as suggested by Säilynoja et
+al. (2022).
...
@@ -138,6 +148,13 @@ Value
A draws object of the same class as x.
+
+ References
+ Teemu Säilynoja, Paul-Christian Bürkner, and Aki Vehtari (2022).
+Graphical test for discrete uniformity and its applications in
+goodness-of-fit evaluation and multiple sample comparison.
+Statistics and Computing. 32, 32. doi:10.1007/s11222-022-10090-6
+
Examples
diff --git a/docs/reference/u_scale.html b/docs/reference/u_scale.html
index e01e2af6..31194a50 100644
--- a/docs/reference/u_scale.html
+++ b/docs/reference/u_scale.html
@@ -20,7 +20,7 @@
diff --git a/docs/reference/variables-set.html b/docs/reference/variables-set.html
new file mode 100644
index 00000000..d410a547
--- /dev/null
+++ b/docs/reference/variables-set.html
@@ -0,0 +1,215 @@
+
+Set variable names in draws objects — variables<- • posterior
+
+
+
+
+
+
+
+
+
+
+
+
+ Set variable names for all variables in a draws object. The
+set_variables() form is useful when using pipe operators.
+
+
+
+ variables(x, ...) <- value
+
+# S3 method for draws_matrix
+variables(x, with_indices = TRUE, ...) <- value
+
+# S3 method for draws_array
+variables(x, with_indices = TRUE, ...) <- value
+
+# S3 method for draws_df
+variables(x, with_indices = TRUE, ...) <- value
+
+# S3 method for draws_list
+variables(x, with_indices = TRUE, ...) <- value
+
+# S3 method for draws_rvars
+variables(x, with_indices = FALSE, ...) <- value
+
+set_variables(x, variables, ...)
+
+
+
+ Arguments
+ - x
+(draws) A draws object or another R object for which the method
+is defined.
+
+
+- ...
+Arguments passed to individual methods (if applicable).
+
+
+- value, variables
+(character vector) new variable names.
+
+
+- with_indices
+(logical) Should indices be included in variable
+names? For example, if the object includes variables named "x[1]" and
+"x[2]", if TRUE, c("x[1]", "x[2]") is returned; if FALSE, only "x"
+is returned. Defaults to TRUE for all formats except draws_rvars().
+
+
+
+ Value
+
+
+Returns a draws object of the same format as x, with
+variables named as specified.
+
+
+ Details
+ variables(x) <- value allows you to modify the vector of variable names,
+similar to how names(x) <- value works for vectors and lists. For renaming
+specific variables, set_variables(x, value) works equivalently, but is more intuitive
+when using the pipe operator.
+For renaming specific variables, rename_variables() may offer a more
+convenient approach.
+
+
+ See also
+
+
+
+
+ Examples
+ x <- example_draws()
+
+variables(x)
+#> [1] "mu" "tau" "theta[1]" "theta[2]" "theta[3]" "theta[4]"
+#> [7] "theta[5]" "theta[6]" "theta[7]" "theta[8]"
+nvariables(x)
+#> [1] 10
+variables(x) <- letters[1:nvariables(x)]
+
+# or equivalently...
+x <- set_variables(x, letters[1:nvariables(x)])
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/docs/reference/variables.html b/docs/reference/variables.html
new file mode 100644
index 00000000..6705a7e1
--- /dev/null
+++ b/docs/reference/variables.html
@@ -0,0 +1,209 @@
+
+Get variable names from draws objects — variables • posterior
+
+
+
+
+
+
+
+
+
+
+
+
+ Get variable names from draws objects.
+
+
+
+ variables(x, ...)
+
+# S3 method for draws_matrix
+variables(x, reserved = FALSE, with_indices = TRUE, ...)
+
+# S3 method for draws_array
+variables(x, reserved = FALSE, with_indices = TRUE, ...)
+
+# S3 method for draws_df
+variables(x, reserved = FALSE, with_indices = TRUE, ...)
+
+# S3 method for draws_list
+variables(x, reserved = FALSE, with_indices = TRUE, ...)
+
+# S3 method for draws_rvars
+variables(x, reserved = FALSE, with_indices = FALSE, ...)
+
+nvariables(x, ...)
+
+
+
+ Arguments
+ - x
+(draws) A draws object or another R object for which the method
+is defined.
+
+
+- ...
+Arguments passed to individual methods (if applicable).
+
+
+- reserved
+(logical) Should reserved variables be included in the
+output? Defaults to FALSE. See reserved_variables for an overview of
+currently reserved variable names.
+
+
+- with_indices
+(logical) Should indices be included in variable
+names? For example, if the object includes variables named "x[1]" and
+"x[2]", if TRUE, c("x[1]", "x[2]") is returned; if FALSE, only "x"
+is returned. Defaults to TRUE for all formats except draws_rvars().
+
+
+
+ Value
+
+
+For variables(), a character vector.
+
+
+For nvariables(), a scalar integer.
+
+
+ Details
+ variables() returns a vector of all variable names, and nvariables()
+returns the number of variables.
+
+
+ See also
+ variables<-, rename_variables, draws-index
+
+
+
+ Examples
+ x <- example_draws()
+
+variables(x)
+#> [1] "mu" "tau" "theta[1]" "theta[2]" "theta[3]" "theta[4]"
+#> [7] "theta[5]" "theta[6]" "theta[7]" "theta[8]"
+nvariables(x)
+#> [1] 10
+variables(x) <- letters[1:nvariables(x)]
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/docs/reference/weight_draws.html b/docs/reference/weight_draws.html
index eac63e84..6eb3931d 100644
--- a/docs/reference/weight_draws.html
+++ b/docs/reference/weight_draws.html
@@ -21,7 +21,7 @@
@@ -119,19 +119,19 @@ Weight draws objects
weight_draws(x, weights, ...)
# S3 method for draws_matrix
-weight_draws(x, weights, log = FALSE, ...)
+weight_draws(x, weights, log = FALSE, pareto_smooth = FALSE, ...)
# S3 method for draws_array
-weight_draws(x, weights, log = FALSE, ...)
+weight_draws(x, weights, log = FALSE, pareto_smooth = FALSE, ...)
# S3 method for draws_df
-weight_draws(x, weights, log = FALSE, ...)
+weight_draws(x, weights, log = FALSE, pareto_smooth = FALSE, ...)
# S3 method for draws_list
-weight_draws(x, weights, log = FALSE, ...)
+weight_draws(x, weights, log = FALSE, pareto_smooth = FALSE, ...)
# S3 method for draws_rvars
-weight_draws(x, weights, log = FALSE, ...)
+weight_draws(x, weights, log = FALSE, pareto_smooth = FALSE, ...)
@@ -153,10 +153,15 @@ Arguments
log
-(logicla) Are the weights passed already on the log scale? The
+
(logical) Are the weights passed already on the log scale? The
default is FALSE, that is, expecting weights to be on the standard
(non-log) scale.
+
+pareto_smooth
+(logical) Should the weights be Pareto-smoothed?
+The default is FALSE.
+
Value
@@ -176,32 +181,36 @@ Examples
# sample some random weights for illustration
wts <- rexp(ndraws(x))
head(wts)
-#> [1] 0.2644506 1.4629054 1.9719158 0.8050840 2.4207868 1.4064785
+#> [1] 0.07943442 1.20975078 0.09003050 0.22634413 0.47107425 0.22054226
# add weights
x <- weight_draws(x, weights = wts)
# extract weights
head(weights(x)) # defaults to normalized weights
-#> [1] 0.0007016715 0.0038815529 0.0052321195 0.0021361438 0.0064231168
-#> [6] 0.0037318345
+#> [1] 0.0002133569 0.0032493306 0.0002418175 0.0006079491 0.0012652821
+#> [6] 0.0005923656
head(weights(x, normalize=FALSE)) # recover original weights
-#> [1] 0.2644506 1.4629054 1.9719158 0.8050840 2.4207868 1.4064785
+#> [1] 0.07943442 1.20975078 0.09003050 0.22634413 0.47107425 0.22054226
head(weights(x, log=TRUE)) # get normalized log-weights
-#> [1] -7.262045 -5.551520 -5.252939 -6.148753 -5.047852 -5.590855
+#> [1] -8.452544 -5.729306 -8.327327 -7.405419 -6.672460 -7.431387
# add weights which are already on the log scale
log_wts <- log(wts)
head(log_wts)
-#> [1] -1.3301008 0.3804244 0.6790056 -0.2168086 0.8840926 0.3410891
+#> [1] -2.5328235 0.1904144 -2.4076068 -1.4856988 -0.7527396 -1.5116659
x <- weight_draws(x, weights = log_wts, log = TRUE)
# extract weights
head(weights(x))
-#> [1] 0.0007016715 0.0038815529 0.0052321195 0.0021361438 0.0064231168
-#> [6] 0.0037318345
+#> [1] 0.0002133569 0.0032493306 0.0002418175 0.0006079491 0.0012652821
+#> [6] 0.0005923656
head(weights(x, log=TRUE, normalize = FALSE)) # recover original log_wts
-#> [1] -1.3301008 0.3804244 0.6790056 -0.2168086 0.8840926 0.3410891
+#> [1] -2.5328235 0.1904144 -2.4076068 -1.4856988 -0.7527396 -1.5116659
+
+# add weights on log scale and Pareto smooth them
+x <- weight_draws(x, weights = log_wts, log = TRUE, pareto_smooth = TRUE)
+#> Pareto k-hat = 0.04.
diff --git a/docs/reference/weights.draws.html b/docs/reference/weights.draws.html
index 51f94eb2..6e4ced68 100644
--- a/docs/reference/weights.draws.html
+++ b/docs/reference/weights.draws.html
@@ -18,7 +18,7 @@
@@ -152,32 +152,36 @@ Examples
# sample some random weights for illustration
wts <- rexp(ndraws(x))
head(wts)
-#> [1] 0.61086446 0.33461334 0.02989975 0.36790580 2.15576680 0.24003173
+#> [1] 0.44924386 0.01132425 0.79627288 0.06250828 0.26065786 0.04660613
# add weights
x <- weight_draws(x, weights = wts)
# extract weights
head(weights(x)) # defaults to normalized weights
-#> [1] 1.574879e-03 8.626719e-04 7.708502e-05 9.485037e-04 5.557816e-03
-#> [6] 6.188295e-04
+#> [1] 1.191764e-03 3.004123e-05 2.112371e-03 1.658234e-04 6.914792e-04
+#> [6] 1.236378e-04
head(weights(x, normalize=FALSE)) # recover original weights
-#> [1] 0.61086446 0.33461334 0.02989975 0.36790580 2.15576680 0.24003173
+#> [1] 0.44924386 0.01132425 0.79627288 0.06250828 0.26065786 0.04660613
head(weights(x, log=TRUE)) # get normalized log-weights
-#> [1] -6.453577 -7.055476 -9.470602 -6.960625 -5.192550 -7.387681
+#> [1] -6.732320 -10.412940 -6.159944 -8.704587 -7.276678 -8.998154
# add weights which are already on the log scale
log_wts <- log(wts)
head(log_wts)
-#> [1] -0.4928802 -1.0947796 -3.5099051 -0.9999283 0.7681465 -1.4269842
+#> [1] -0.8001894 -4.4808090 -0.2278133 -2.7724563 -1.3445466 -3.0660231
x <- weight_draws(x, weights = log_wts, log = TRUE)
# extract weights
head(weights(x))
-#> [1] 1.574879e-03 8.626719e-04 7.708502e-05 9.485037e-04 5.557816e-03
-#> [6] 6.188295e-04
+#> [1] 1.191764e-03 3.004123e-05 2.112371e-03 1.658234e-04 6.914792e-04
+#> [6] 1.236378e-04
head(weights(x, log=TRUE, normalize = FALSE)) # recover original log_wts
-#> [1] -0.4928802 -1.0947796 -3.5099051 -0.9999283 0.7681465 -1.4269842
+#> [1] -0.8001894 -4.4808090 -0.2278133 -2.7724563 -1.3445466 -3.0660231
+
+# add weights on log scale and Pareto smooth them
+x <- weight_draws(x, weights = log_wts, log = TRUE, pareto_smooth = TRUE)
+#> Pareto k-hat = 0.08.
diff --git a/docs/reference/z_scale.html b/docs/reference/z_scale.html
index e2f83d8e..bbb6b5e9 100644
--- a/docs/reference/z_scale.html
+++ b/docs/reference/z_scale.html
@@ -20,7 +20,7 @@
diff --git a/docs/sitemap.xml b/docs/sitemap.xml
index 6546f2e1..7ff5a983 100644
--- a/docs/sitemap.xml
+++ b/docs/sitemap.xml
@@ -3,6 +3,9 @@
https://mc-stan.org/posterior/404.html
+
+ https://mc-stan.org/posterior/CONTRIBUTING.html
+
https://mc-stan.org/posterior/LICENSE-text.html
@@ -12,6 +15,9 @@
https://mc-stan.org/posterior/articles/index.html
+
+ https://mc-stan.org/posterior/articles/pareto_diagnostics.html
+
https://mc-stan.org/posterior/articles/posterior.html
@@ -120,6 +126,9 @@
https://mc-stan.org/posterior/reference/extract_variable.html
+
+ https://mc-stan.org/posterior/reference/extract_variable_array.html
+
https://mc-stan.org/posterior/reference/extract_variable_matrix.html
@@ -195,6 +204,18 @@
https://mc-stan.org/posterior/reference/print.rvar.html
+
+ https://mc-stan.org/posterior/reference/ps_convergence_rate.html
+
+
+ https://mc-stan.org/posterior/reference/ps_khat_threshold.html
+
+
+ https://mc-stan.org/posterior/reference/ps_min_ss.html
+
+
+ https://mc-stan.org/posterior/reference/ps_tail_length.html
+
https://mc-stan.org/posterior/reference/quantile2.html
@@ -300,6 +321,12 @@
https://mc-stan.org/posterior/reference/u_scale.html
+
+ https://mc-stan.org/posterior/reference/variables-set.html
+
+
+ https://mc-stan.org/posterior/reference/variables.html
+
https://mc-stan.org/posterior/reference/vctrs-compat.html
diff --git a/inst/CITATION b/inst/CITATION
index 8f9f75dc..f92fc2ae 100644
--- a/inst/CITATION
+++ b/inst/CITATION
@@ -21,7 +21,64 @@ for assessing convergence of MCMC (with discussion)",
person("Daniel", "Simpson"),
person("Bob", "Carpenter"),
person("Paul-Christian", "Bürkner")),
- year = "2021",
journal = "Bayesian Analysis",
- header = "To cite the MCMC convergence diagnostics:"
+ year = "2021",
+ volume = "16",
+ number = "2",
+ pages = "667-718",
+ header = "To cite the MCMC convergence diagnostics (`rhat`, `ess_bulk`, `ess_tail`,
+ `ess_median`, `ess_quantile`, `mcse_median`, and `mcse_quantile`):"
+)
+
+bibentry(
+ title = "Nested Rhat: Assessing the convergence of Markov chain Monte Carlo when running many short chains",
+ bibtype = "Article",
+ author = c(
+ person("Charles C.", "Margossian"),
+ person("Matthew D.", "Hoffman"),
+ person("Pavel", "Sountsov"),
+ person("Lionel", "Riou-Durand"),
+ person("Aki", "Vehtari"),
+ person("Andrew", "Gelman")
+ ),
+ journal = "Bayesian Analysis",
+ year = 2024,
+ doi = "10.1214/24-BA1453",
+ header = "To cite MCMC convergence diagnostic `nested_rhat`:"
+)
+
+bibentry(
+ title = "Rstar: A robust MCMC convergence diagnostic with uncertainty using decision tree classifiers",
+ bibtype = "Article",
+ author = c(
+ person("Ben", "Lambert"),
+ person("Aki", "Vehtari")
+ ),
+ journal = "Bayesian Analysis",
+ year = 2022,
+ volume = 17,
+ number = 2,
+ pages = "353-379",
+ doi = "10.1214/20-BA1252",
+ header = "To cite MCMC convergence diagnostic `rstar`:"
+)
+
+bibentry(
+ title = "Pareto smoothed importance sampling",
+ bibtype = "Article",
+ author = c(
+ person("Aki", "Vehtari"),
+ person("Daniel", "Simpson"),
+ person("Andrew", "Gelman"),
+ person("Yuling", "Yao"),
+ person("Jonah", "Gabry")
+ ),
+ journal = "Journal of Machine Learning Research",
+ year = 2024,
+ volume = 25,
+ number = 72,
+ pages = "1-58",
+ header = "To cite Pareto-k diagnostics and Pareto smoothing (`pareto_khat`, `pareto_min_ss`,
+`pareto_convergence_rate`, `khat_threshold`, `pareto_diags`, and
+`pareto_smooth`):"
)
diff --git a/man-roxygen/args-pareto.R b/man-roxygen/args-pareto.R
index 8a4d92b4..488e60bf 100644
--- a/man-roxygen/args-pareto.R
+++ b/man-roxygen/args-pareto.R
@@ -12,10 +12,10 @@
#' `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)).
+#' al. (2024)).
#' @param r_eff (numeric) relative effective sample size estimate. If
#' `r_eff` is NULL, it will be calculated assuming the draws are
-#' from MCMC. Default is 1.
+#' from MCMC. Default is NULL.
#' @param verbose (logical) Should diagnostic messages be printed? If
#' `TRUE`, messages related to Pareto diagnostics will be
#' printed. Default is `FALSE`.
diff --git a/man-roxygen/ref-margossian-nestedrhat-2023.R b/man-roxygen/ref-margossian-nestedrhat-2024.R
similarity index 55%
rename from man-roxygen/ref-margossian-nestedrhat-2023.R
rename to man-roxygen/ref-margossian-nestedrhat-2024.R
index 474881e2..da9d00b2 100644
--- a/man-roxygen/ref-margossian-nestedrhat-2023.R
+++ b/man-roxygen/ref-margossian-nestedrhat-2024.R
@@ -1,5 +1,5 @@
#' @references
#' Charles C. Margossian, Matthew D. Hoffman, Pavel Sountsov, Lionel
-#' Riou-Durand, Aki Vehtari and Andrew Gelman (2023). Nested R-hat:
+#' Riou-Durand, Aki Vehtari and Andrew Gelman (2024). Nested R-hat:
#' Assessing the convergence of Markov chain Monte Carlo when running
-#' many short chains. arxiv:arXiv:2110.13017 (version 4)
+#' many short chains. *Bayesian Analysis*. doi:10.1214/24-BA1453
diff --git a/man-roxygen/ref-vehtari-paretosmooth-2022.R b/man-roxygen/ref-vehtari-paretosmooth-2022.R
index 30f526fc..addcf4c9 100644
--- a/man-roxygen/ref-vehtari-paretosmooth-2022.R
+++ b/man-roxygen/ref-vehtari-paretosmooth-2022.R
@@ -1,4 +1,5 @@
#' @references
#' Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao and
-#' Jonah Gabry (2022). Pareto Smoothed Importance Sampling.
-#' arxiv:arXiv:1507.02646 (version 8)
+#' Jonah Gabry (2024). Pareto Smoothed Importance Sampling.
+#' *Journal of Machine Learning Research*, 25(72):1-58.
+#' [PDF](https://jmlr.org/papers/v25/19-556.html)
diff --git a/man-roxygen/ref-vehtari-rhat-2021.R b/man-roxygen/ref-vehtari-rhat-2021.R
index 91ffea91..baaedd17 100644
--- a/man-roxygen/ref-vehtari-rhat-2021.R
+++ b/man-roxygen/ref-vehtari-rhat-2021.R
@@ -2,5 +2,5 @@
#' Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and
#' Paul-Christian Bürkner (2021). Rank-normalization, folding, and
#' localization: An improved R-hat for assessing convergence of
-#' MCMC (with discussion). *Bayesian Data Analysis*. 16(2), 667-–718.
+#' MCMC (with discussion). *Bayesian Analysis*. 16(2), 667-–718.
#' doi:10.1214/20-BA1221
diff --git a/man/diagnostics.Rd b/man/diagnostics.Rd
index f2ba4998..64b5b8fc 100644
--- a/man/diagnostics.Rd
+++ b/man/diagnostics.Rd
@@ -16,11 +16,14 @@ A list of available diagnostics and links to their individual help pages.
\code{\link[=ess_basic]{ess_basic()}} \tab Basic version of effective sample size \cr
\code{\link[=ess_bulk]{ess_bulk()}} \tab Bulk effective sample size \cr
\code{\link[=ess_tail]{ess_tail()}} \tab Tail effective sample size \cr
+ \code{\link[=ess_mean]{ess_mean()}} \tab Effective sample sizes for the mean \cr
+ \code{\link[=ess_median]{ess_median()}} \tab Effective sample sizes for the median \cr
\code{\link[=ess_quantile]{ess_quantile()}} \tab Effective sample sizes for quantiles \cr
- \code{\link[=ess_sd]{ess_sd()}} \tab Effective sample sizes for standard deviations \cr
+ \code{\link[=ess_sd]{ess_sd()}} \tab Effective sample sizes for the standard deviation \cr
\code{\link[=mcse_mean]{mcse_mean()}} \tab Monte Carlo standard error for the mean \cr
+ \code{\link[=mcse_median]{mcse_median()}} \tab Monte Carlo standard error for the median \cr
\code{\link[=mcse_quantile]{mcse_quantile()}} \tab Monte Carlo standard error for quantiles \cr
- \code{\link[=mcse_sd]{mcse_sd()}} \tab Monte Carlo standard error for standard deviations \cr
+ \code{\link[=mcse_sd]{mcse_sd()}} \tab Monte Carlo standard error for the standard deviation \cr
\code{\link[=pareto_khat]{pareto_khat()}} \tab Pareto khat diagnostic for tail(s) \cr
\code{\link[=pareto_diags]{pareto_diags()}} \tab Additional diagnostics related to Pareto khat \cr
\code{\link[=rhat_basic]{rhat_basic()}} \tab Basic version of Rhat \cr
diff --git a/man/draws_array.Rd b/man/draws_array.Rd
index 359703ca..941cf85e 100755
--- a/man/draws_array.Rd
+++ b/man/draws_array.Rd
@@ -74,10 +74,10 @@ str(x2)
}
\seealso{
Other formats:
+\code{\link{draws}},
\code{\link{draws_df}()},
\code{\link{draws_list}()},
\code{\link{draws_matrix}()},
-\code{\link{draws_rvars}()},
-\code{\link{draws}}
+\code{\link{draws_rvars}()}
}
\concept{formats}
diff --git a/man/draws_df.Rd b/man/draws_df.Rd
index a34daffb..bbfc8206 100755
--- a/man/draws_df.Rd
+++ b/man/draws_df.Rd
@@ -96,10 +96,10 @@ print(xnew)
}
\seealso{
Other formats:
+\code{\link{draws}},
\code{\link{draws_array}()},
\code{\link{draws_list}()},
\code{\link{draws_matrix}()},
-\code{\link{draws_rvars}()},
-\code{\link{draws}}
+\code{\link{draws_rvars}()}
}
\concept{formats}
diff --git a/man/draws_list.Rd b/man/draws_list.Rd
index f45526a8..b696350e 100755
--- a/man/draws_list.Rd
+++ b/man/draws_list.Rd
@@ -76,10 +76,10 @@ str(x2)
}
\seealso{
Other formats:
+\code{\link{draws}},
\code{\link{draws_array}()},
\code{\link{draws_df}()},
\code{\link{draws_matrix}()},
-\code{\link{draws_rvars}()},
-\code{\link{draws}}
+\code{\link{draws_rvars}()}
}
\concept{formats}
diff --git a/man/draws_matrix.Rd b/man/draws_matrix.Rd
index 1b548412..432e3f9b 100755
--- a/man/draws_matrix.Rd
+++ b/man/draws_matrix.Rd
@@ -74,10 +74,10 @@ str(x2)
}
\seealso{
Other formats:
+\code{\link{draws}},
\code{\link{draws_array}()},
\code{\link{draws_df}()},
\code{\link{draws_list}()},
-\code{\link{draws_rvars}()},
-\code{\link{draws}}
+\code{\link{draws_rvars}()}
}
\concept{formats}
diff --git a/man/draws_rvars.Rd b/man/draws_rvars.Rd
index 0e4b614a..28786629 100755
--- a/man/draws_rvars.Rd
+++ b/man/draws_rvars.Rd
@@ -77,10 +77,10 @@ str(x2)
}
\seealso{
Other formats:
+\code{\link{draws}},
\code{\link{draws_array}()},
\code{\link{draws_df}()},
\code{\link{draws_list}()},
-\code{\link{draws_matrix}()},
-\code{\link{draws}}
+\code{\link{draws_matrix}()}
}
\concept{formats}
diff --git a/man/ess_basic.Rd b/man/ess_basic.Rd
index 867076ca..86208d10 100755
--- a/man/ess_basic.Rd
+++ b/man/ess_basic.Rd
@@ -8,7 +8,7 @@
\usage{
ess_basic(x, ...)
-\method{ess_basic}{default}(x, split = TRUE, ...)
+\method{ess_basic}{default}(x, split = TRUE, weights = NULL, ...)
\method{ess_basic}{rvar}(x, split = TRUE, ...)
}
@@ -24,6 +24,9 @@ ess_basic(x, ...)
\item{split}{(logical) Should the estimate be computed on split chains? The
default is \code{TRUE}.}
+
+\item{weights}{(numeric vector) A vector of weights of length \code{ndraws(x)},
+or \code{NULL} to not have weights. Note that if \code{x} is an \code{rvar} with weights, they are used instead of this argument.}
}
\value{
If the input is an array, returns a single numeric value. If any of the draws
@@ -47,6 +50,8 @@ al. (2021). For practical applications, we strongly
recommend the improved ESS convergence diagnostics implemented in
\code{\link[=ess_bulk]{ess_bulk()}} and \code{\link[=ess_tail]{ess_tail()}}. See Vehtari (2021) for an in-depth
comparison of different effective sample size estimators.
+If computed on a weighted \code{rvar}, weights will be
+taken into account.
}
\examples{
mu <- extract_variable_matrix(example_draws(), "mu")
@@ -64,7 +69,7 @@ Hall/CRC.
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and
Paul-Christian Bürkner (2021). Rank-normalization, folding, and
localization: An improved R-hat for assessing convergence of
-MCMC (with discussion). \emph{Bayesian Data Analysis}. 16(2), 667-–718.
+MCMC (with discussion). \emph{Bayesian Analysis}. 16(2), 667-–718.
doi:10.1214/20-BA1221
Aki Vehtari (2021). Comparison of MCMC effective sample size estimators.
@@ -81,9 +86,9 @@ Other diagnostics:
\code{\link{mcse_sd}()},
\code{\link{pareto_diags}()},
\code{\link{pareto_khat}()},
+\code{\link{rhat}()},
\code{\link{rhat_basic}()},
\code{\link{rhat_nested}()},
-\code{\link{rhat}()},
\code{\link{rstar}()}
}
\concept{diagnostics}
diff --git a/man/ess_bulk.Rd b/man/ess_bulk.Rd
index c1456be3..74f62cbf 100755
--- a/man/ess_bulk.Rd
+++ b/man/ess_bulk.Rd
@@ -8,7 +8,7 @@
\usage{
ess_bulk(x, ...)
-\method{ess_bulk}{default}(x, ...)
+\method{ess_bulk}{default}(x, weights = NULL, ...)
\method{ess_bulk}{rvar}(x, ...)
}
@@ -21,6 +21,9 @@ ess_bulk(x, ...)
}}
\item{...}{Arguments passed to individual methods (if applicable).}
+
+\item{weights}{(numeric vector) A vector of weights of length \code{ndraws(x)},
+or \code{NULL} to not have weights. Note that if \code{x} is an \code{rvar} with weights, they are used instead of this argument.}
}
\value{
If the input is an array, returns a single numeric value. If any of the draws
@@ -44,6 +47,8 @@ the bulk of the posterior. It is defined as the effective sample size for
rank normalized values using split chains. For the tail effective sample size
see \code{\link[=ess_tail]{ess_tail()}}. See Vehtari (2021) for an in-depth
comparison of different effective sample size estimators.
+If computed on a weighted \code{rvar}, weights will be
+taken into account.
}
\examples{
mu <- extract_variable_matrix(example_draws(), "mu")
@@ -57,7 +62,7 @@ ess_bulk(d$Sigma)
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and
Paul-Christian Bürkner (2021). Rank-normalization, folding, and
localization: An improved R-hat for assessing convergence of
-MCMC (with discussion). \emph{Bayesian Data Analysis}. 16(2), 667-–718.
+MCMC (with discussion). \emph{Bayesian Analysis}. 16(2), 667-–718.
doi:10.1214/20-BA1221
Aki Vehtari (2021). Comparison of MCMC effective sample size estimators.
@@ -74,9 +79,9 @@ Other diagnostics:
\code{\link{mcse_sd}()},
\code{\link{pareto_diags}()},
\code{\link{pareto_khat}()},
+\code{\link{rhat}()},
\code{\link{rhat_basic}()},
\code{\link{rhat_nested}()},
-\code{\link{rhat}()},
\code{\link{rstar}()}
}
\concept{diagnostics}
diff --git a/man/ess_mean.Rd b/man/ess_mean.Rd
index cd157a8a..b595fa54 100755
--- a/man/ess_mean.Rd
+++ b/man/ess_mean.Rd
@@ -2,11 +2,14 @@
% Please edit documentation in R/convergence.R
\name{ess_mean}
\alias{ess_mean}
+\alias{ess_mean.default}
\alias{ess_mean.rvar}
\title{Effective sample size for the mean}
\usage{
ess_mean(x, ...)
+\method{ess_mean}{default}(x, weights = NULL, ...)
+
\method{ess_mean}{rvar}(x, ...)
}
\arguments{
@@ -18,6 +21,9 @@ ess_mean(x, ...)
}}
\item{...}{Arguments passed to individual methods (if applicable).}
+
+\item{weights}{(numeric vector) A vector of weights of length \code{ndraws(x)},
+or \code{NULL} to not have weights. Note that if \code{x} is an \code{rvar} with weights, they are used instead of this argument.}
}
\value{
If the input is an array, returns a single numeric value. If any of the draws
diff --git a/man/ess_quantile.Rd b/man/ess_quantile.Rd
index aa85c909..d58ae7a6 100755
--- a/man/ess_quantile.Rd
+++ b/man/ess_quantile.Rd
@@ -5,18 +5,15 @@
\alias{ess_quantile.default}
\alias{ess_quantile.rvar}
\alias{ess_median}
-\alias{ess_mean.default}
\title{Effective sample sizes for quantiles}
\usage{
ess_quantile(x, probs = c(0.05, 0.95), ...)
-\method{ess_quantile}{default}(x, probs = c(0.05, 0.95), names = TRUE, ...)
+\method{ess_quantile}{default}(x, probs = c(0.05, 0.95), names = TRUE, weights = NULL, ...)
\method{ess_quantile}{rvar}(x, probs = c(0.05, 0.95), names = TRUE, ...)
ess_median(x, ...)
-
-\method{ess_mean}{default}(x, ...)
}
\arguments{
\item{x}{(multiple options) One of:
@@ -33,6 +30,9 @@ ess_median(x, ...)
\item{names}{(logical) Should the result have a \code{names} attribute? The
default is \code{TRUE}, but use \code{FALSE} for improved speed if there are many
values in \code{probs}.}
+
+\item{weights}{(numeric vector) A vector of weights of length \code{ndraws(x)},
+or \code{NULL} to not have weights. Note that if \code{x} is an \code{rvar} with weights, they are used instead of this argument.}
}
\value{
If the input is an array,
@@ -54,8 +54,9 @@ result indexes the input probabilities; i.e. the result has dimension
\code{c(length(probs), dim(x))}.
}
\description{
-Compute effective sample size estimates for quantile estimates of a single
-variable.
+Compute effective sample size estimates for quantile estimates of a
+single variable. If computed on a weighted \code{rvar}, weights will be
+taken into account.
}
\examples{
mu <- extract_variable_matrix(example_draws(), "mu")
@@ -69,7 +70,7 @@ ess_quantile(d$mu, probs = c(0.1, 0.9))
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and
Paul-Christian Bürkner (2021). Rank-normalization, folding, and
localization: An improved R-hat for assessing convergence of
-MCMC (with discussion). \emph{Bayesian Data Analysis}. 16(2), 667-–718.
+MCMC (with discussion). \emph{Bayesian Analysis}. 16(2), 667-–718.
doi:10.1214/20-BA1221
}
\seealso{
@@ -83,9 +84,9 @@ Other diagnostics:
\code{\link{mcse_sd}()},
\code{\link{pareto_diags}()},
\code{\link{pareto_khat}()},
+\code{\link{rhat}()},
\code{\link{rhat_basic}()},
\code{\link{rhat_nested}()},
-\code{\link{rhat}()},
\code{\link{rstar}()}
}
\concept{diagnostics}
diff --git a/man/ess_sd.Rd b/man/ess_sd.Rd
index 38475d2a..f9bf0d26 100755
--- a/man/ess_sd.Rd
+++ b/man/ess_sd.Rd
@@ -8,7 +8,7 @@
\usage{
ess_sd(x, ...)
-\method{ess_sd}{default}(x, ...)
+\method{ess_sd}{default}(x, weights = NULL, ...)
\method{ess_sd}{rvar}(x, ...)
}
@@ -21,6 +21,9 @@ ess_sd(x, ...)
}}
\item{...}{Arguments passed to individual methods (if applicable).}
+
+\item{weights}{(numeric vector) A vector of weights of length \code{ndraws(x)},
+or \code{NULL} to not have weights. Note that if \code{x} is an \code{rvar} with weights, they are used instead of this argument.}
}
\value{
If the input is an array, returns a single numeric value. If any of the draws
@@ -41,6 +44,8 @@ passing the draws array for that element of the \code{\link{rvar}} to this funct
Compute an effective sample size estimate for the standard deviation (SD)
estimate of a single variable. This is defined as the effective sample size
estimate for the absolute deviation from mean.
+If computed on a weighted \code{rvar}, weights will be
+taken into account.
}
\examples{
mu <- extract_variable_matrix(example_draws(), "mu")
@@ -54,7 +59,7 @@ ess_sd(d$Sigma)
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and
Paul-Christian Bürkner (2021). Rank-normalization, folding, and
localization: An improved R-hat for assessing convergence of
-MCMC (with discussion). \emph{Bayesian Data Analysis}. 16(2), 667-–718.
+MCMC (with discussion). \emph{Bayesian Analysis}. 16(2), 667-–718.
doi:10.1214/20-BA1221
}
\seealso{
@@ -68,9 +73,9 @@ Other diagnostics:
\code{\link{mcse_sd}()},
\code{\link{pareto_diags}()},
\code{\link{pareto_khat}()},
+\code{\link{rhat}()},
\code{\link{rhat_basic}()},
\code{\link{rhat_nested}()},
-\code{\link{rhat}()},
\code{\link{rstar}()}
}
\concept{diagnostics}
diff --git a/man/ess_tail.Rd b/man/ess_tail.Rd
index 8f959718..f20f7ecb 100755
--- a/man/ess_tail.Rd
+++ b/man/ess_tail.Rd
@@ -8,7 +8,7 @@
\usage{
ess_tail(x, ...)
-\method{ess_tail}{default}(x, ...)
+\method{ess_tail}{default}(x, weights = NULL, ...)
\method{ess_tail}{rvar}(x, ...)
}
@@ -21,6 +21,9 @@ ess_tail(x, ...)
}}
\item{...}{Arguments passed to individual methods (if applicable).}
+
+\item{weights}{(numeric vector) A vector of weights of length \code{ndraws(x)},
+or \code{NULL} to not have weights. Note that if \code{x} is an \code{rvar} with weights, they are used instead of this argument.}
}
\value{
If the input is an array, returns a single numeric value. If any of the draws
@@ -44,6 +47,8 @@ the tails of the posterior. It is defined as the minimum of the effective
sample sizes for 5\% and 95\% quantiles. For the bulk effective sample
size see \code{\link[=ess_bulk]{ess_bulk()}}. See Vehtari (2021) for an in-depth
comparison of different effective sample size estimators.
+If computed on a weighted \code{rvar}, weights will be
+taken into account.
}
\examples{
mu <- extract_variable_matrix(example_draws(), "mu")
@@ -57,7 +62,7 @@ ess_tail(d$Sigma)
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and
Paul-Christian Bürkner (2021). Rank-normalization, folding, and
localization: An improved R-hat for assessing convergence of
-MCMC (with discussion). \emph{Bayesian Data Analysis}. 16(2), 667-–718.
+MCMC (with discussion). \emph{Bayesian Analysis}. 16(2), 667-–718.
doi:10.1214/20-BA1221
Aki Vehtari (2021). Comparison of MCMC effective sample size estimators.
@@ -74,9 +79,9 @@ Other diagnostics:
\code{\link{mcse_sd}()},
\code{\link{pareto_diags}()},
\code{\link{pareto_khat}()},
+\code{\link{rhat}()},
\code{\link{rhat_basic}()},
\code{\link{rhat_nested}()},
-\code{\link{rhat}()},
\code{\link{rstar}()}
}
\concept{diagnostics}
diff --git a/man/extract_variable_array.Rd b/man/extract_variable_array.Rd
index 348a1a74..2c1eae76 100644
--- a/man/extract_variable_array.Rd
+++ b/man/extract_variable_array.Rd
@@ -45,7 +45,7 @@ str(Sigma)
}
\seealso{
Other variable extraction methods:
-\code{\link{extract_variable_matrix}()},
-\code{\link{extract_variable}()}
+\code{\link{extract_variable}()},
+\code{\link{extract_variable_matrix}()}
}
\concept{variable extraction methods}
diff --git a/man/extract_variable_matrix.Rd b/man/extract_variable_matrix.Rd
index 1b9c97c1..dedb49c1 100644
--- a/man/extract_variable_matrix.Rd
+++ b/man/extract_variable_matrix.Rd
@@ -47,7 +47,7 @@ rhat(mu)
}
\seealso{
Other variable extraction methods:
-\code{\link{extract_variable_array}()},
-\code{\link{extract_variable}()}
+\code{\link{extract_variable}()},
+\code{\link{extract_variable_array}()}
}
\concept{variable extraction methods}
diff --git a/man/for_each_draw.Rd b/man/for_each_draw.Rd
index 4c84f31e..531303d8 100644
--- a/man/for_each_draw.Rd
+++ b/man/for_each_draw.Rd
@@ -13,7 +13,8 @@ is defined.}
\item{expr}{(expression) A bare expression that can contain references to
variables in \code{x} by name. This expression will be executed once per draw
of \code{x}, where references to variables in \code{x} resolve to the value of that
-variable in that draw. The expression supports \link{quasiquotation}.}
+variable in that draw. The expression supports
+\link[rlang:topic-inject]{quasiquotation}.}
}
\value{
As \code{for_each_draw()} is used primarily for its side effects (the expression
diff --git a/man/mcse_mean.Rd b/man/mcse_mean.Rd
index c75935b1..96263d29 100755
--- a/man/mcse_mean.Rd
+++ b/man/mcse_mean.Rd
@@ -8,7 +8,7 @@
\usage{
mcse_mean(x, ...)
-\method{mcse_mean}{default}(x, ...)
+\method{mcse_mean}{default}(x, weights = NULL, ...)
\method{mcse_mean}{rvar}(x, ...)
}
@@ -21,6 +21,9 @@ mcse_mean(x, ...)
}}
\item{...}{Arguments passed to individual methods (if applicable).}
+
+\item{weights}{(numeric vector) A vector of weights of length \code{ndraws(x)},
+or \code{NULL} to not have weights. Note that if \code{x} is an \code{rvar} with weights, they are used instead of this argument.}
}
\value{
If the input is an array, returns a single numeric value. If any of the draws
@@ -39,7 +42,8 @@ passing the draws array for that element of the \code{\link{rvar}} to this funct
}
\description{
Compute the Monte Carlo standard error for the mean (expectation) of a
-single variable.
+single variable. If computed on a weighted \code{rvar}, weights will be
+taken into account.
}
\examples{
mu <- extract_variable_matrix(example_draws(), "mu")
@@ -65,9 +69,9 @@ Other diagnostics:
\code{\link{mcse_sd}()},
\code{\link{pareto_diags}()},
\code{\link{pareto_khat}()},
+\code{\link{rhat}()},
\code{\link{rhat_basic}()},
\code{\link{rhat_nested}()},
-\code{\link{rhat}()},
\code{\link{rstar}()}
}
\concept{diagnostics}
diff --git a/man/mcse_quantile.Rd b/man/mcse_quantile.Rd
index 2d05f626..ad5c6025 100755
--- a/man/mcse_quantile.Rd
+++ b/man/mcse_quantile.Rd
@@ -9,7 +9,7 @@
\usage{
mcse_quantile(x, probs = c(0.05, 0.95), ...)
-\method{mcse_quantile}{default}(x, probs = c(0.05, 0.95), names = TRUE, ...)
+\method{mcse_quantile}{default}(x, probs = c(0.05, 0.95), names = TRUE, weights = NULL, ...)
\method{mcse_quantile}{rvar}(x, probs = c(0.05, 0.95), names = TRUE, ...)
@@ -30,6 +30,9 @@ mcse_median(x, ...)
\item{names}{(logical) Should the result have a \code{names} attribute? The
default is \code{TRUE}, but use \code{FALSE} for improved speed if there are many
values in \code{probs}.}
+
+\item{weights}{(numeric vector) A vector of weights of length \code{ndraws(x)},
+or \code{NULL} to not have weights. Note that if \code{x} is an \code{rvar} with weights, they are used instead of this argument.}
}
\value{
If the input is an array,
@@ -52,7 +55,8 @@ result indexes the input probabilities; i.e. the result has dimension
}
\description{
Compute Monte Carlo standard errors for quantile estimates of a
-single variable.
+single variable. If computed on a weighted \code{rvar}, weights will be
+taken into account.
}
\examples{
mu <- extract_variable_matrix(example_draws(), "mu")
@@ -66,7 +70,7 @@ mcse_quantile(d$mu)
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and
Paul-Christian Bürkner (2021). Rank-normalization, folding, and
localization: An improved R-hat for assessing convergence of
-MCMC (with discussion). \emph{Bayesian Data Analysis}. 16(2), 667-–718.
+MCMC (with discussion). \emph{Bayesian Analysis}. 16(2), 667-–718.
doi:10.1214/20-BA1221
}
\seealso{
@@ -80,9 +84,9 @@ Other diagnostics:
\code{\link{mcse_sd}()},
\code{\link{pareto_diags}()},
\code{\link{pareto_khat}()},
+\code{\link{rhat}()},
\code{\link{rhat_basic}()},
\code{\link{rhat_nested}()},
-\code{\link{rhat}()},
\code{\link{rstar}()}
}
\concept{diagnostics}
diff --git a/man/mcse_sd.Rd b/man/mcse_sd.Rd
index 671ef249..4b109d9d 100755
--- a/man/mcse_sd.Rd
+++ b/man/mcse_sd.Rd
@@ -8,7 +8,7 @@
\usage{
mcse_sd(x, ...)
-\method{mcse_sd}{default}(x, ...)
+\method{mcse_sd}{default}(x, weights = NULL, ...)
\method{mcse_sd}{rvar}(x, ...)
}
@@ -21,6 +21,9 @@ mcse_sd(x, ...)
}}
\item{...}{Arguments passed to individual methods (if applicable).}
+
+\item{weights}{(numeric vector) A vector of weights of length \code{ndraws(x)},
+or \code{NULL} to not have weights. Note that if \code{x} is an \code{rvar} with weights, they are used instead of this argument.}
}
\value{
If the input is an array, returns a single numeric value. If any of the draws
@@ -54,7 +57,7 @@ mcse_sd(d$Sigma)
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and
Paul-Christian Bürkner (2021). Rank-normalization, folding, and
localization: An improved R-hat for assessing convergence of
-MCMC (with discussion). \emph{Bayesian Data Analysis}. 16(2), 667-–718.
+MCMC (with discussion). \emph{Bayesian Analysis}. 16(2), 667-–718.
doi:10.1214/20-BA1221
J. F. Kenney & E. S. Keeping (1951). \emph{Mathematics of Statistics, Vol. II.}
@@ -70,9 +73,9 @@ Other diagnostics:
\code{\link{mcse_quantile}()},
\code{\link{pareto_diags}()},
\code{\link{pareto_khat}()},
+\code{\link{rhat}()},
\code{\link{rhat_basic}()},
\code{\link{rhat_nested}()},
-\code{\link{rhat}()},
\code{\link{rstar}()}
}
\concept{diagnostics}
diff --git a/man/pareto_diags.Rd b/man/pareto_diags.Rd
index 46370c49..dcd0194a 100644
--- a/man/pareto_diags.Rd
+++ b/man/pareto_diags.Rd
@@ -68,13 +68,13 @@ The default is \code{"both"}.}
\item{r_eff}{(numeric) relative effective sample size estimate. If
\code{r_eff} is NULL, it will be calculated assuming the draws are
-from MCMC. Default is 1.}
+from MCMC. Default is NULL.}
\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)).}
+al. (2024)).}
\item{verbose}{(logical) Should diagnostic messages be printed? If
\code{TRUE}, messages related to Pareto diagnostics will be
@@ -109,19 +109,19 @@ 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).
+Section 3.2.3, Equation 11 in Vehtari et al. (2024).
\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).
+in Vehtari et al. (2024).
\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).
+theorem convergence rate. See Appendix B in Vehtari et al. (2024).
}
}
\examples{
@@ -133,12 +133,15 @@ 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 (version 8)
+Jonah Gabry (2024). Pareto Smoothed Importance Sampling.
+\emph{Journal of Machine Learning Research}, 25(72):1-58.
+\href{https://jmlr.org/papers/v25/19-556.html}{PDF}
}
\seealso{
-\code{\link{pareto_khat}} for only calculating khat, and
-\code{\link{pareto_smooth}} for Pareto smoothed draws.
+\code{\link{pareto_khat}}, \code{\link{pareto_min_ss}},
+\code{\link{pareto_khat_threshold}}, and \code{\link{pareto_convergence_rate}} for
+individual diagnostics; and \code{\link{pareto_smooth}} for Pareto smoothing
+draws.
Other diagnostics:
\code{\link{ess_basic}()},
@@ -150,9 +153,9 @@ Other diagnostics:
\code{\link{mcse_quantile}()},
\code{\link{mcse_sd}()},
\code{\link{pareto_khat}()},
+\code{\link{rhat}()},
\code{\link{rhat_basic}()},
\code{\link{rhat_nested}()},
-\code{\link{rhat}()},
\code{\link{rstar}()}
}
\concept{diagnostics}
diff --git a/man/pareto_khat.Rd b/man/pareto_khat.Rd
index a4f91707..c419b1d2 100644
--- a/man/pareto_khat.Rd
+++ b/man/pareto_khat.Rd
@@ -41,13 +41,13 @@ The default is \code{"both"}.}
\item{r_eff}{(numeric) relative effective sample size estimate. If
\code{r_eff} is NULL, it will be calculated assuming the draws are
-from MCMC. Default is 1.}
+from MCMC. Default is NULL.}
\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)).}
+al. (2024)).}
\item{verbose}{(logical) Should diagnostic messages be printed? If
\code{TRUE}, messages related to Pareto diagnostics will be
@@ -58,13 +58,25 @@ printed. Default is \code{FALSE}.}
draws are log weights, and only right tail will be smoothed.}
}
\value{
-\code{khat} estimated Generalized Pareto Distribution shape parameter k
+If the input is an array, returns a single numeric value. If any of the draws
+is non-finite, that is, \code{NA}, \code{NaN}, \code{Inf}, or \code{-Inf}, the returned output
+will be (numeric) \code{NA}. Also, if all draws within any of the chains of a
+variable are the same (constant), the returned output will be (numeric) \code{NA}
+as well. The reason for the latter is that, for constant draws, we cannot
+distinguish between variables that are supposed to be constant (e.g., a
+diagonal element of a correlation matrix is always 1) or variables that just
+happened to be constant because of a failure of convergence or other problems
+in the sampling process.
+
+If the input is an \code{\link{rvar}}, returns an array of the same dimensions as the
+\code{\link{rvar}}, where each element is equal to the value that would be returned by
+passing the draws array for that element of the \code{\link{rvar}} to this function.
}
\description{
Estimate Pareto k value by fitting a Generalized Pareto
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).
+diagnostics. For further details see Vehtari et al. (2024).
}
\examples{
mu <- extract_variable_matrix(example_draws(), "mu")
@@ -75,8 +87,9 @@ 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 (version 8)
+Jonah Gabry (2024). Pareto Smoothed Importance Sampling.
+\emph{Journal of Machine Learning Research}, 25(72):1-58.
+\href{https://jmlr.org/papers/v25/19-556.html}{PDF}
}
\seealso{
\code{\link{pareto_diags}} for additional related diagnostics, and
@@ -92,9 +105,9 @@ Other diagnostics:
\code{\link{mcse_quantile}()},
\code{\link{mcse_sd}()},
\code{\link{pareto_diags}()},
+\code{\link{rhat}()},
\code{\link{rhat_basic}()},
\code{\link{rhat_nested}()},
-\code{\link{rhat}()},
\code{\link{rstar}()}
}
\concept{diagnostics}
diff --git a/man/pareto_smooth.Rd b/man/pareto_smooth.Rd
index 24273139..48ae2175 100644
--- a/man/pareto_smooth.Rd
+++ b/man/pareto_smooth.Rd
@@ -13,7 +13,7 @@ pareto_smooth(x, ...)
\method{pareto_smooth}{default}(
x,
tail = c("both", "right", "left"),
- r_eff = 1,
+ r_eff = NULL,
ndraws_tail = NULL,
return_k = FALSE,
extra_diags = FALSE,
@@ -53,13 +53,13 @@ The default is \code{"both"}.}
\item{r_eff}{(numeric) relative effective sample size estimate. If
\code{r_eff} is NULL, it will be calculated assuming the draws are
-from MCMC. Default is 1.}
+from MCMC. Default is NULL.}
\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)).}
+al. (2024)).}
\item{verbose}{(logical) Should diagnostic messages be printed? If
\code{TRUE}, messages related to Pareto diagnostics will be
@@ -72,17 +72,25 @@ draws are log weights, and only right tail will be smoothed.}
\value{
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: * \code{khat}: estimated
-Pareto k shape parameter, and optionally * \code{min_ss}: minimum
-sample size for reliable Pareto smoothed estimate *
-\code{khat_threshold}: khat-threshold for reliable Pareto smoothed
-estimates * \code{convergence_rate}: Relative convergence rate for
+containing numeric values:
+\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}: sample size specific khat threshold for
+reliable Pareto smoothed estimates
+\item \code{convergence_rate}: Relative convergence rate for
Pareto smoothed estimates
}
+
+If any of the draws is non-finite, that is, \code{NA}, \code{NaN}, \code{Inf}, or
+\code{-Inf}, Pareto smoothing will not be performed, and the original
+draws will be returned and and diagnostics will be \code{NA} (numeric).
+}
\description{
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).
+tail(s). For further details see Vehtari et al. (2024).
}
\examples{
mu <- extract_variable_matrix(example_draws(), "mu")
@@ -93,8 +101,9 @@ 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 (version 8)
+Jonah Gabry (2024). Pareto Smoothed Importance Sampling.
+\emph{Journal of Machine Learning Research}, 25(72):1-58.
+\href{https://jmlr.org/papers/v25/19-556.html}{PDF}
}
\seealso{
\code{\link{pareto_khat}} for only calculating khat, and
diff --git a/man/pit.Rd b/man/pit.Rd
new file mode 100644
index 00000000..6294502b
--- /dev/null
+++ b/man/pit.Rd
@@ -0,0 +1,78 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pit.R
+\name{pit}
+\alias{pit}
+\alias{pit.default}
+\alias{pit.draws_matrix}
+\alias{pit.rvar}
+\title{Probability integral transform}
+\usage{
+pit(x, y, ...)
+
+\method{pit}{default}(x, y, weights = NULL, log = FALSE, ...)
+
+\method{pit}{draws_matrix}(x, y, weights = NULL, log = FALSE, ...)
+
+\method{pit}{rvar}(x, y, weights = NULL, log = FALSE, ...)
+}
+\arguments{
+\item{x}{(draws) A \code{\link{draws_matrix}} object or one coercible to a
+\code{draws_matrix} object, or an \code{\link{rvar}} object.}
+
+\item{y}{(observations) A 1D vector, or an array of dim(x), if x is \code{rvar}.
+Each element of \code{y} corresponds to a variable in \code{x}.}
+
+\item{...}{Arguments passed to individual methods (if applicable).}
+
+\item{weights}{A matrix of weights for each draw and variable. \code{weights}
+should have one column per variable in \code{x}, and \code{ndraws(x)} rows.}
+
+\item{log}{(logical) Are the weights passed already on the log scale? The
+default is \code{FALSE}, that is, expecting \code{weights} to be on the standard
+(non-log) scale.}
+}
+\value{
+A numeric vector of length \code{length(y)} containing the PIT values, or
+an array of shape \code{dim(y)}, if \code{x} is an \code{rvar}.
+}
+\description{
+Probability integral transform (PIT). LOO-PIT is given by a weighted sample.
+}
+\details{
+The \code{pit()} function computes the probability integral transform of
+\code{y} using the empirical cumulative distribution computed from the samples
+in \code{x}. For continuous valued \code{y} and \code{x}, the PIT for the elements of \code{y}
+is computed as the empirical cumulative distribution value:
+
+\if{html}{\out{}}\preformatted{PIT(y_i) = Pr(x_i < y_i),
+}\if{html}{\out{}}
+
+where x_i, is the corresponding set of draws in \code{x}. For \code{draws} objects,
+this corresponds to the draws of the \emph{i}th variable, and for \code{rvar}
+the elements of \code{y} and \code{x} are matched.
+
+The draws in \code{x} can further be provided (log-)weights in
+
+If \code{y} and \code{x} are discrete, randomisation is used to obtain continuous PIT
+values. (see, e.g., Czado, C., Gneiting, T., Held, L.: Predictive model
+assessment for count data. Biometrics 65(4), 1254–1261 (2009).)
+}
+\examples{
+# PIT for a draws object
+x <- example_draws()
+# Create a vector of observations
+y <- rnorm(nvariables(x), 5, 5)
+pit(x, y)
+
+# Compute weighted PIT (for example LOO-PIT)
+weights <- matrix(runif(length(x)), ncol = nvariables(x))
+
+pit(x, y, weights)
+
+# PIT for an rvar
+x <- rvar(example_draws())
+# Create an array of observations with the same dimensions as x.
+y_arr <- array(rnorm(length(x), 5, 5), dim = dim(x))
+pit(x, y_arr)
+
+}
diff --git a/man/posterior-package.Rd b/man/posterior-package.Rd
index 5396bf16..e28a88c4 100644
--- a/man/posterior-package.Rd
+++ b/man/posterior-package.Rd
@@ -100,6 +100,7 @@ Other contributors:
\item Ozan Adıgüzel [contributor]
\item Jacob Socolar [contributor]
\item Noa Kallioinen [contributor]
+ \item Teemu Säilynoja [contributor]
}
}
diff --git a/man/ps_convergence_rate.Rd b/man/ps_convergence_rate.Rd
new file mode 100644
index 00000000..535f0ebc
--- /dev/null
+++ b/man/ps_convergence_rate.Rd
@@ -0,0 +1,32 @@
+% 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 convergence rate}
+\usage{
+ps_convergence_rate(k, ndraws, ...)
+}
+\arguments{
+\item{k}{pareto-k values}
+
+\item{ndraws}{number of draws}
+
+\item{...}{unused}
+}
+\value{
+convergence rate
+}
+\description{
+Given number of draws and scalar or array of k's, compute the
+relative convergence rate of PSIS estimate RMSE. See Appendix B of
+Vehtari et al. (2024). This function is exported to be usable by
+other packages. For user-facing diagnostic functions, see
+\code{\link{pareto_convergence_rate}} and \code{\link{pareto_diags}}.
+}
+\seealso{
+Other helper-functions:
+\code{\link{ps_khat_threshold}()},
+\code{\link{ps_min_ss}()},
+\code{\link{ps_tail_length}()}
+}
+\concept{helper-functions}
diff --git a/man/ps_khat_threshold.Rd b/man/ps_khat_threshold.Rd
new file mode 100644
index 00000000..d82ebcb9
--- /dev/null
+++ b/man/ps_khat_threshold.Rd
@@ -0,0 +1,31 @@
+% 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 k-hat threshold}
+\usage{
+ps_khat_threshold(ndraws, ...)
+}
+\arguments{
+\item{ndraws}{number of draws}
+
+\item{...}{unused}
+}
+\value{
+threshold
+}
+\description{
+Given number of draws, computes khat threshold for reliable Pareto
+smoothed estimate (to have small probability of large error). See
+section 3.2.4, equation (13) of Vehtari et al. (2024). This
+function is exported to be usable by other packages. For
+user-facing diagnostic functions, see \code{\link{pareto_khat_threshold}} and
+\code{\link{pareto_diags}}.
+}
+\seealso{
+Other helper-functions:
+\code{\link{ps_convergence_rate}()},
+\code{\link{ps_min_ss}()},
+\code{\link{ps_tail_length}()}
+}
+\concept{helper-functions}
diff --git a/man/ps_min_ss.Rd b/man/ps_min_ss.Rd
new file mode 100644
index 00000000..533cfa34
--- /dev/null
+++ b/man/ps_min_ss.Rd
@@ -0,0 +1,30 @@
+% 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). See
+section 3.2.3 in Vehtari et al. (2024). This function is exported
+to be usable by other packages. For user-facing diagnostic functions, see
+\code{\link{pareto_min_ss}} and \code{\link{pareto_diags}}.
+}
+\seealso{
+Other helper-functions:
+\code{\link{ps_convergence_rate}()},
+\code{\link{ps_khat_threshold}()},
+\code{\link{ps_tail_length}()}
+}
+\concept{helper-functions}
diff --git a/man/ps_tail_length.Rd b/man/ps_tail_length.Rd
new file mode 100644
index 00000000..036afe30
--- /dev/null
+++ b/man/ps_tail_length.Rd
@@ -0,0 +1,30 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pareto_smooth.R
+\name{ps_tail_length}
+\alias{ps_tail_length}
+\title{Pareto tail length}
+\usage{
+ps_tail_length(ndraws, r_eff, ...)
+}
+\arguments{
+\item{ndraws}{number of draws}
+
+\item{r_eff}{relative efficiency}
+
+\item{...}{unused}
+}
+\value{
+tail length
+}
+\description{
+Calculate the tail length from number of draws and relative efficiency
+r_eff. See Appendix H in Vehtari et al. (2024). This function is
+used internally and is exported to be available for other packages.
+}
+\seealso{
+Other helper-functions:
+\code{\link{ps_convergence_rate}()},
+\code{\link{ps_khat_threshold}()},
+\code{\link{ps_min_ss}()}
+}
+\concept{helper-functions}
diff --git a/man/quantile2.Rd b/man/quantile2.Rd
index 3ca81977..45fa5d90 100755
--- a/man/quantile2.Rd
+++ b/man/quantile2.Rd
@@ -6,9 +6,16 @@
\alias{quantile2.rvar}
\title{Compute Quantiles}
\usage{
-quantile2(x, probs = c(0.05, 0.95), na.rm = FALSE, ...)
+quantile2(x, ...)
-\method{quantile2}{default}(x, probs = c(0.05, 0.95), na.rm = FALSE, names = TRUE, ...)
+\method{quantile2}{default}(
+ x,
+ probs = c(0.05, 0.95),
+ na.rm = FALSE,
+ names = TRUE,
+ weights = NULL,
+ ...
+)
\method{quantile2}{rvar}(x, probs = c(0.05, 0.95), na.rm = FALSE, names = TRUE, ...)
}
@@ -20,17 +27,20 @@ quantile2(x, probs = c(0.05, 0.95), na.rm = FALSE, ...)
\item An \code{\link{rvar}}.
}}
+\item{...}{Arguments passed to individual methods (if applicable) and then
+on to \code{\link[stats:quantile]{stats::quantile()}}.}
+
\item{probs}{(numeric vector) Probabilities in \verb{[0, 1]}.}
\item{na.rm}{(logical) Should \code{NA} and \code{NaN} values be removed from \code{x} prior
to computing quantiles? The default is \code{FALSE}.}
-\item{...}{Arguments passed to individual methods (if applicable) and then
-on to \code{\link[stats:quantile]{stats::quantile()}}.}
-
\item{names}{(logical) Should the result have a \code{names} attribute? The
default is \code{TRUE}, but use \code{FALSE} for improved speed if there are many
values in \code{probs}.}
+
+\item{weights}{(numeric vector) A vector of weights of length \code{ndraws(x)},
+or \code{NULL} to not have weights. Note that if \code{x} is an \code{rvar} with weights, they are used instead of this argument.}
}
\value{
A numeric vector of length \code{length(probs)}. If \code{names = TRUE}, it has a
diff --git a/man/rdo.Rd b/man/rdo.Rd
index f3031d5c..b75ca40c 100755
--- a/man/rdo.Rd
+++ b/man/rdo.Rd
@@ -8,7 +8,7 @@ rdo(expr, dim = NULL, ndraws = NULL)
}
\arguments{
\item{expr}{(expression) A bare expression that can (optionally) contain
-\code{\link{rvar}}s. The expression supports \link{quasiquotation}.}
+\code{\link{rvar}}s. The expression supports \link[rlang:topic-inject]{quasiquotation}.}
\item{dim}{(integer vector) One or more integers giving the maximal indices
in each dimension to override the dimensions of the \code{\link{rvar}} to be created
diff --git a/man/rhat.Rd b/man/rhat.Rd
index 263561cd..6af3b66b 100755
--- a/man/rhat.Rd
+++ b/man/rhat.Rd
@@ -54,7 +54,7 @@ rhat(d$Sigma)
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and
Paul-Christian Bürkner (2021). Rank-normalization, folding, and
localization: An improved R-hat for assessing convergence of
-MCMC (with discussion). \emph{Bayesian Data Analysis}. 16(2), 667-–718.
+MCMC (with discussion). \emph{Bayesian Analysis}. 16(2), 667-–718.
doi:10.1214/20-BA1221
}
\seealso{
diff --git a/man/rhat_basic.Rd b/man/rhat_basic.Rd
index 16ffd332..04f47433 100755
--- a/man/rhat_basic.Rd
+++ b/man/rhat_basic.Rd
@@ -62,7 +62,7 @@ Hall/CRC.
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and
Paul-Christian Bürkner (2021). Rank-normalization, folding, and
localization: An improved R-hat for assessing convergence of
-MCMC (with discussion). \emph{Bayesian Data Analysis}. 16(2), 667-–718.
+MCMC (with discussion). \emph{Bayesian Analysis}. 16(2), 667-–718.
doi:10.1214/20-BA1221
}
\seealso{
@@ -77,8 +77,8 @@ Other diagnostics:
\code{\link{mcse_sd}()},
\code{\link{pareto_diags}()},
\code{\link{pareto_khat}()},
-\code{\link{rhat_nested}()},
\code{\link{rhat}()},
+\code{\link{rhat_nested}()},
\code{\link{rstar}()}
}
\concept{diagnostics}
diff --git a/man/rhat_nested.Rd b/man/rhat_nested.Rd
index f2536efd..71f51195 100644
--- a/man/rhat_nested.Rd
+++ b/man/rhat_nested.Rd
@@ -45,7 +45,7 @@ passing the draws array for that element of the \code{\link{rvar}} to this funct
}
\description{
Compute the nested Rhat convergence diagnostic for a single
-variable as proposed in Margossian et al. (2023).
+variable as proposed in Margossian et al. (2024).
}
\details{
Nested Rhat is a convergence diagnostic useful when
@@ -56,8 +56,8 @@ point.
Note that there is a slight difference in the calculation of Rhat
and nested Rhat, as nested Rhat is lower bounded by 1. This means
that nested Rhat with one chain per superchain will not be
-exactly equal to basic Rhat (see Footnote 1 in Margossian et
-al. (2023)).
+exactly equal to basic Rhat (see Footnote 3 in Margossian et
+al. (2024)).
}
\examples{
mu <- extract_variable_matrix(example_draws(), "mu")
@@ -69,9 +69,9 @@ rhat_nested(d$Sigma, superchain_ids = c(1, 1, 2, 2))
}
\references{
Charles C. Margossian, Matthew D. Hoffman, Pavel Sountsov, Lionel
-Riou-Durand, Aki Vehtari and Andrew Gelman (2023). Nested R-hat:
+Riou-Durand, Aki Vehtari and Andrew Gelman (2024). Nested R-hat:
Assessing the convergence of Markov chain Monte Carlo when running
-many short chains. arxiv:arXiv:2110.13017 (version 4)
+many short chains. \emph{Bayesian Analysis}. doi:10.1214/24-BA1453
}
\seealso{
Other diagnostics:
@@ -85,8 +85,8 @@ Other diagnostics:
\code{\link{mcse_sd}()},
\code{\link{pareto_diags}()},
\code{\link{pareto_khat}()},
-\code{\link{rhat_basic}()},
\code{\link{rhat}()},
+\code{\link{rhat_basic}()},
\code{\link{rstar}()}
}
\concept{diagnostics}
diff --git a/man/rstar.Rd b/man/rstar.Rd
index c9479902..f2e74898 100644
--- a/man/rstar.Rd
+++ b/man/rstar.Rd
@@ -117,8 +117,8 @@ Other diagnostics:
\code{\link{mcse_sd}()},
\code{\link{pareto_diags}()},
\code{\link{pareto_khat}()},
+\code{\link{rhat}()},
\code{\link{rhat_basic}()},
-\code{\link{rhat_nested}()},
-\code{\link{rhat}()}
+\code{\link{rhat_nested}()}
}
\concept{diagnostics}
diff --git a/man/subset_draws.Rd b/man/subset_draws.Rd
index 7d040bcf..acd92dec 100644
--- a/man/subset_draws.Rd
+++ b/man/subset_draws.Rd
@@ -22,6 +22,7 @@ subset_draws(x, ...)
regex = FALSE,
unique = TRUE,
exclude = FALSE,
+ scalar = FALSE,
...
)
@@ -34,6 +35,7 @@ subset_draws(x, ...)
regex = FALSE,
unique = TRUE,
exclude = FALSE,
+ scalar = FALSE,
...
)
@@ -46,6 +48,7 @@ subset_draws(x, ...)
regex = FALSE,
unique = TRUE,
exclude = FALSE,
+ scalar = FALSE,
...
)
@@ -58,6 +61,7 @@ subset_draws(x, ...)
regex = FALSE,
unique = TRUE,
exclude = FALSE,
+ scalar = FALSE,
...
)
@@ -70,6 +74,7 @@ subset_draws(x, ...)
regex = FALSE,
unique = TRUE,
exclude = FALSE,
+ scalar = FALSE,
...
)
@@ -108,6 +113,10 @@ how often they appear in the respective selecting arguments.}
If \code{FALSE} (the default) only the selected subset will be
returned. If \code{TRUE} everything but the selected subset will be
returned.}
+
+\item{scalar}{(logical) Should only scalar variables be selected?
+If \code{FALSE} (the default), all variables with matching names and
+\emph{arbitrary} indices will be selected (see examples).}
}
\value{
A \code{draws} object of the same class as \code{x}.
@@ -132,4 +141,7 @@ subset_draws(x, chain = c(1, 1), unique = FALSE)
# extract all elements of 'theta'
subset_draws(x, variable = "theta")
+# trying to extract only a scalar 'theta' will fail
+# subset_draws(x, variable = "theta", scalar = TRUE)
+
}
diff --git a/man/u_scale.Rd b/man/u_scale.Rd
index 87d4edf9..6424636a 100644
--- a/man/u_scale.Rd
+++ b/man/u_scale.Rd
@@ -19,7 +19,9 @@ and dimension as the input.
\description{
Compute rank uniformization for a numeric array. First replace each value by
its rank. Average rank for ties are used to conserve the number of unique
-values of discrete quantities. Second, uniformize ranks to the scale
-\verb{[1/(2S), 1-1/(2S)]}, where \code{S} is the number of values.
+values of discrete quantities. Second, uniformize ranks using formula
+\code{(r - c) / (S - 2 * c + 1)}, where \code{r} is a rank, \code{S} is the number of
+values, and \code{c} is a fractional offset which defaults to c = 3/8 as
+recommend by Blom (1958).
}
\keyword{internal}
diff --git a/tests/testthat/test-as_draws.R b/tests/testthat/test-as_draws.R
index f33720ae..b5d55791 100644
--- a/tests/testthat/test-as_draws.R
+++ b/tests/testthat/test-as_draws.R
@@ -127,6 +127,28 @@ test_that("arrays can be transformed to draws_array objects", {
expect_equal(nchains(y), 4)
})
+test_that("lists of matrices can be transformed to draws_array objects", {
+ x <- round(rnorm(200), 2)
+ x <- matrix(x, nrow = 50)
+ colnames(x) <- paste0("theta", 1:4)
+
+ # one chain
+ z1 <- list(x)
+ y <- as_draws(z1)
+ expect_is(y, "draws_array")
+ expect_equal(variables(y), colnames(z1[[1]]))
+ expect_equal(niterations(y), 50)
+ expect_equal(nchains(y), 1)
+
+ # multiple chains
+ z3 <- list(x, x, x)
+ y <- as_draws(z3)
+ expect_is(y, "draws_array")
+ expect_equal(variables(y), colnames(z3[[1]]))
+ expect_equal(niterations(y), 50)
+ expect_equal(nchains(y), 3)
+})
+
test_that("data.frames can be transformed to draws_df objects", {
x <- data.frame(
v1 = rnorm(100),
diff --git a/tests/testthat/test-convergence.R b/tests/testthat/test-convergence.R
index 1e9034a9..534933d7 100644
--- a/tests/testthat/test-convergence.R
+++ b/tests/testthat/test-convergence.R
@@ -19,7 +19,7 @@ test_that("ess diagnostics return reasonable values", {
expect_true(ess > 250 & ess < 310)
ess <- ess_sd(tau)
- expect_true(ess > 250 & ess < 310)
+ expect_true(ess > 230 & ess < 280)
ess <- ess_bulk(tau)
expect_true(ess > 230 & ess < 280)
diff --git a/tests/testthat/test-log_sum_exp.R b/tests/testthat/test-log_sum_exp.R
new file mode 100644
index 00000000..3e343938
--- /dev/null
+++ b/tests/testthat/test-log_sum_exp.R
@@ -0,0 +1,21 @@
+# Ensure that log_sum_exp agrees with expected behaviour from log(sum(exp(x)))
+test_that("log_sum_exp of for x containing Inf is Inf", {
+ x <- c(Inf, 1)
+ expect_equal(log_sum_exp(x), Inf)
+})
+
+test_that("log_sum_exp of -Inf is -Inf", {
+ x <- c(-Inf, -Inf)
+ expect_equal(log_sum_exp(x), -Inf)
+})
+
+test_that("log_sum_exp of empty input is -Inf", {
+ x <- numeric(0)
+ expect_equal(log_sum_exp(x), -Inf)
+})
+
+test_that("log_sum_exp works", {
+ set.seed(1)
+ x <- log(runif(10))
+ expect_equal(log_sum_exp(x), log(sum(exp(x))))
+})
diff --git a/tests/testthat/test-pareto_smooth.R b/tests/testthat/test-pareto_smooth.R
index abf540a6..78f550f9 100644
--- a/tests/testthat/test-pareto_smooth.R
+++ b/tests/testthat/test-pareto_smooth.R
@@ -1,11 +1,18 @@
-test_that("pareto_khat returns expected reasonable values", {
- tau <- extract_variable_matrix(example_draws(), "tau")
+test_that("pareto_khat handles constant tail correctly", {
+
+ # left tail is constant, so khat should be NA, but for "both" it
+ # should be the same as the right tail
+ x <- c(rep(-100, 10), sort(rnorm(100)))
- pk <- pareto_khat(tau)
- expect_true(names(pk) == "khat")
+ expect_true(is.na(pareto_khat(x, tail = "left", ndraws_tail = 10)))
+ expect_equal(
+ pareto_khat(x, tail = "right", ndraws_tail = 10),
+ pareto_khat(x, tail = "both", ndraws_tail = 10)
+ )
})
+
test_that("pareto_khat handles tail argument", {
# as tau is bounded (0, Inf) the left pareto k should be lower than
@@ -14,8 +21,8 @@ test_that("pareto_khat handles tail argument", {
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)
+ expect_true(pkl < pkr)
+ expect_equal(pkr, pkb)
})
test_that("pareto_khat handles ndraws_tail argument", {
@@ -23,12 +30,12 @@ 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_true(pk10 > pk25)
expect_warning(pareto_khat(tau, tail = "both", ndraws_tail = 201),
- "Number of tail draws cannot be more than half ",
+ paste0("Number of tail draws cannot be more than half ",
"the total number of draws if both tails are fit, ",
- "changing to 200.")
+ "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.")
@@ -41,7 +48,7 @@ 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)
+ expect_true(pk1 < pk0.6)
})
@@ -180,14 +187,27 @@ test_that("pareto_khat functions work with rvars with and without chains", {
})
-test_that("pareto_smooth returns x with smoothed tail", {
- tau <- extract_variable_matrix(example_draws(), "tau")
+test_that("pareto_smooth returns x with smoothed tail(s)", {
+ mu <- extract_variable_matrix(example_draws(), "mu")
+
+ mu_smoothed_right <- pareto_smooth(mu, ndraws_tail = 10, tail = "right", return_k = TRUE)$x
+
+ mu_smoothed_left <- pareto_smooth(mu, ndraws_tail = 10, tail = "left", return_k = TRUE)$x
+
+ mu_smoothed_both <- pareto_smooth(mu, ndraws_tail = 10, tail = "both", return_k = TRUE)$x
+
+ expect_equal(sort(mu)[1:390], sort(mu_smoothed_right)[1:390])
+ expect_equal(sort(mu_smoothed_both)[11:400], sort(mu_smoothed_right)[11:400])
- tau_smoothed <- pareto_smooth(tau, ndraws_tail = 10, tail = "right", return_k = TRUE)$x
+ expect_equal(sort(mu)[11:400], sort(mu_smoothed_left)[11:400])
+ expect_equal(sort(mu_smoothed_both)[1:390], sort(mu_smoothed_left)[1:390])
- expect_equal(sort(tau)[1:390], sort(tau_smoothed)[1:390])
+ expect_false(isTRUE(all.equal(sort(mu), sort(mu_smoothed_left))))
+ expect_false(isTRUE(all.equal(sort(mu), sort(mu_smoothed_right))))
+ expect_false(isTRUE(all.equal(sort(mu), sort(mu_smoothed_both))))
- expect_false(isTRUE(all.equal(sort(tau), sort(tau_smoothed))))
+ expect_false(isTRUE(all.equal(sort(mu_smoothed_both), sort(mu_smoothed_left))))
+ expect_false(isTRUE(all.equal(sort(mu_smoothed_both), sort(mu_smoothed_right))))
})
@@ -222,9 +242,9 @@ test_that("pareto khat works for weighted rvars", {
w2 <- w2 / sum(w2)
xw2 <- weight_draws(xr, weights = w2)
- k <- pareto_khat(xw2)$khat
- kw <- pareto_khat(w2, are_log_weights = TRUE)$khat
- kp <- pareto_khat(draws_of(xw2) * w2)$khat
+ k <- pareto_khat(xw2)
+ kw <- pareto_khat(w2, are_log_weights = TRUE)
+ kp <- pareto_khat(draws_of(xw2) * w2)
expect_true(k > 0.7)
expect_equal(k, max(kw, kp))
diff --git a/tests/testthat/test-pit.R b/tests/testthat/test-pit.R
new file mode 100644
index 00000000..762639b4
--- /dev/null
+++ b/tests/testthat/test-pit.R
@@ -0,0 +1,175 @@
+test_that("validate_y handles numeric inputs correctly", {
+ expect_equal(
+ validate_y(c(1, 2, 3), as_draws_matrix(matrix(1, ncol = 3))),
+ c(1, 2, 3)
+ )
+})
+
+test_that("validate_y raises an error for non-numeric 'y'", {
+ expect_error(validate_y(c("a", "b", "c")), "`y` must be numeric.")
+})
+
+test_that("validate_y raises an error for NAs in 'y'", {
+ expect_error(validate_y(c(1, NA, 3)), "NAs not allowed in `y`.")
+})
+
+test_that("validate_y checks 'y' length against nvariables(x)", {
+ x <- example_draws() |> as_draws_matrix()
+ expect_error(
+ validate_y(1:3, x),
+ "`y` must be a vector of length `nvariables(x)`.",
+ fixed = TRUE
+ )
+})
+
+test_that("validate_y verifies dimensions for array inputs", {
+ x <- posterior::example_draws() |> rvar()
+ y <- posterior::example_draws() |> rvar() |> mean()
+ expect_equal(validate_y(y, x), y)
+
+ y_mismatched <- array(1:8, dim = c(2, 2, 2))
+ expect_error(
+ validate_y(y_mismatched, x),
+ "`dim(y)` must match `dim(x)`.",
+ fixed = TRUE
+ )
+})
+
+# test normalize_log_weights
+test_that("normalize_log_weights returns log-normalized columns", {
+ set.seed(1)
+ x <- matrix(log(runif(200)), ncol = 10)
+ result <- colSums(exp(normalize_log_weights(x)))
+ expect_equal(result, rep(1, 10))
+})
+
+# tests for pit.default
+test_that("pit works without weights", {
+ x <- matrix(c(1, 2, 3, 5), nrow = 2)
+ y <- c(3, 4)
+ result <- pit(x, y)
+
+ expect_length(result, ncol(x))
+ expect_equal(unname(result), c(1, 0.5))
+})
+
+test_that("pit.default works with log-weights", {
+ x <- matrix(c(1, 2, 3, 5), nrow = 2)
+ y <- c(3, 4)
+ weights <- matrix(log(c(1, 1, 1, 1)), nrow = 2)
+ result <- pit.default(x, y, weights = weights, log = TRUE)
+
+ expect_length(result, ncol(x))
+ expect_equal(unname(result), c(1, 0.5))
+})
+
+test_that("pit.default works with non-log weights", {
+ x <- matrix(c(1, 2, 3, 5), nrow = 2)
+ y <- c(3, 4)
+ weights <- matrix(c(1, 1, 1, 1), nrow = 2)
+ result <- pit.default(x, y, weights = weights)
+
+ expect_length(result, ncol(x))
+ expect_equal(unname(result), c(1, 0.5))
+})
+
+test_that("pit.default handles randomized PIT with no weights", {
+ x <- matrix(c(rep(c(0, 1), 1000), rep(c(3, 3), 1000)), nrow = 2)
+ y <- rep(c(1, 3), each = 1000)
+ result <- pit.default(x, y)
+
+ expect_length(result, ncol(x))
+ expect_true(all(result[1:1000] >= .5 & result[1:1000] <= 1))
+ expect_true(all(result[1001:2000] >= 0 & result[1001:2000] <= 1))
+})
+
+test_that("pit.default handles randomized PIT with weights", {
+ x <- matrix(c(rep(c(0, 1), 1000), rep(c(3, 3), 1000)), nrow = 2)
+ y <- rep(c(1, 3), each = 1000)
+ weights <- matrix(.5, nrow = 2, ncol = 2000)
+ result <- pit.default(x, y, weights)
+
+ expect_length(result, ncol(x))
+ expect_true(all(result[1:1000] >= .5 & result[1:1000] <= 1))
+ expect_true(all(result[1001:2000] >= 0 & result[1001:2000] <= 1))
+})
+
+test_that("pit.default handles randomized PIT with log-weights", {
+ x <- matrix(c(rep(c(0, 1), 1000), rep(c(3, 3), 1000)), nrow = 2)
+ y <- rep(c(1, 3), each = 1000)
+ weights <- matrix(log(.5), nrow = 2, ncol = 2000)
+ result <- pit.default(x, y, weights, log = TRUE)
+
+ expect_length(result, ncol(x))
+ expect_true(all(result[1:1000] >= .5 & result[1:1000] <= 1))
+ expect_true(all(result[1001:2000] >= 0 & result[1001:2000] <= 1))
+})
+
+test_that("pit works with draws objects", {
+ set.seed(1)
+ n_chains <- 4
+ n_vars <- 3
+ n_draws <- 100
+
+ test_array <- array(
+ rnorm(n_draws * n_vars * n_chains),
+ dim = c(n_draws, n_chains, n_vars)
+ )
+
+ y <- rnorm(n_vars)
+
+ pit_true <- sapply(1:n_vars, \(v) mean(test_array[, , v] < y[v]))
+
+ x <- as_draws(test_array)
+
+ expect_equal(unname(pit(x, y)), pit_true)
+ expect_equal(unname(pit(as_draws_df(test_array), y)), pit_true)
+ expect_equal(unname(pit(as_draws_matrix(test_array), y)), pit_true)
+ expect_equal(unname(pit(as_draws_list(test_array), y)), pit_true)
+ expect_equal(unname(pit(as_draws_rvars(test_array), y)), pit_true)
+})
+
+test_that("pit works with rvars", {
+ set.seed(1)
+ n_col <- 4
+ n_row <- 3
+ n_draws <- 100
+
+ test_array <- array(
+ rnorm(n_draws * n_col * n_row),
+ dim = c(n_draws, n_row, n_col)
+ )
+
+ weights <- matrix(
+ runif(n_draws * n_col * n_row),
+ ncol = n_row * n_col
+ )
+
+ y <- array(rnorm(n_col * n_row), dim = c(n_row, n_col))
+
+ result <- pit(rvar(test_array), y)
+ expect_equal(dim(result), dim(y))
+ expect_true(all(
+ result ==
+ array(
+ pit(array(test_array, dim = c(n_draws, n_col * n_row)), c(y)),
+ dim(y)
+ )
+ ))
+
+ result <- pit(rvar(test_array), y, weights)
+ expect_equal(dim(result), dim(y))
+ expect_true(all(result == array(
+ pit(
+ array(test_array, dim = c(n_draws, n_col * n_row)),
+ c(y),
+ array(weights, dim = c(n_draws, n_col * n_row))
+ ),
+ dim(y)
+ )))
+})
+
+test_that("pit doesn't error for empty draws", {
+ expect_numeric(pit(empty_draws_array(), numeric(0)))
+ expect_numeric(pit(rvar(), numeric(0)))
+})
diff --git a/tests/testthat/test-rstar.R b/tests/testthat/test-rstar.R
index 65a47517..be5169df 100644
--- a/tests/testthat/test-rstar.R
+++ b/tests/testthat/test-rstar.R
@@ -29,6 +29,7 @@ test_that("rstar works with draws_df example", {
test_that("rstar with uncertainty returns vectors of correct length", {
skip_if_not_installed("caret")
+ skip_if_not_installed("gbm")
x <- example_draws()
val <- rstar(x, method = "gbm", uncertainty = T, verbose = F)
expect_equal(length(val), 1000)
@@ -45,6 +46,7 @@ test_that("incorrect nsimulations values throws error", {
test_that("rstar with uncertainty returns reasonable values", {
skip_if_not_installed("caret")
+ skip_if_not_installed("gbm")
x <- example_draws()
val <- rstar(x, method = "gbm", uncertainty = T, verbose = F)
expect_true(max(val) > 0.3 & min(val) < 10)
@@ -52,6 +54,7 @@ test_that("rstar with uncertainty returns reasonable values", {
test_that("rstar accepts different classifiers", {
skip_if_not_installed("caret")
+ skip_if_not_installed("gbm")
x <- example_draws()
val <- rstar(x, method = "gbm", verbose=F)
expect_true(is.numeric(val))
@@ -61,6 +64,7 @@ test_that("rstar accepts different classifiers", {
test_that("rstar accepts different hyperparameters", {
skip_if_not_installed("caret")
+ skip_if_not_installed("gbm")
x <- example_draws()
# use fast hyperparameters
diff --git a/tests/testthat/test-subset_draws.R b/tests/testthat/test-subset_draws.R
index b5e44602..a1d500dc 100644
--- a/tests/testthat/test-subset_draws.R
+++ b/tests/testthat/test-subset_draws.R
@@ -24,7 +24,10 @@ test_that("subset_draws works correctly for draws_matrix objects", {
x_sub <- subset_draws(x, iteration = c(1, 2, 3), exclude = TRUE)
expect_equal(ndraws(x) - 3 * nchains(x), ndraws(x_sub))
-
+ expect_error(
+ subset_draws(x, variable = "theta", scalar = TRUE),
+ "The following variables are missing in the draws object"
+ )
})
test_that("subset_draws works correctly for draws_array objects", {
diff --git a/tests/testthat/test-weight_draws.R b/tests/testthat/test-weight_draws.R
index cd94a7aa..86568c47 100644
--- a/tests/testthat/test-weight_draws.R
+++ b/tests/testthat/test-weight_draws.R
@@ -164,8 +164,10 @@ test_that("weights must match draws", {
}
})
-test_that("weights must be a vector, not array/matrix", {
- x <- example_draws()
- w <- seq_len(ndraws(x))
- expect_error(weight_draws(x, matrix(w)), "Must be.*vector.*not.*matrix")
-})
+
+## pit now allows for weights to be a draws matrix, so this no longer passes
+## test_that("weights must be a vector, not array/matrix", {
+## x <- example_draws()
+## w <- seq_len(ndraws(x))
+## expect_error(weight_draws(x, matrix(w)), "Must be.*vector.*not.*matrix")
+## })
diff --git a/vignettes/pareto_diagnostics.Rmd b/vignettes/pareto_diagnostics.Rmd
new file mode 100644
index 00000000..868e9f80
--- /dev/null
+++ b/vignettes/pareto_diagnostics.Rmd
@@ -0,0 +1,244 @@
+---
+title: "Pareto-khat diagnostics"
+author: "Aki Vehtari"
+output:
+ rmarkdown::html_vignette:
+ toc: true
+ toc_depth: 3
+vignette: >
+ %\VignetteIndexEntry{Pareto-khat diagnostics}
+ %\VignetteEngine{knitr::rmarkdown}
+ %\VignetteEncoding{UTF-8}
+---
+
+## Introduction
+
+The paper
+
+* Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao, and
+ Jonah Gabry (2024). Pareto smoothed importance sampling.
+ *Journal of Machine Learning Research*, 25(72):1-58.
+
+presents Pareto smoothed importance sampling, but also
+Pareto-$\hat{k}$ diagnostic that can be used when estimating any
+expectation based on finite sample. This vignette illustrates the use of
+these diagnostics.
+The individual diagnostic functions are `pareto_khat()`, `pareto_min_ss()`, `pareto_convergence_rate()` and `pareto_khat_threshold()`. The function `pareto_diags()` will return all of these.
+
+Additionally, the `pareto_smooth()` function can be used to transform draws by smoothing the tail(s).
+## Example
+
+```{r setup}
+library(posterior)
+library(dplyr)
+options(pillar.neg = FALSE, pillar.subtle=FALSE, pillar.sigfig=2)
+```
+
+### Simulated data
+
+Generate `xn` a simulated MCMC sample with 4 chains each with 1000
+iterations using AR process with marginal normal(0,1)
+
+```{r simulate-data-1}
+N <- 1000
+phi <- 0.3
+set.seed(6534)
+dr <- array(data=replicate(4,as.numeric(arima.sim(n = N,
+ list(ar = c(phi)),
+ sd = sqrt((1-phi^2))))),
+ dim=c(N,4,1)) %>%
+ as_draws_df() %>%
+ set_variables('xn')
+```
+
+Transform `xn` via cdf-inverse-cdf so that we have variables that
+have marginally distributions $t_3$, $t_{2.5}$, $t_2$, $t_{1.5}$,
+and $t_1$. These all have thick tails. In addition $t_2$,
+$t_{1.5}$, and $t_1$ have infinite variance, and $t_1$ (aka Cauchy)
+has infinite mean.
+
+```{r simulate-data-2}
+drt <- dr %>%
+ mutate_variables(xt3=qt(pnorm(xn), df=3),
+ xt2_5=qt(pnorm(xn), df=2.5),
+ xt2=qt(pnorm(xn), df=2),
+ xt1_5=qt(pnorm(xn), df=1.5),
+ xt1=qt(pnorm(xn), df=1))
+```
+
+### MCMC convergence diagnostics
+
+We examine the draws with the default `summarise_draws()`.
+
+```{r summarise_draws}
+drt %>%
+ summarise_draws()
+```
+
+All the usual convergence diagnostics $\widehat{R}$, Bulk-ESS, and
+Tail-ESS look good, which is fine as they have been designed to
+work also with infinite variance (Vehtari et al., 2020).
+
+If these variables would present variables of interest for which
+we would like to estimate means, we would be also interested in
+Monte Carlo standard error (MCSE, see case study [How many
+iterations to run and how many digits to
+report](https://users.aalto.fi/~ave/casestudies/Digits/digits.html)).
+
+```{r summarise_draws-mcse}
+drt %>%
+ summarise_draws(mean, sd, mcse_mean, ess_bulk, ess_basic)
+```
+
+Here MCSE for mean is based on standard deviation and Basic-ESS,
+but these assume finite variance. We did sample also from
+distributions with infinite variance, but given a finite sample size, the
+empirical variance estimates are always finite, and thus we get
+overoptimistic MCSE.
+
+## Pareto-$\hat{k}$
+
+To diagnose whether our variables of interest may have infinite
+variance and even infinite mean, we can use Pareto-$\hat{k}$
+diagnostic.
+
+```{r summarise_draws-pareto_khat}
+drt %>%
+ summarise_draws(mean, sd, mcse_mean, ess_basic, pareto_khat)
+```
+
+$\hat{k} \leq 0$ indicates that all moments exist, and the inverse
+of positive $\hat{k}$ tells estimate for the number of finite (fractional)
+moments. Thus, $\hat{k}\geq 1/2$ indicates infinite variance,
+and $\hat{k}\geq 1$ indicates infinite mean. Sometimes very thick
+distribution tails may affect also sampling, but assuming sampling
+did go well, and we would be interested only in quantiles, infinite
+variance and mean are not a problem. But if we are interested in mean,
+then we need to care about the number of (fractional) moments. Here we
+see $\hat{k} \geq 1/2$ for $t_2$, $t_{1.5}$, and $t_{1}$, and
+we should not trust their `mcse_mean` values. Without trustworthy MCSE
+estimate we don't have good estimate of how accurate the mean estimate is.
+Furthermore, as $\hat{k} \geq 1$ for $t_{1}$, the mean is not finite and
+the mean estimate is not valid.
+
+## Pareto smoothing
+
+If we really do need those mean estimates, we can improve
+trustworthiness by Pareto smoothing, which replaces
+extreme tail draws with expected ordered statistics of Pareto
+distribution fitted to the tails of the distribution. Pareto
+smoothed mean estimate (computed using Pareto smoothed draws) has
+finite variance with a cost of some bias which we know when it is
+negligible. As a thumb rule when $\hat{k}<0.7$, the bias is
+negligible.
+
+We do Pareto smoothing for all the variables.
+
+```{r pareto_smooth}
+drts <- drt %>%
+ mutate_variables(xt3_s=pareto_smooth(xt3),
+ xt2_5_s=pareto_smooth(xt2_5),
+ xt2_s=pareto_smooth(xt2),
+ xt1_5_s=pareto_smooth(xt1_5),
+ xt1_s=pareto_smooth(xt1)) %>%
+ subset_draws(variable="_s", regex=TRUE)
+```
+
+Now the `mcse_mean` values are more trustworthy when $\hat{k} < 0.7$.
+When $\hat{k}>0.7$ both bias and variance grow so fast that Pareto smoothing
+rarely helps (see more details in the paper).
+
+```{r summarise_draws-pareto_khat-2}
+drts %>%
+ summarise_draws(mean, mcse_mean, ess_basic, pareto_khat)
+```
+
+
+## Minimum sample size required
+
+The bias and variance depend on the sample size, and we can
+use additional diagnostic `min_ss` which tells the minimum sample size needed
+so that `mcse_mean` can be trusted.
+
+```{r summarise_draws-min_ss}
+drt %>%
+ summarise_draws(mean, mcse_mean, ess_basic,
+ pareto_khat, min_ss=pareto_min_ss)
+```
+
+Here required `min_ss` is smaller than `ess_basic` for all except $t_1$, for
+which there is no hope.
+
+## Convergence rate
+
+Given finite variance, the central limit theorem states that to halve
+MCSE we need four times bigger sample size. With Pareto smoothing,
+we can go further, but the convergence rate decreases when $\hat{k}$ increases.
+
+```{r summarise_draws-conv_rate}
+drt %>%
+ summarise_draws(mean, mcse_mean, ess_basic,
+ pareto_khat, min_ss=pareto_min_ss,
+ conv_rate=pareto_convergence_rate)
+```
+
+We see that with $t_2$, $t_{1.5}$, and $t_1$ we need $4^{1/0.86}\approx 5$,
+$4^{1/0.60}\approx 10$, and $4^{1/0}\approx \infty$ times bigger sample sizes to
+halve MCSE for mean.
+
+## Pareto-$\hat{k}$-threshold
+
+The final Pareto diagnostic, $\hat{k}$-threshold, is useful for smaller sample sizes. Here we select only 100 iterations per chain to get total of 400 draws.
+
+```{r summarise_draws-khat_threshold}
+drt %>%
+ subset_draws(iteration=1:100) %>%
+ summarise_draws(mean, mcse_mean, ess_basic,
+ pareto_khat, min_ss=pareto_min_ss,
+ khat_thres=pareto_khat_threshold,
+ conv_rate=pareto_convergence_rate)
+```
+
+With only 400 draws, we can trust the Pareto smoothed result only when
+$\hat{k}<0.62$. For $t_{1.5}$ $\hat{k}\approx 0.64$, and `min_ss` reveals
+we would probably need more than 560 draws to be on the safe side.
+
+## Pareto diagnostics
+
+We can get all these diagnostics with `pareto_diags()`, and it's
+easy to use it also for derived quantities.
+
+```{r summarise_draws-pareto_diags}
+drt %>%
+ mutate_variables(xt2_5_sq=xt2_5^2) %>%
+ subset_draws(variable="xt2_5_sq") %>%
+ summarise_draws(mean, mcse_mean,
+ pareto_diags)
+```
+
+## Discussion
+
+All these diagnostics are presented in Section 3 and summarized in
+Table 1 in PSIS paper (Vehtari et al., 2024).
+
+If you don't need to estimate means of thick tailed distributions,
+and there are no sampling issues due to thick tails, then you don't
+need to check existence of finite variance, and thus there is no
+need to check Pareto-$\hat{k}$ for all the parameters and derived
+quantities.
+
+It is possible that the distribution has finite variance, but
+pre-asymptotically given a finite sample size the behavior can be
+similar to infinite variance. Thus the diagnostic is useful even in
+cases where theory guarantees finite variance.
+
+## Reference
+
+Vehtari, A., Simpson, D., Gelman, A., Yao, Y., & Gabry, J. (2024).
+Pareto smoothed importance sampling.
+*Journal of Machine Learning Research*, 25(72):1-58.
+
+Vehtari A., Gelman A., Simpson D., Carpenter B., & Bürkner P. C. (2020).
+Rank-normalization, folding, and localization: An improved Rhat for assessing
+convergence of MCMC. *Bayesian Analysis*, 16(2):667-718.
+
\ No newline at end of file
diff --git a/vignettes/posterior.Rmd b/vignettes/posterior.Rmd
index 85584539..a1e2cb75 100644
--- a/vignettes/posterior.Rmd
+++ b/vignettes/posterior.Rmd
@@ -396,11 +396,11 @@ with `draws` objects, please open an issue at
## References
-Gelman A., Carlin J. B., Stern H. S., David B. Dunson D. B., Aki Vehtari A.,
+Gelman A., Carlin J. B., Stern H. S., David B. Dunson D. B., Vehtari, A.,
& Rubin D. B. (2013). *Bayesian Data Analysis, Third Edition*. Chapman and
Hall/CRC.
Vehtari A., Gelman A., Simpson D., Carpenter B., & Bürkner P. C. (2020).
Rank-normalization, folding, and localization: An improved Rhat for assessing
-convergence of MCMC. *Bayesian Analysis*.
+convergence of MCMC. *Bayesian Analysis*, 16(2):667-718.