Tools for working with posterior (and prior) distributions

Description

https://mc-stan.org/posterior/

The posterior package is intended to provide 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:

• Efficiently convert between many different useful formats of draws (samples) from posterior or prior distributions.

• Provide consistent methods for operations commonly performed on draws, for example, subsetting, binding, or mutating draws.

• Provide various summaries of draws in convenient formats.

• Provide lightweight implementations of state of the art posterior inference diagnostics.

Package options

The following options are used to format and print draws objects, as in print.draws_array(), print.draws_df(), print.draws_list(), print.draws_matrix(), and print.draws_rvars():

• posterior.max_draws: Maximum number of draws to print.

• posterior.max_iterations: Maximum number of iterations to print.

• posterior.max_chains: Maximum number of chains to print.

• posterior.max_variables: Maximum number of variables to print.

The following options are used for formatting the output of summarize_draws:

• posterior.num_args: Arguments passed to num() for pretty printing of summaries.

The following options are used to format and print rvar objects, as in print.rvar() and print.draws_rvars():

• posterior.rvar_summary: What style of summary to display: "mean_sd" displays ⁠mean ± sd⁠, "median_mad" displays ⁠median ± mad⁠.

• posterior.digits: How many significant digits are displayed. This defaults to a smaller value (2) than getOption("digits") because rvars print two numbers (point summary and uncertainty) next to each other.

The following option is used to construct new rvar objects, as in rfun() and rdo():

• posterior.rvar_ndraws: The number of draws used to construct new random variables when this number cannot be determined from existing arguments (e.g., other rvars passed to a function).

The following options are used to control warning messages:

• posterior.warn_on_merge_chains: (logical) Some operations will 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.

Author(s)

Maintainer: Paul-Christian Bürkner paul.buerkner@gmail.com

Authors:

• Jonah Gabry jsg2201@columbia.edu

• Matthew Kay mjskay@northwestern.edu

• Aki Vehtari Aki.Vehtari@aalto.fi

Other contributors:

• Måns Magnusson [contributor]

• Rok Češnovar [contributor]

• Ben Lambert [contributor]

• Ozan Adıgüzel [contributor]

• Jacob Socolar [contributor]

• Noa Kallioinen [contributor]

See Also

Useful links:

• https://mc-stan.org/posterior/

• https://discourse.mc-stan.org/

• Report bugs at https://github.com/stan-dev/posterior/issues

Extract parts of a draws_array object

Description

Extract parts of a draws_array object. They are strictly defined as arrays of 3 dimensions (iteration x chain x variable) so dropping any of the dimensions breaks the expected structure of the object. Accordingly, no dropping of dimensions is done by default even if the extracted slices are of length 1. If drop is manually set to TRUE and any of the dimensions is actually dropped, this will lead to dropping the "draws_array" class as well.

Usage

## S3 method for class 'draws_array'
x[i, j, ..., drop = FALSE]

Arguments

Value

An object of class "draws_array" unless any of the dimensions was dropped during the extraction.

Extract parts of a draws_matrix object

Description

Extract parts of a draws_matrix object. They are strictly defined as matrices (draws x variable) so dropping any of the dimensions breaks the expected structure of the object. Accordingly, no dropping of dimensions is done by default even if the extracted slices are of length 1. If drop is manually set to TRUE and any of the dimensions is actually dropped, this will lead to dropping the "draws_matrix" class as well.

Usage

## S3 method for class 'draws_matrix'
x[i, j, ..., drop = FALSE]

Arguments

Value

An object of class "draws_matrix" unless any of the dimensions was dropped during the extraction.

Coerce to a random variable

Description

Convert x to an rvar object.

Usage

as_rvar(x, dim = NULL, dimnames = NULL, nchains = NULL)

as_rvar_numeric(x, dim = NULL, dimnames = NULL, nchains = NULL)

as_rvar_integer(x, dim = NULL, dimnames = NULL, nchains = NULL)

as_rvar_logical(x, dim = NULL, dimnames = NULL, nchains = NULL)

Arguments

Details

For objects that are already rvars, returns them (with modified dimensions if dim is not NULL).

For numeric or logical vectors or arrays, returns an rvar with a single draw and the same dimensions as x. This is in contrast to the rvar() constructor, which treats the first dimension of x as the draws dimension. As a result, as_rvar() is useful for creating constants.

While as_rvar() attempts to pick the most suitable subtype of rvar based on the type of x (possibly returning an rvar_factor or rvar_ordered), as_rvar_numeric(), as_rvar_integer(), and as_rvar_logical() always coerce the draws of the output rvar to be numeric, integer, or logical (respectively), and always return a base rvar, never a subtype.

Value

An object of class "rvar" (or one of its subtypes) representing a random variable.

See Also

rvar() to construct rvars directly. See rdo(), rfun(), and rvar_rng() for higher-level interfaces for creating rvars.

Examples

# You can use as_rvar() to create "constant" rvars (having only one draw):
x <- as_rvar(1)
x

# Such constants can be of arbitrary shape:
as_rvar(1:4)
as_rvar(matrix(1:10, nrow = 5))
as_rvar(array(1:12, dim = c(2, 3, 2)))

# as_rvar_numeric() coerces subtypes of rvar to the base rvar type
y <- as_rvar_factor(c("a", "b", "c"))
y
as_rvar_numeric(y)

Coerce to a factor random variable

Description

Convert x to an rvar_factor or rvar_ordered object.

Usage

as_rvar_factor(x, dim = NULL, dimnames = NULL, nchains = NULL, ...)

as_rvar_ordered(x, dim = NULL, dimnames = NULL, nchains = NULL, ...)

Arguments

Details

For objects that are already rvars, returns them (with modified dimensions if dim is not NULL), possibly adding levels using the unique values of the draws of the rvar (if the object is not already factor-like).

For numeric, logical, factor, or character vectors or arrays, returns an rvar_factor or rvar_ordered with a single draw and the same dimensions as x. This is in contrast to the rvar_factor() and rvar_ordered() constructors, which treats the first dimension of x as the draws dimension. As a result, as_rvar_factor() and as_rvar_ordered() are useful for creating constants.

Value

An object of class "rvar_factor" or "rvar_ordered" representing a random variable.

See Also

rvar(), rvar_factor(), and rvar_ordered() to construct rvars directly. See rdo(), rfun(), and rvar_rng() for higher-level interfaces for creating rvars.

Examples

# You can use as_rvar_factor() to create "constant" rvars (having only one draw):
x <- as_rvar_factor("a")
x

# Such constants can be of arbitrary shape:
as_rvar_factor(letters[1:4])
as_rvar_ordered(matrix(letters[1:10], nrow = 5))
as_rvar_factor(array(letters[1:12], dim = c(2, 3, 2)))

Autocorrelation estimates

Description

Compute autocorrelation estimates for every lag for the specified input sequence using a fast Fourier transform approach. The estimate for lag t is scaled by N-t where N is the length of the sequence.

Usage

autocorrelation(x)

Arguments

Value

A numeric vector of autocorrelations at every lag (scaled by N-lag).

Autocovariance estimates

Description

Compute autocovariance estimates for every lag for the specified input sequence using a fast Fourier transform approach. The estimate for lag t is scaled by N-t where N is the length of the sequence.

Usage

autocovariance(x)

Arguments

Value

A numeric vector of autocovariances at every lag (scaled by N-lag).

Bind draws objects together

Description

Bind multiple draws objects together to form a single draws object.

Usage

bind_draws(x, ...)

## S3 method for class 'draws_matrix'
bind_draws(x, ..., along = "variable")

## S3 method for class 'draws_array'
bind_draws(x, ..., along = "variable")

## S3 method for class 'draws_df'
bind_draws(x, ..., along = "variable")

## S3 method for class 'draws_list'
bind_draws(x, ..., along = "variable")

## S3 method for class 'draws_rvars'
bind_draws(x, ..., along = "variable")

Arguments

Value

A draws object of the same class as x.

Examples

x1 <- draws_matrix(alpha = rnorm(5), beta = rnorm(5))
x2 <- draws_matrix(alpha = rnorm(5), beta = rnorm(5))
ndraws(x1)
ndraws(x2)
x3 <- bind_draws(x1, x2, along = "draw")
ndraws(x3)

x4 <- draws_matrix(theta = rexp(5))
x5 <- bind_draws(x1, x4, along = "variable")
variables(x5)

Cholesky decomposition of random matrix

Description

Cholesky decomposition of an rvar containing a matrix.

Usage

## S3 method for class 'rvar'
chol(x, ...)

Arguments

Value

An rvar containing the upper triangular factor of the Cholesky decomposition, i.e., the matrix $R$ such that $R'R = x$.

Matrix diagonals (including for random variables)

Description

Extract the diagonal of a matrix or construct a matrix, including random matrices (2-dimensional rvars). Makes base::diag() generic.

Usage

## S4 method for signature 'rvar'
diag(x = 1, nrow, ncol, names = TRUE)

Arguments

Details

Makes base::diag() into a generic function. See that function's documentation for usage with numerics and for usage of diag<-, which is also supported by rvar.

Value

For rvars, has two modes:

1. x is a matrix-like rvar: it returns the diagonal as a vector-like rvar

2. x is a vector-like rvar: it returns a matrix-like rvar with x as the diagonal and zero for off-diagonal entries.

See Also

base::diag()

Examples

# Sigma is a 3x3 covariance matrix
Sigma <- as_draws_rvars(example_draws("multi_normal"))$Sigma
Sigma

diag(Sigma)

diag(Sigma) <- 1:3
Sigma

diag(as_rvar(1:3))

List of available convergence diagnostics

Description

A list of available diagnostics and links to their individual help pages.

Details

Value

See individual functions for a description of return types.

Dissention

Description

Dissention, for measuring dispersion in draws from ordinal distributions.

Usage

dissent(x)

## Default S3 method:
dissent(x)

## S3 method for class 'rvar'
dissent(x)

Arguments

Details

Calculates Tastle and Wierman's (2007) dissention measure:

$-\sum_{i = 1}^{n} p_i \log_2 \left(1 - \frac{|x_i - \mathrm{E}(x)| }{\max(x) - \min(x)} \right)$

This ranges from 0 (all probability in one category) through 0.5 (uniform) to 1 (bimodal: all probability split equally between the first and last category).

Value

If x is a factor or numeric, returns a length-1 numeric vector with a value between 0 and 1 (inclusive) giving the dissention of x.

If x is an rvar, returns an array of the same shape as x, where each cell is the dissention of the draws in the corresponding cell of x.

References

William J. Tastle, Mark J. Wierman (2007). Consensus and dissention: A measure of ordinal dispersion. International Journal of Approximate Reasoning. 45(3), 531–545. doi:10.1016/j.ijar.2006.06.024.

Examples

set.seed(1234)

levels <- c("lowest", "low", "neutral", "high", "highest")

# a bimodal distribution: high dissention
x <- ordered(
  sample(levels, 4000, replace = TRUE, prob = c(0.45, 0.04, 0.02, 0.04, 0.45)),
  levels = levels
)
dissent(x)

# a unimodal distribution: low dissention
y <- ordered(
  sample(levels, 4000, replace = TRUE, prob = c(0.95, 0.02, 0.015, 0.01, 0.005)),
  levels = levels
)
dissent(y)

# both together, as an rvar
xy <- c(rvar(x), rvar(y))
xy
dissent(xy)

Transform to draws objects

Description

Try to transform an R object to a format supported by the posterior package.

Usage

as_draws(x, ...)

is_draws(x)

Arguments

Details

The class "draws" is the parent class of all supported formats, which also have their own subclasses of the form "draws_{format}" (e.g. "draws_array").

Value

If possible, a draws object in the closest supported format to x. The formats are linked to in the See Also section below.

See Also

Other formats: draws_array(), draws_df(), draws_list(), draws_matrix(), draws_rvars()

Examples

# create some random draws
x <- matrix(rnorm(30), nrow = 10)
colnames(x) <- c("a", "b", "c")
str(x)

# transform to a draws object
y <- as_draws(x)
str(y)

# remove the draws classes from the object
class(y) <- class(y)[-(1:2)]
str(y)

The draws_array format

Description

The as_draws_array() methods convert objects to the draws_array format. The draws_array() function creates an object of the draws_array format based on a set of numeric vectors. See Details.

Usage

as_draws_array(x, ...)

## Default S3 method:
as_draws_array(x, ...)

## S3 method for class 'draws_array'
as_draws_array(x, ...)

## S3 method for class 'draws_matrix'
as_draws_array(x, ...)

## S3 method for class 'draws_df'
as_draws_array(x, ...)

## S3 method for class 'draws_list'
as_draws_array(x, ...)

## S3 method for class 'draws_rvars'
as_draws_array(x, ...)

## S3 method for class 'mcmc'
as_draws_array(x, ...)

## S3 method for class 'mcmc.list'
as_draws_array(x, ...)

draws_array(..., .nchains = 1)

is_draws_array(x)

Arguments

Details

Objects of class "draws_array" are 3-D arrays with dimensions "iteration", "chain", and "variable". See Examples.

Value

A draws_array object, which has classes c("draws_array", "draws", "array").

See Also

Other formats: draws, draws_df(), draws_list(), draws_matrix(), draws_rvars()

Examples

x1 <- as_draws_array(example_draws())
class(x1)
print(x1)
str(x1)

x2 <- draws_array(a = rnorm(10), b = rnorm(10), c = 1)
class(x2)
print(x2)
str(x2)

The draws_df format

Description

The as_draws_df() methods convert objects to the draws_df format. The draws_df() function creates an object of the draws_df format based on a set of numeric vectors. See Details.

Usage

as_draws_df(x, ...)

## Default S3 method:
as_draws_df(x, ...)

## S3 method for class 'data.frame'
as_draws_df(x, ...)

## S3 method for class 'draws_df'
as_draws_df(x, ...)

## S3 method for class 'draws_matrix'
as_draws_df(x, ...)

## S3 method for class 'draws_array'
as_draws_df(x, ...)

## S3 method for class 'draws_list'
as_draws_df(x, ...)

## S3 method for class 'draws_rvars'
as_draws_df(x, ...)

## S3 method for class 'mcmc'
as_draws_df(x, ...)

## S3 method for class 'mcmc.list'
as_draws_df(x, ...)

draws_df(..., .nchains = 1)

is_draws_df(x)

Arguments

Details

Objects of class "draws_df" are tibble data frames. They have one column per variable as well as additional metadata columns ".iteration", ".chain", and ".draw". The difference between the ".iteration" and ".draw" columns is that the former is relative to the MCMC chain while the latter ignores the chain information and has all unique values. See Examples.

If a data.frame-like object is supplied to as_draws_df that contains columns named ".iteration" or ".chain", they will be treated as iteration and chain indices, respectively. See Examples.

Value

A draws_df object, which has classes c("draws_df", "draws", class(tibble::tibble())).

See Also

Other formats: draws, draws_array(), draws_list(), draws_matrix(), draws_rvars()

Examples

x1 <- as_draws_df(example_draws())
class(x1)
print(x1)
str(x1)

x2 <- draws_df(a = rnorm(10), b = rnorm(10), c = 1)
class(x2)
print(x2)
str(x2)

# the difference between iteration and draw is clearer when contrasting
# the head and tail of the data frame
print(head(x1), reserved = TRUE, max_variables = 2)
print(tail(x1), reserved = TRUE, max_variables = 2)

# manually supply chain information
xnew <- data.frame(mu = rnorm(10), .chain = rep(1:2, each = 5))
xnew <- as_draws_df(xnew)
print(xnew)

The draws_list format

Description

The as_draws_list() methods convert objects to the draws_list format. The draws_list() function creates an object of the draws_list format based on a set of numeric vectors. See Details.

Usage

as_draws_list(x, ...)

## Default S3 method:
as_draws_list(x, ...)

## S3 method for class 'draws_list'
as_draws_list(x, ...)

## S3 method for class 'draws_matrix'
as_draws_list(x, ...)

## S3 method for class 'draws_array'
as_draws_list(x, ...)

## S3 method for class 'draws_df'
as_draws_list(x, ...)

## S3 method for class 'draws_rvars'
as_draws_list(x, ...)

## S3 method for class 'mcmc'
as_draws_list(x, ...)

## S3 method for class 'mcmc.list'
as_draws_list(x, ...)

draws_list(..., .nchains = 1)

is_draws_list(x)

Arguments

Details

Objects of class "draws_list" are lists with one element per MCMC chain. Each of these elements is itself a named list of numeric vectors with one vector per variable. The length of each vector is equal to the number of saved iterations per chain. See Examples.

Value

A draws_list object, which has classes c("draws_list", "draws", "list").

See Also

Other formats: draws, draws_array(), draws_df(), draws_matrix(), draws_rvars()

Examples

x1 <- as_draws_list(example_draws())
class(x1)
print(x1)
str(x1)

x2 <- draws_list(a = rnorm(10), b = rnorm(10), c = 1)
class(x2)
print(x2)
str(x2)

The draws_matrix format

Description

The as_draws_matrix() methods convert objects to the draws_matrix format. The draws_matrix() function creates an object of the draws_matrix format based on a set of numeric vectors. See Details.

Usage

as_draws_matrix(x, ...)

## Default S3 method:
as_draws_matrix(x, ...)

## S3 method for class 'draws_matrix'
as_draws_matrix(x, ...)

## S3 method for class 'draws_array'
as_draws_matrix(x, ...)

## S3 method for class 'draws_df'
as_draws_matrix(x, ...)

## S3 method for class 'draws_list'
as_draws_matrix(x, ...)

## S3 method for class 'draws_rvars'
as_draws_matrix(x, ...)

## S3 method for class 'mcmc'
as_draws_matrix(x, ...)

## S3 method for class 'mcmc.list'
as_draws_matrix(x, ...)

draws_matrix(..., .nchains = 1)

is_draws_matrix(x)

Arguments

Details

Objects of class "draws_matrix" are matrices (2-D arrays) with dimensions "draw" and "variable". See Examples.

Value

A draws_matrix object, which has classes c("draws_matrix", "draws", "matrix").

See Also

Other formats: draws, draws_array(), draws_df(), draws_list(), draws_rvars()

Examples

x1 <- as_draws_matrix(example_draws())
class(x1)
print(x1)
str(x1)

x2 <- draws_matrix(a = rnorm(10), b = rnorm(10), c = 1)
class(x2)
print(x2)
str(x2)

Get/set array of draws underlying a random variable

Description

Gets/sets the array-representation that backs an rvar. Should be used rarely.

Usage

draws_of(x, with_chains = FALSE)

draws_of(x, with_chains = FALSE) <- value

Arguments

Details

While rvars implement fast versions of basic math operations (including matrix multiplication), sometimes you may need to bypass the rvar abstraction to do what you need to do more efficiently. draws_of() allows you to get / set the underlying array of draws in order to do that.

rvars represent draws internally using arrays of arbitrary dimension, which is returned by draws_of(x) and can be set using draws_of(x) <- value. The first dimension of these arrays is the index of the draws. If with_chains = TRUE, then the dimensions of the returned array are modified so that the first dimension is the index of the iterations and the second dimension is the index of the chains.

Value

If with_chains = FALSE, an array with dimensions c(ndraws(x), dim(x)).

If with_chains = TRUE, an array with dimensions c(niterations(x), nchains(x), dim(x)).

Examples

x <- rvar(1:10, nchains = 2)
x

# draws_of() without arguments will return the array of draws without
# chain information (first dimension is draw)
draws_of(x)

# draws_of() with with_chains = TRUE will reshape the returned array to
# include chain information in the second dimension
draws_of(x, with_chains = TRUE)

# you can also set draws using draws_of(). When with_chains = FALSE the
# existing chain information will be retained ...
draws_of(x) <- 2:11
x

# when with_chains = TRUE the chain information will be set by the
# second dimension of the assigned array
draws_of(x, with_chains = TRUE) <- array(2:11, dim = c(2,5))
x

The draws_rvars format

Description

The as_draws_rvars() methods convert objects to the draws_rvars format. The draws_rvars() function creates an object of the draws_rvars format based on a set of numeric vectors. See Details.

Usage

as_draws_rvars(x, ...)

## Default S3 method:
as_draws_rvars(x, ...)

## S3 method for class 'draws_rvars'
as_draws_rvars(x, ...)

## S3 method for class 'list'
as_draws_rvars(x, ...)

## S3 method for class 'draws_matrix'
as_draws_rvars(x, ...)

## S3 method for class 'draws_array'
as_draws_rvars(x, ...)

## S3 method for class 'draws_df'
as_draws_rvars(x, ...)

## S3 method for class 'draws_list'
as_draws_rvars(x, ...)

## S3 method for class 'mcmc'
as_draws_rvars(x, ...)

## S3 method for class 'mcmc.list'
as_draws_rvars(x, ...)

draws_rvars(..., .nchains = 1)

is_draws_rvars(x)

Arguments

Details

Objects of class "draws_rvars" are lists of rvar objects. See Examples.

Value

A draws_rvars object, which has classes c("draws_rvars", "draws", "list").

See Also

Other formats: draws, draws_array(), draws_df(), draws_list(), draws_matrix()

Examples

x1 <- as_draws_rvars(example_draws())
class(x1)
print(x1)
str(x1)

x2 <- draws_rvars(a = rnorm(10), b = rnorm(10), c = 1)
class(x2)
print(x2)
str(x2)

Summaries of draws objects

Description

The summarise_draws() (and summarize_draws()) methods provide a quick way to get a table of summary statistics and diagnostics. These methods will convert an object to a draws object if it isn't already. For convenience, a summary() method for draws and rvar objects are also provided as an alias for summarise_draws() if the input object is a draws or rvar object.

Usage

summarise_draws(.x, ...)

summarize_draws(.x, ...)

## S3 method for class 'draws'
summarise_draws(
  .x,
  ...,
  .args = list(),
  .num_args = getOption("posterior.num_args", list()),
  .cores = 1
)

## S3 method for class 'draws'
summary(object, ...)

## S3 method for class 'rvar'
summarise_draws(.x, ...)

## S3 method for class 'rvar'
summary(object, ...)

default_summary_measures()

default_convergence_measures()

default_mcse_measures()

Arguments

Details

The default summary functions used are the ones specified by default_summary_measures() and default_convergence_measures():

default_summary_measures()

• mean()

• median()

• sd()

• mad()

• quantile2()

default_convergence_measures()

• rhat()

• ess_bulk()

• ess_tail()

The var() function should not be used to compute variances due to its inconsistent behavior with matrices. Instead, please use distributional::variance().

Value

The summarise_draws() methods return a tibble data frame. The first column ("variable") contains the variable names and the remaining columns contain summary statistics and diagnostics.

The functions default_summary_measures(), default_convergence_measures(), and default_mcse_measures() return character vectors of names of the default measures.

See Also

diagnostics for a list of available diagnostics and links to their individual help pages.

Examples

x <- example_draws("eight_schools")
class(x)
str(x)

summarise_draws(x)
summarise_draws(x, "mean", "median")
summarise_draws(x, mean, mcse = mcse_mean)
summarise_draws(x, ~quantile(.x, probs = c(0.4, 0.6)))

# using default_*_meaures()
summarise_draws(x, default_summary_measures())
summarise_draws(x, default_convergence_measures())
summarise_draws(x, default_mcse_measures())

# compute variance of variables
summarise_draws(x, var = distributional::variance)

# illustrate use of '.args'
ws <- rexp(ndraws(x))
summarise_draws(x, weighted.mean, .args = list(w = ws))

# adjust how numerical summaries are printed
summarise_draws(x, .num_args = list(sigfig = 2, notation = "dec"))

Index draws objects

Description

Index iterations, chains, and draws of draws objects.

Usage

iteration_ids(x)

chain_ids(x)

draw_ids(x)

niterations(x)

nchains(x)

ndraws(x)

Arguments

Details

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.

Value

For iteration_ids(), chain_ids(), and draw_ids(), an integer vector.

For niterations(), nchains(), and ndraws(), a scalar integer.

See Also

variables, rename_variables

Examples

x <- example_draws()

iteration_ids(x)
niterations(x)

chain_ids(x)
nchains(x)

draw_ids(x)
ndraws(x)

Drop redundant dimensions

Description

Delete the dimensions of an rvar which are of size one. See base::drop()

Usage

## S4 method for signature 'rvar'
drop(x)

Arguments

Value

An rvar with the same length as x, but where any entry equal to 1 in dim(x) has been removed. The exception is if dim(x) == 1, in which case dim(drop(x)) == 1 as well (this is because rvars, unlike numerics, never have NULL dimensions).

Examples

# Sigma is a 3x3 covariance matrix
Sigma <- as_draws_rvars(example_draws("multi_normal"))$Sigma
Sigma

Sigma[1, ]

drop(Sigma[1, ])

# equivalently ...
Sigma[1, drop = TRUE]

Normalized entropy

Description

Normalized entropy, for measuring dispersion in draws from categorical distributions.

Usage

entropy(x)

## Default S3 method:
entropy(x)

## S3 method for class 'rvar'
entropy(x)

Arguments

Details

Calculates the normalized Shannon entropy of the draws in x. This value is the entropy of x divided by the maximum entropy of a distribution with n categories, where n is length(unique(x)) for numeric vectors and length(levels(x)) for factors:

$-\frac{\sum_{i = 1}^{n} p_i \log(p_i)}{\log(n)}$

This scales the output to be between 0 (all probability in one category) and 1 (uniform). This form of normalized entropy is referred to as $H_\mathrm{REL}$ in Wilcox (1967).

Value

If x is a factor or numeric, returns a length-1 numeric vector with a value between 0 and 1 (inclusive) giving the normalized Shannon entropy of x.

If x is an rvar, returns an array of the same shape as x, where each cell is the normalized Shannon entropy of the draws in the corresponding cell of x.

References

Allen R. Wilcox (1967). Indices of Qualitative Variation (No. ORNL-TM-1919). Oak Ridge National Lab., Tenn.

Examples

set.seed(1234)

levels <- c("a", "b", "c", "d", "e")

# a uniform distribution: high normalized entropy
x <- factor(
  sample(levels, 4000, replace = TRUE, prob = c(0.2, 0.2, 0.2, 0.2, 0.2)),
  levels = levels
)
entropy(x)

# a unimodal distribution: low normalized entropy
y <- factor(
  sample(levels, 4000, replace = TRUE, prob = c(0.95, 0.02, 0.015, 0.01, 0.005)),
  levels = levels
)
entropy(y)

# both together, as an rvar
xy <- c(rvar(x), rvar(y))
xy
entropy(xy)

Basic version of the effective sample size

Description

Compute the basic effective sample size (ESS) estimate for a single variable as described in Gelman et al. (2013) with some changes according to Vehtari et al. (2021). For practical applications, we strongly recommend the improved ESS convergence diagnostics implemented in ess_bulk() and ess_tail(). See Vehtari (2021) for an in-depth comparison of different effective sample size estimators.

Usage

ess_basic(x, ...)

## Default S3 method:
ess_basic(x, split = TRUE, ...)

## S3 method for class 'rvar'
ess_basic(x, split = TRUE, ...)

Arguments

Value

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

Andrew Gelman, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari and Donald B. Rubin (2013). Bayesian Data Analysis, Third Edition. Chapman and 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). Bayesian Analysis. 16(2), 667-–718. doi:10.1214/20-BA1221

Aki Vehtari (2021). Comparison of MCMC effective sample size estimators. Retrieved from https://avehtari.github.io/rhat_ess/ess_comparison.html

See Also

Other diagnostics: ess_bulk(), ess_quantile(), ess_sd(), ess_tail(), mcse_mean(), mcse_quantile(), mcse_sd(), pareto_diags(), pareto_khat(), rhat(), rhat_basic(), rhat_nested(), rstar()

Examples

mu <- extract_variable_matrix(example_draws(), "mu")
ess_basic(mu)

d <- as_draws_rvars(example_draws("multi_normal"))
ess_basic(d$Sigma)

Bulk effective sample size (bulk-ESS)

Description

Compute a bulk effective sample size estimate (bulk-ESS) for a single variable. Bulk-ESS is useful as a diagnostic for the sampling efficiency in 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 ess_tail(). See Vehtari (2021) for an in-depth comparison of different effective sample size estimators.

Usage

ess_bulk(x, ...)

## Default S3 method:
ess_bulk(x, ...)

## S3 method for class 'rvar'
ess_bulk(x, ...)

Arguments

Value

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, 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 Analysis. 16(2), 667-–718. doi:10.1214/20-BA1221

Aki Vehtari (2021). Comparison of MCMC effective sample size estimators. Retrieved from https://avehtari.github.io/rhat_ess/ess_comparison.html

See Also

Other diagnostics: ess_basic(), ess_quantile(), ess_sd(), ess_tail(), mcse_mean(), mcse_quantile(), mcse_sd(), pareto_diags(), pareto_khat(), rhat(), rhat_basic(), rhat_nested(), rstar()

Examples

mu <- extract_variable_matrix(example_draws(), "mu")
ess_bulk(mu)

d <- as_draws_rvars(example_draws("multi_normal"))
ess_bulk(d$Sigma)

Effective sample size for the mean

Description

Compute an effective sample size estimate for a mean (expectation) estimate of a single variable.

Usage

ess_mean(x, ...)

## Default S3 method:
ess_mean(x, ...)

## S3 method for class 'rvar'
ess_mean(x, ...)

Arguments

Value

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

Andrew Gelman, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari and Donald B. Rubin (2013). Bayesian Data Analysis, Third Edition. Chapman and Hall/CRC.

Examples

mu <- extract_variable_matrix(example_draws(), "mu")
ess_mean(mu)

d <- as_draws_rvars(example_draws("multi_normal"))
ess_mean(d$Sigma)

Effective sample sizes for quantiles

Description

Compute effective sample size estimates for quantile estimates of a single variable.

Usage

ess_quantile(x, probs = c(0.05, 0.95), ...)

## Default S3 method:
ess_quantile(x, probs = c(0.05, 0.95), names = TRUE, ...)

## S3 method for class 'rvar'
ess_quantile(x, probs = c(0.05, 0.95), names = TRUE, ...)

ess_median(x, ...)

Arguments

Value

If the input is an array, returns a numeric vector with one element per quantile. If any of the draws is non-finite, that is, NA, NaN, Inf, or -Inf, the returned output will be a vector of (numeric) NA values. Also, if all draws of a variable are the same (constant), the returned output will be a vector of (numeric) NA values 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 and length(probs) == 1, 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. If length(probs) > 1, the first dimension of the result indexes the input probabilities; i.e. the result has dimension c(length(probs), dim(x)).

References

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 Analysis. 16(2), 667-–718. doi:10.1214/20-BA1221

See Also

Other diagnostics: ess_basic(), ess_bulk(), ess_sd(), ess_tail(), mcse_mean(), mcse_quantile(), mcse_sd(), pareto_diags(), pareto_khat(), rhat(), rhat_basic(), rhat_nested(), rstar()

Examples

mu <- extract_variable_matrix(example_draws(), "mu")
ess_quantile(mu, probs = c(0.1, 0.9))

d <- as_draws_rvars(example_draws("multi_normal"))
ess_quantile(d$mu, probs = c(0.1, 0.9))

Effective sample size for the standard deviation

Description

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.

Usage

ess_sd(x, ...)

## Default S3 method:
ess_sd(x, ...)

## S3 method for class 'rvar'
ess_sd(x, ...)

Arguments

Value

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, 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 Analysis. 16(2), 667-–718. doi:10.1214/20-BA1221

See Also

Other diagnostics: ess_basic(), ess_bulk(), ess_quantile(), ess_tail(), mcse_mean(), mcse_quantile(), mcse_sd(), pareto_diags(), pareto_khat(), rhat(), rhat_basic(), rhat_nested(), rstar()

Examples

mu <- extract_variable_matrix(example_draws(), "mu")
ess_sd(mu)

d <- as_draws_rvars(example_draws("multi_normal"))
ess_sd(d$Sigma)

Tail effective sample size (tail-ESS)

Description

Compute a tail effective sample size estimate (tail-ESS) for a single variable. Tail-ESS is useful as a diagnostic for the sampling efficiency in 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 ess_bulk(). See Vehtari (2021) for an in-depth comparison of different effective sample size estimators.

Usage

ess_tail(x, ...)

## Default S3 method:
ess_tail(x, ...)

## S3 method for class 'rvar'
ess_tail(x, ...)

Arguments

Value

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, 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 Analysis. 16(2), 667-–718. doi:10.1214/20-BA1221

Aki Vehtari (2021). Comparison of MCMC effective sample size estimators. Retrieved from https://avehtari.github.io/rhat_ess/ess_comparison.html

See Also

Other diagnostics: ess_basic(), ess_bulk(), ess_quantile(), ess_sd(), mcse_mean(), mcse_quantile(), mcse_sd(), pareto_diags(), pareto_khat(), rhat(), rhat_basic(), rhat_nested(), rstar()

Examples

mu <- extract_variable_matrix(example_draws(), "mu")
ess_tail(mu)

d <- as_draws_rvars(example_draws("multi_normal"))
ess_tail(d$Sigma)

Example draws objects

Description

Objects for use in examples, vignettes, and tests.

Usage

example_draws(example = "eight_schools")

Arguments

Details

The following example draws objects are available.

eight_schools: A draws_array object with 100 iterations from each of 4 Markov chains obtained by fitting the eight schools model described in Gelman et al. (2013) with Stan. The variables are:

• mu: Overall mean of the eight schools

• tau: Standard deviation between schools

• theta: Individual means of each of the eight schools

multi_normal: A draws_array object with 100 iterations from each of the 4 Markov chains obtained by fitting a 3-dimensional multivariate normal model to 100 simulated observations. The variables are:

• mu: Mean parameter vector of length 3

• Sigma: Covariance matrix of dimension 3 x 3

Value

A draws object.

Note

These objects are only intended to be used in demonstrations and tests. They contain fewer iterations and chains than recommended for performing actual inference.

References

Andrew Gelman, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari and Donald B. Rubin (2013). Bayesian Data Analysis, Third Edition. Chapman and Hall/CRC.

Examples

draws_eight_schools <- example_draws("eight_schools")
summarise_draws(draws_eight_schools)

draws_multi_normal <- example_draws("multi_normal")
summarise_draws(draws_multi_normal)

Extract draws of a single variable

Description

Extract a vector of draws of a single variable.

Usage

extract_variable(x, variable, ...)

## Default S3 method:
extract_variable(x, variable, ...)

## S3 method for class 'draws'
extract_variable(x, variable, ...)

## S3 method for class 'draws_df'
extract_variable(x, variable, ...)

## S3 method for class 'draws_list'
extract_variable(x, variable, ...)

## S3 method for class 'draws_rvars'
extract_variable(x, variable, ...)

Arguments

Value

A vector of length equal to the number of draws.

See Also

Other variable extraction methods: extract_variable_array(), extract_variable_matrix()

Examples

x <- example_draws()
mu <- extract_variable(x, variable = "mu")
str(mu)

Extract array of a single (possibly indexed) variable

Description

Extract an array of draws of a single variable, including any dimensions of variables with indices.

Usage

extract_variable_array(x, variable, ...)

## Default S3 method:
extract_variable_array(x, variable, ...)

## S3 method for class 'draws'
extract_variable_array(x, variable, ...)

Arguments

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)

mu1 <- extract_variable_array(x, variable = "mu[1]")
str(mu1)

Sigma <- extract_variable_array(x, variable = "Sigma")
str(Sigma)

Extract matrix of a single variable

Description

Extract an iterations x chains matrix of draws of a single variable. This is primarily used for convergence diagnostic functions such as rhat().

Usage

extract_variable_matrix(x, variable, ...)

## Default S3 method:
extract_variable_matrix(x, variable, ...)

## S3 method for class 'draws'
extract_variable_matrix(x, variable, ...)

## S3 method for class 'draws_df'
extract_variable_matrix(x, variable, ...)

## S3 method for class 'draws_list'
extract_variable_matrix(x, variable, ...)

## S3 method for class 'draws_rvars'
extract_variable_matrix(x, variable, ...)

Arguments

Value

A matrix with dimension iterations x chains.

See Also

Other variable extraction methods: extract_variable(), extract_variable_array()

Examples

x <- example_draws()
mu <- extract_variable_matrix(x, variable = "mu")
dim(mu)
rhat(mu)

Loop over draws

Description

Executes an expression once for every draw in a draws object. Used primarily for its side effects and returns the input x invisibly.

Usage

for_each_draw(x, expr)

Arguments

Details

If x is not in the draws_rvars format, it is first converted to that format. This allows the variables in x to include their dimensions (i.e, to act as R vectors and arrays) when being referred to in expr.

Within expr, use .draw to refer to the draw index, which will be a value between 1 and ndraws(x). expr is executed in the calling environment of for_each_draw(), so it can use variables in that environment (however, due to the use of data masking, to modify variables in that environment, one must use ⁠<<-⁠.)

Value

As for_each_draw() is used primarily for its side effects (the expression executed for each draw of x), it returns the input x invisibly.

Examples

eight_schools <- as_draws_rvars(example_draws())


# 1. A simple example --- looping over draws and printing each draw
# NOTE: You probably don't want to do this in practice! This example is
# just intended to show what for_each_draw() is doing. If you just want to
# print the draws of an rvar, it is probably better to use draws_of()
for_each_draw(eight_schools, {
  print(mu)
})


# 2. A more complex example --- building a parallel coordinates plot
# First, construct the plot bounds
plot(1, type = "n",
  xlim = c(1, length(eight_schools$theta)),
  ylim = range(range(eight_schools$theta)),
  xlab = "school", ylab = "theta"
)

# Then, use for_each_draw() to make a parallel coordinates plot of all draws
# of eight_schools$theta. Use resample_draws(eight_schools, n = ...)
# in place of eight_schools if a smaller sample is desired for the plot.
for_each_draw(eight_schools, {
  lines(seq_along(theta), theta, col = rgb(1, 0, 0, 0.05))
})

# Finally, add means and 90% intervals
lines(seq_along(eight_schools$theta), mean(eight_schools$theta))
with(summarise_draws(eight_schools$theta),
  segments(seq_along(eight_schools$theta), y0 = q5, y1 = q95)
)

Is x a random variable?

Description

Test if x is an rvar.

Usage

is_rvar(x)

Arguments

Value

TRUE if x is an rvar, FALSE otherwise.

See Also

as_rvar() to convert objects to rvars.

Is x a factor random variable?

Description

Test if x is an rvar_factor or rvar_ordered.

Usage

is_rvar_factor(x)

is_rvar_ordered(x)

Arguments

Value

TRUE if x is an rvar_factor or rvar_ordered, FALSE otherwise.

See Also

as_rvar_factor() and as_rvar_ordered() to convert objects to rvar_factors and rvar_ordereds.

Value Matching

Description

Generic version of base::match(). For base vectors, returns a vector of the positions of (first) matches of its first argument in its second. For rvars, returns an rvar of the matches.

Usage

match(x, table, ...)

## Default S3 method:
match(x, ...)

## S3 method for class 'rvar'
match(x, ...)

x %in% table

Arguments

Details

For more information on how match behaves with base vectors, see base::match().

When x is an rvar, the draws of x are matched against table using base::match(), and the result is returned as an rvar.

The implementation of %in% here is identical to base::%in%, except it uses the generic version of match() so that non-base vectors (such as rvars) are supported.

Value

When x is a base vector, a vector of the same length as x.

When x is an rvar, an rvar the same shape as x.

Examples

x <- rvar(c("a","b","b","c","d"))
x %in% c("b","d")

# for additional examples, see base::match()

Monte Carlo standard error for the mean

Description

Compute the Monte Carlo standard error for the mean (expectation) of a single variable.

Usage

mcse_mean(x, ...)

## Default S3 method:
mcse_mean(x, ...)

## S3 method for class 'rvar'
mcse_mean(x, ...)

Arguments

Value

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

Andrew Gelman, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari and Donald B. Rubin (2013). Bayesian Data Analysis, Third Edition. Chapman and Hall/CRC.

See Also

Other diagnostics: ess_basic(), ess_bulk(), ess_quantile(), ess_sd(), ess_tail(), mcse_quantile(), mcse_sd(), pareto_diags(), pareto_khat(), rhat(), rhat_basic(), rhat_nested(), rstar()

Examples

mu <- extract_variable_matrix(example_draws(), "mu")
mcse_mean(mu)

d <- as_draws_rvars(example_draws("multi_normal"))
mcse_mean(d$Sigma)

Monte Carlo standard error for quantiles

Description

Compute Monte Carlo standard errors for quantile estimates of a single variable.

Usage

mcse_quantile(x, probs = c(0.05, 0.95), ...)

## Default S3 method:
mcse_quantile(x, probs = c(0.05, 0.95), names = TRUE, ...)

## S3 method for class 'rvar'
mcse_quantile(x, probs = c(0.05, 0.95), names = TRUE, ...)

mcse_median(x, ...)

Arguments

Value

If the input is an array, returns a numeric vector with one element per quantile. If any of the draws is non-finite, that is, NA, NaN, Inf, or -Inf, the returned output will be a vector of (numeric) NA values. Also, if all draws of a variable are the same (constant), the returned output will be a vector of (numeric) NA values 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 and length(probs) == 1, 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. If length(probs) > 1, the first dimension of the result indexes the input probabilities; i.e. the result has dimension c(length(probs), dim(x)).

References

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 Analysis. 16(2), 667-–718. doi:10.1214/20-BA1221

See Also

Other diagnostics: ess_basic(), ess_bulk(), ess_quantile(), ess_sd(), ess_tail(), mcse_mean(), mcse_sd(), pareto_diags(), pareto_khat(), rhat(), rhat_basic(), rhat_nested(), rstar()

Examples

mu <- extract_variable_matrix(example_draws(), "mu")
mcse_quantile(mu, probs = c(0.1, 0.9))

d <- as_draws_rvars(example_draws("multi_normal"))
mcse_quantile(d$mu)

Monte Carlo standard error for the standard deviation

Description

Compute the Monte Carlo standard error for the standard deviation (SD) of a single variable without assuming normality using moments of moments and first order Taylor series approximation (Kenney and Keeping, 1951, p. 141).

Usage

mcse_sd(x, ...)

## Default S3 method:
mcse_sd(x, ...)

## S3 method for class 'rvar'
mcse_sd(x, ...)

Arguments

Value

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, 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 Analysis. 16(2), 667-–718. doi:10.1214/20-BA1221

J. F. Kenney & E. S. Keeping (1951). Mathematics of Statistics, Vol. II.

See Also

Other diagnostics: ess_basic(), ess_bulk(), ess_quantile(), ess_sd(), ess_tail(), mcse_mean(), mcse_quantile(), pareto_diags(), pareto_khat(), rhat(), rhat_basic(), rhat_nested(), rstar()

Examples

mu <- extract_variable_matrix(example_draws(), "mu")
mcse_sd(mu)

d <- as_draws_rvars(example_draws("multi_normal"))
mcse_sd(d$Sigma)

Merge chains of draws objects

Description

Merge chains of draws objects into a single chain. Some operations will trigger an automatic merging of chains, for example, because chains do not match between two objects involved in a binary operation. By default, no warning will be issued when this happens but you can activate one via options(posterior.warn_on_merge_chains = TRUE).

Usage

merge_chains(x, ...)

## S3 method for class 'draws_matrix'
merge_chains(x, ...)

## S3 method for class 'draws_array'
merge_chains(x, ...)

## S3 method for class 'draws_df'
merge_chains(x, ...)

## S3 method for class 'draws_list'
merge_chains(x, ...)

## S3 method for class 'rvar'
merge_chains(x, ...)

## S3 method for class 'draws_rvars'
merge_chains(x, ...)

Arguments

Value

A draws object of the same class as x.

Examples

x <- example_draws()

# draws_array with 4 chains, 100 iters each
str(x)

# draws_array with 1 chain of 400 iterations
str(merge_chains(x))

Modal category

Description

Modal category of a vector.

Usage

modal_category(x)

## Default S3 method:
modal_category(x)

## S3 method for class 'rvar'
modal_category(x)

Arguments

Details

Finds the modal category (i.e., most frequent value) in x. In the case of ties, returns the first tie.

Value

If x is a factor or numeric, returns a length-1 vector containing the modal value.

If x is an rvar, returns an array of the same shape as x, where each cell is the modal value of the draws in the corresponding cell of x.

Examples

x <- factor(c("a","b","b","c","d"))
modal_category(x)

# in the case of ties, the first tie is returned
y <- factor(c("a","c","c","d","d"))
modal_category(y)

# both together, as an rvar
xy <- c(rvar(x), rvar(y))
xy
modal_category(xy)

Mutate variables in draws objects

Description

Mutate variables in a draws object.

Usage

mutate_variables(.x, ...)

## S3 method for class 'draws_matrix'
mutate_variables(.x, ...)

## S3 method for class 'draws_array'
mutate_variables(.x, ...)

## S3 method for class 'draws_df'
mutate_variables(.x, ...)

## S3 method for class 'draws_list'
mutate_variables(.x, ...)

## S3 method for class 'draws_rvars'
mutate_variables(.x, ...)

Arguments

Details

In order to mutate variables in draws_matrix and draws_array objects, they are transformed to draws_df objects first and then transformed back after mutation. As those transformations are quite expensive for larger number of draws, we recommend using mutate_variables on draws_df and draws_list objects if speed is an issue.

In draws_rvars objects, the output of each expression in ... is coerced to an rvar object if it is not already one using as_rvar().

Value

Returns a draws object of the same format as .x, with variables mutated according to the expressions provided in ....

See Also

variables, rename_variables

Examples

x <- as_draws_df(example_draws())
x <- subset(x, variable = c("mu", "tau"))

mutate_variables(x, tau2 = tau^2)
mutate_variables(x, scale = 1.96 * tau, lower = mu - scale)

Order draws objects

Description

Order draws objects according to iteration and chain number. By default, draws objects are ordered but subsetting or extracting parts of them may leave them in an unordered state.

Usage

order_draws(x, ...)

## S3 method for class 'draws_matrix'
order_draws(x, ...)

## S3 method for class 'draws_array'
order_draws(x, ...)

## S3 method for class 'draws_df'
order_draws(x, ...)

## S3 method for class 'draws_list'
order_draws(x, ...)

## S3 method for class 'draws_rvars'
order_draws(x, ...)

## S3 method for class 'rvar'
order_draws(x, ...)

Arguments

Value

A draws object of the same class as x.

See Also

repair_draws()

Examples

x <- as_draws_array(example_draws())
dimnames(x[10:5, 4:3, ])
dimnames(order_draws(x[10:5, 4:3, ]))

Pareto smoothing diagnostics

Description

Compute diagnostics for Pareto smoothing the tail draws of x by replacing tail draws by order statistics of a generalized Pareto distribution fit to the tail(s).

Usage

pareto_diags(x, ...)

## Default S3 method:
pareto_diags(
  x,
  tail = c("both", "right", "left"),
  r_eff = NULL,
  ndraws_tail = NULL,
  verbose = FALSE,
  are_log_weights = FALSE,
  ...
)

## S3 method for class 'rvar'
pareto_diags(x, ...)

pareto_khat_threshold(x, ...)

## Default S3 method:
pareto_khat_threshold(x, ...)

## S3 method for class 'rvar'
pareto_khat_threshold(x, ...)

pareto_min_ss(x, ...)

## Default S3 method:
pareto_min_ss(x, ...)

## S3 method for class 'rvar'
pareto_min_ss(x, ...)

pareto_convergence_rate(x, ...)

## Default S3 method:
pareto_convergence_rate(x, ...)

## S3 method for class 'rvar'
pareto_convergence_rate(x, ...)

Arguments

Details

When the fitted Generalized Pareto Distribution is used to smooth the tail values and these smoothed values are used to compute expectations, the following diagnostics can give further information about the reliability of these estimates.

• min_ss: Minimum sample size for reliable Pareto smoothed estimate. If the actual sample size is greater than min_ss, then Pareto smoothed estimates can be considered reliable. If the actual sample size is lower than min_ss, increasing the sample size might result in more reliable estimates. For further details, see Section 3.2.3, Equation 11 in Vehtari et al. (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. (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. (2024).

Value

List of Pareto smoothing diagnostics:

• khat: estimated Pareto k shape parameter,

• min_ss: minimum sample size for reliable Pareto smoothed estimate,

• khat_threshold: khat-threshold for reliable Pareto smoothed estimate,

• convergence_rate: Pareto smoothed estimate RMSE convergence rate.

References

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. 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()

Examples

mu <- extract_variable_matrix(example_draws(), "mu")
pareto_diags(mu)

d <- as_draws_rvars(example_draws("multi_normal"))
pareto_diags(d$Sigma)

Pareto khat diagnostic

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. (2024).

Usage

pareto_khat(x, ...)

## Default S3 method:
pareto_khat(
  x,
  tail = c("both", "right", "left"),
  r_eff = NULL,
  ndraws_tail = NULL,
  verbose = FALSE,
  are_log_weights = FALSE,
  ...
)

## S3 method for class 'rvar'
pareto_khat(x, ...)

Arguments

Value

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 (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)

d <- as_draws_rvars(example_draws("multi_normal"))
pareto_khat(d$Sigma)

Pareto smoothing

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. (2024).

Usage

pareto_smooth(x, ...)

## S3 method for class 'rvar'
pareto_smooth(x, return_k = FALSE, extra_diags = FALSE, ...)

## Default S3 method:
pareto_smooth(
  x,
  tail = c("both", "right", "left"),
  r_eff = NULL,
  ndraws_tail = NULL,
  return_k = FALSE,
  extra_diags = FALSE,
  verbose = TRUE,
  are_log_weights = FALSE,
  ...
)

Arguments

Value

Either a vector x of smoothed values or a named list 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

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 (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)

d <- as_draws_rvars(example_draws("multi_normal"))
pareto_smooth(d$Sigma)

Print draws_array objects

Description

Pretty printing for draws_array objects.

Usage

## S3 method for class 'draws_array'
print(
  x,
  digits = 2,
  max_iterations = getOption("posterior.max_iterations", 5),
  max_chains = getOption("posterior.max_chains", 8),
  max_variables = getOption("posterior.max_variables", 4),
  reserved = FALSE,
  ...
)

Arguments

Value

A draws object of the same class as x.

Examples

x <- as_draws_array(example_draws())
print(x)

Print draws_df objects

Description

Pretty printing for draws_df objects.

Usage

## S3 method for class 'draws_df'
print(
  x,
  digits = 2,
  max_draws = getOption("posterior.max_draws", 10),
  max_variables = getOption("posterior.max_variables", 8),
  reserved = FALSE,
  ...
)

Arguments

Value

A draws object of the same class as x.

Examples

x <- as_draws_df(example_draws())
print(x)

Print draws_list objects

Description

Pretty printing for draws_list objects.

Usage

## S3 method for class 'draws_list'
print(
  x,
  digits = 2,
  max_iterations = getOption("posterior.max_iterations", 10),
  max_chains = getOption("posterior.max_chains", 2),
  max_variables = getOption("posterior.max_variables", 4),
  reserved = FALSE,
  ...
)

Arguments

Value

A draws object of the same class as x.

Examples

x <- as_draws_list(example_draws())
print(x)

Print draws_matrix objects

Description

Pretty printing for draws_matrix objects.

Usage

## S3 method for class 'draws_matrix'
print(
  x,
  digits = 2,
  max_draws = getOption("posterior.max_draws", 10),
  max_variables = getOption("posterior.max_variables", 8),
  reserved = FALSE,
  ...
)

Arguments

Value

A draws object of the same class as x.

Examples

x <- as_draws_matrix(example_draws())
print(x)

Print draws_rvars objects

Description

Pretty printing for draws_rvars objects.

Usage

## S3 method for class 'draws_rvars'
print(
  x,
  digits = 2,
  max_variables = getOption("posterior.max_variables", 8),
  summary = getOption("posterior.rvar_summary", "mean_sd"),
  reserved = FALSE,
  ...
)

Arguments

Value

A draws object of the same class as x.

Examples

x <- as_draws_rvars(example_draws())
print(x)

Print summaries of draws objects

Description

Print output from summarise_draws().

Usage

## S3 method for class 'draws_summary'
print(x, ..., num_args = NULL)

Arguments

Value

An invisible version of the input object.

Examples

x <- example_draws("eight_schools")

# adjust how summaries are printed when calling summarise_draws()...
summarise_draws(x, .num_args = list(sigfig = 2, notation = "dec"))

# ... or when printing
s <- summarise_draws(x)
print(s, num_args = list(sigfig = 2, notation = "dec"))
print(s, num_args = list(digits = 3))

Print or format a random variable

Description

Printing and formatting methods for rvars.

Usage

## S3 method for class 'rvar'
print(
  x,
  ...,
  summary = NULL,
  digits = NULL,
  color = TRUE,
  width = getOption("width")
)

## S3 method for class 'rvar'
format(x, ..., summary = NULL, digits = NULL, color = FALSE)

## S3 method for class 'rvar'
str(
  object,
  ...,
  summary = NULL,
  vec.len = NULL,
  indent.str = paste(rep.int(" ", max(0, nest.lev + 1)), collapse = ".."),
  nest.lev = 0,
  give.attr = TRUE
)

Arguments

Details

print() and str() print out rvar objects by summarizing each element in the random variable with either its mean±sd or median±mad, depending on the value of summary. Both functions use the format() implementation for rvar objects under the hood, which returns a character vector in the mean±sd or median±mad form.

Value

For print(), an invisible version of the input object.

For str(), nothing; i.e. invisible(NULL).

For format(), a character vector of the same dimensions as x where each entry is of the form "mean±sd" or "median±mad", depending on the value of summary.

References

William J. Tastle, Mark J. Wierman (2007). Consensus and dissention: A measure of ordinal dispersion. International Journal of Approximate Reasoning. 45(3), 531–545. doi:10.1016/j.ijar.2006.06.024.

Examples

set.seed(5678)
x = rbind(
  cbind(rvar(rnorm(1000, 1)), rvar(rnorm(1000, 2))),
  cbind(rvar(rnorm(1000, 3)), rvar(rnorm(1000, 4)))
)

print(x)
print(x, summary = "median_mad")

str(x)

format(x)

Pareto 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 pareto_convergence_rate and pareto_diags.

Usage

ps_convergence_rate(k, ndraws, ...)

Arguments

Value

convergence rate

See Also

Other helper-functions: ps_khat_threshold(), ps_min_ss(), ps_tail_length()

Pareto k-hat 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 pareto_khat_threshold and pareto_diags.

Usage

ps_khat_threshold(ndraws, ...)

Arguments

Value

threshold

See Also

Other helper-functions: ps_convergence_rate(), ps_min_ss(), ps_tail_length()

Pareto-smoothing 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 pareto_min_ss and pareto_diags.

Usage

ps_min_ss(k, ...)

Arguments

Value

minimum sample size

See Also

Other helper-functions: ps_convergence_rate(), ps_khat_threshold(), ps_tail_length()

Pareto 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.

Usage

ps_tail_length(ndraws, r_eff, ...)

Arguments

Value

tail length

See Also

Other helper-functions: ps_convergence_rate(), ps_khat_threshold(), ps_min_ss()

Compute Quantiles

Description

Compute quantiles of a sample and return them in a format consistent with other summary functions in the posterior package.

Usage

quantile2(x, probs = c(0.05, 0.95), na.rm = FALSE, ...)

## Default S3 method:
quantile2(x, probs = c(0.05, 0.95), na.rm = FALSE, names = TRUE, ...)

## S3 method for class 'rvar'
quantile2(x, probs = c(0.05, 0.95), na.rm = FALSE, names = TRUE, ...)

Arguments

Value

A numeric vector of length length(probs). If names = TRUE, it has a names attribute with names like "q5", "q95", etc, based on the values of probs.

Examples

mu <- extract_variable_matrix(example_draws(), "mu")
quantile2(mu)

Rank values

Description

Compute ranks for a numeric array, that is, replace each value by its rank. Average rank for ties are used to conserve the number of unique values of discrete quantities.

Usage

r_scale(x)

Arguments

Value

A numeric array of ranked values with the same size and dimension as the input.

Execute expressions of random variables

Description

Execute (nearly) arbitrary R expressions that may include rvars, producing a new rvar.

Usage

rdo(expr, dim = NULL, ndraws = NULL)

Arguments

Details

This function evaluates expr possibly multiple times, once for each draw of the rvars it contains, then returns a new rvar representing the output of those expressions. To identify rvars, rdo() searches the calling environment for any variables named in expr for which is_rvar() evaluates to TRUE. If expr contains no rvars, then it will be executed ndraws times and an rvar with that many draws returned.

rdo() is not necessarily fast (in fact in some cases it may be very slow), but it has the advantage of allowing a nearly arbitrary R expression to be executed against rvars simply by wrapping it with rdo( ... ). This makes it especially useful as a prototyping tool. If you create code with rdo() and it is unacceptably slow for your application, consider rewriting it using math operations directly on rvars (which should be fast), using rvar_rng(), and/or using operations directly on the arrays that back the rvars (via draws_of()).

Value

An rvar.

See Also

Other rfun: rfun(), rvar_rng()

Examples

mu <- rdo(rnorm(10, mean = 1:10, sd = 1))
sigma <- rdo(rgamma(1, shape = 1, rate = 1))
x <- rdo(rnorm(10, mu, sigma))
x

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

distributional

cdf, variance

Rename variables in draws objects

Description

Rename variables in a draws object.

Usage

rename_variables(.x, ...)

## S3 method for class 'draws'
rename_variables(.x, ...)

Arguments

Value

Returns a draws object of the same format as .x, with variables renamed according to the expressions provided in ....

See Also

variables, variables<-, mutate_variables

Examples

x <- as_draws_df(example_draws())
variables(x)

x <- rename_variables(x, mean = mu, sigma = tau)
variables(x)

x <- rename_variables(x, b = `theta[1]`) # or b  = "theta[1]"
variables(x)

# rename all elements of 'theta' at once
x <- rename_variables(x, alpha = theta)
variables(x)

Repair indices of draws objects

Description

Repair indices of draws objects so that iterations, chains, and draws are continuously and consistently numbered.

Usage

repair_draws(x, order = TRUE, ...)

## S3 method for class 'draws_matrix'
repair_draws(x, order = TRUE, ...)

## S3 method for class 'draws_array'
repair_draws(x, order = TRUE, ...)

## S3 method for class 'draws_df'
repair_draws(x, order = TRUE, ...)

## S3 method for class 'draws_list'
repair_draws(x, order = TRUE, ...)

## S3 method for class 'draws_rvars'
repair_draws(x, order = TRUE, ...)

## S3 method for class 'rvar'
repair_draws(x, order = TRUE, ...)

Arguments

Value

A draws object of the same class as x.

See Also

order_draws()

Examples

x <- as_draws_array(example_draws())
(x <- x[10:5, 3:4, ])
repair_draws(x)

Resample draws objects

Description

Resample draws objects according to provided weights, for example weights obtained through importance sampling.

Usage

resample_draws(x, ...)

## S3 method for class 'draws'
resample_draws(x, weights = NULL, method = "stratified", ndraws = NULL, ...)

## S3 method for class 'rvar'
resample_draws(x, ...)

Arguments

Details

Upon usage of resample_draws(), chains will automatically be merged due to subsetting of individual draws (see subset_draws for details). Also, weights stored in the draws object will be removed in the process, as resampling invalidates existing weights.

Value

A draws object of the same class as x.

References

Kitagawa, G., Monte Carlo Filter and Smoother for Non-Gaussian Nonlinear ' State Space Models, Journal of Computational and Graphical Statistics, 5(1):1-25, 1996.

See Also

resample_draws()

Examples

x <- as_draws_df(example_draws())

# random weights for justr for demonstration
w <- runif(ndraws(x), 0, 10)

# use default stratified sampling
x_rs <- resample_draws(x, weights = w)
summarise_draws(x_rs, default_summary_measures())

# use simple random sampling
x_rs <- resample_draws(x, weights = w, method = "simple")
summarise_draws(x_rs, default_summary_measures())

Reserved variables

Description

Get names of reserved variables from objects in the posterior package.

Usage

reserved_variables(x, ...)

## Default S3 method:
reserved_variables(x, ...)

## S3 method for class 'draws_matrix'
reserved_variables(x, ...)

## S3 method for class 'draws_array'
reserved_variables(x, ...)

## S3 method for class 'draws_df'
reserved_variables(x, ...)

## S3 method for class 'draws_list'
reserved_variables(x, ...)

## S3 method for class 'draws_rvars'
reserved_variables(x, ...)

Arguments

Details

reserved_variables() returns the names of reserved variables in use by an object.

The following variables names are currently reserved for special use cases in all draws formats:

• .log_weight: Log weights per draw (see weight_draws).

Further, specific for the draws_df format, there are three additional reserved variables:

• .chain: Chain index per draw

• .iteration: Iteration index within each chain

• .draw: Draw index across chains

More reserved variables may be added in the future.

Value

A character vector of reserved variables used in x.

Examples

x <- example_draws()
reserved_variables(x)

# if we add weights, the `.log_weight` reserved variable is used
x <- weight_draws(x, rexp(ndraws(x)))
reserved_variables(x)

Create functions of random variables

Description

Function that create functions that can accept and/or produce rvars.

Usage

rfun(.f, rvar_args = NULL, rvar_dots = TRUE, ndraws = NULL)

Arguments

Details

This function wraps an existing function (.f) such that it returns rvars containing whatever type of data .f would normally return.

The returned function, when called, executes .f possibly multiple times, once for each draw of the rvars passed to it, then returns a new rvar representing the output of those function evaluations. If the arguments contain no rvars, then .f will be executed ndraws times and an rvar with that many draws returned.

Functions created by rfun() are not necessarily fast (in fact in some cases they may be very slow), but they have the advantage of allowing a nearly arbitrary R functions to be executed against rvars simply by wrapping them with rfun(). This makes it especially useful as a prototyping tool. If you create code with rfun() and it is unacceptably slow for your application, consider rewriting it using math operations directly on rvars (which should be fast), using rvar_rng(), and/or using operations directly on the arrays that back the rvars (via draws_of()).

Value

A function with the same argument specification as .f, but which can accept and return rvars.

See Also

Other rfun: rdo(), rvar_rng()

Examples

rvar_norm <- rfun(rnorm)
rvar_gamma <- rfun(rgamma)

mu <- rvar_norm(10, mean = 1:10, sd = 1)
sigma <- rvar_gamma(1, shape = 1, rate = 1)
x <- rvar_norm(10, mu, sigma)
x

Rhat convergence diagnostic

Description

Compute the Rhat convergence diagnostic for a single variable as the maximum of rank normalized split-Rhat and rank normalized folded-split-Rhat as proposed in Vehtari et al. (2021).

Usage

rhat(x, ...)

## Default S3 method:
rhat(x, ...)

## S3 method for class 'rvar'
rhat(x, ...)

Arguments

Value

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, 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 Analysis. 16(2), 667-–718. doi:10.1214/20-BA1221

See Also

Other diagnostics: ess_basic(), ess_bulk(), ess_quantile(), ess_sd(), ess_tail(), mcse_mean(), mcse_quantile(), mcse_sd(), pareto_diags(), pareto_khat(), rhat_basic(), rhat_nested(), rstar()