Title:  Tools for Working with Posterior Distributions 

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: (a) Efficiently convert between many different useful formats of draws (samples) from posterior or prior distributions. (b) Provide consistent methods for operations commonly performed on draws, for example, subsetting, binding, or mutating draws. (c) Provide various summaries of draws in convenient formats. (d) Provide lightweight implementations of state of the art posterior inference diagnostics. References: Vehtari et al. (2021) <doi:10.1214/20BA1221>. 
Authors:  PaulChristian Bürkner [aut, cre], Jonah Gabry [aut], Matthew Kay [aut], Aki Vehtari [aut], Måns Magnusson [ctb], Rok Češnovar [ctb], Ben Lambert [ctb], Ozan Adıgüzel [ctb], Jacob Socolar [ctb], Noa Kallioinen [ctb] 
Maintainer:  PaulChristian Bürkner <[email protected]> 
License:  BSD_3_clause + file LICENSE 
Version:  1.6.0 
Built:  20240806 08:19:09 UTC 
Source:  https://github.com/standev/posterior 
https://mcstan.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.
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
rvar
s 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 rvar
s 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.
Maintainer: PaulChristian Bürkner [email protected]
Authors:
Jonah Gabry [email protected]
Matthew Kay [email protected]
Aki Vehtari [email protected]
Other contributors:
Måns Magnusson [contributor]
Rok Češnovar [contributor]
Ben Lambert [contributor]
Ozan Adıgüzel [contributor]
Jacob Socolar [contributor]
Noa Kallioinen [contributor]
Useful links:
Report bugs at https://github.com/standev/posterior/issues
Convert x
to an rvar
object.
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)
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)
x 
(multiple options) An object that can be converted to an 
dim 
(integer vector) One or more integers giving the maximal indices
in each dimension to override the dimensions of the 
dimnames 
(list) Character vectors giving the names in each dimension
to override the names of the dimensions of the 
nchains 
(positive integer) The number of chains. The default is 
For objects that are already rvar
s, 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.
An object of class "rvar"
(or one of its subtypes) representing a random variable.
rvar()
to construct rvar
s directly. See rdo()
, rfun()
, and
rvar_rng()
for higherlevel interfaces for creating rvar
s.
# 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)
# 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)
Convert x
to an rvar_factor
or rvar_ordered
object.
as_rvar_factor(x, dim = NULL, dimnames = NULL, nchains = NULL, ...) as_rvar_ordered(x, dim = NULL, dimnames = NULL, nchains = NULL, ...)
as_rvar_factor(x, dim = NULL, dimnames = NULL, nchains = NULL, ...) as_rvar_ordered(x, dim = NULL, dimnames = NULL, nchains = NULL, ...)
x 
(multiple options) An object that can be converted to an 
dim 
(integer vector) One or more integers giving the maximal indices
in each dimension to override the dimensions of the 
dimnames 
(list) Character vectors giving the names in each dimension
to override the names of the dimensions of the 
nchains 
(positive integer) The number of chains. The default is 
... 
Arguments passed on to

For objects that are already rvar
s, 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 factorlike).
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.
An object of class "rvar_factor"
or "rvar_ordered"
representing a random variable.
rvar()
, rvar_factor()
, and rvar_ordered()
to construct rvar
s directly.
See rdo()
, rfun()
, and rvar_rng()
for higherlevel interfaces for creating rvar
s.
# 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)))
# 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)))
draws
objects togetherBind multiple draws
objects together to form a single draws
object.
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")
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")
x 
(draws) A 
... 
(draws) Additional 
along 
(string) The dimension along which draws objects should be bound
together. Possible values are 
A draws
object of the same class as x
.
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)
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 an rvar
containing a matrix.
## S3 method for class 'rvar' chol(x, ...)
## S3 method for class 'rvar' chol(x, ...)
x 
(rvar) A 2dimensional 
... 
Additional parameters passed on to 
An rvar
containing the upper triangular factor of the Cholesky
decomposition, i.e., the matrix $R$
such that $R'R = x$
.
Extract the diagonal of a matrix or construct a matrix, including random
matrices (2dimensional rvar
s). Makes base::diag()
generic.
## S4 method for signature 'rvar' diag(x = 1, nrow, ncol, names = TRUE)
## S4 method for signature 'rvar' diag(x = 1, nrow, ncol, names = TRUE)
x 
(numeric,rvar) a matrix, vector, 1D array, missing, or a 1 or
2dimensional 
nrow , ncol

optional dimensions for the result when 
names 
(when 
Makes base::diag()
into a generic function. See that function's documentation
for usage with numeric
s and for usage of diag<
, which is also supported
by rvar
.
For rvar
s, has two modes:
x
is a matrixlike rvar
: it returns the diagonal as a vectorlike rvar
x
is a vectorlike rvar
: it returns a matrixlike rvar
with x
as
the diagonal and zero for offdiagonal entries.
# 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))
# 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))
A list of available diagnostics and links to their individual help pages.
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, rankbased version of Rhat 
rhat_nested() 
Rhat for use with many short chains 
rstar() 
R* diagnostic 
See individual functions for a description of return types.
Dissention, for measuring dispersion in draws from ordinal distributions.
dissent(x) ## Default S3 method: dissent(x) ## S3 method for class 'rvar' dissent(x)
dissent(x) ## Default S3 method: dissent(x) ## S3 method for class 'rvar' dissent(x)
x 
(multiple options) A vector to be interpreted as draws from an ordinal distribution, such as:

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).
If x
is a factor or numeric, returns a length1 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
.
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.
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)
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)
draws
objectsTry to transform an R object to a format supported by the posterior package.
as_draws(x, ...) is_draws(x)
as_draws(x, ...) is_draws(x)
x 
(draws) A 
... 
Arguments passed to individual methods (if applicable). 
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"
).
If possible, a draws
object in the closest supported format to x
.
The formats are linked to in the See Also section below.
Other formats:
draws_array()
,
draws_df()
,
draws_list()
,
draws_matrix()
,
draws_rvars()
# 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)
# 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)
draws_array
formatThe 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.
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)
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)
x 
An object to convert to a 
... 
For 
.nchains 
(positive integer) The number of chains. The default is 
Objects of class "draws_array"
are 3D arrays with dimensions
"iteration"
, "chain"
, and "variable"
. See Examples.
A draws_array
object, which has classes
c("draws_array", "draws", "array")
.
Other formats:
draws
,
draws_df()
,
draws_list()
,
draws_matrix()
,
draws_rvars()
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)
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)
draws_df
formatThe 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.
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)
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)
x 
An object to convert to a 
... 
For 
.nchains 
(positive integer) The number of chains. The default is 
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.
A draws_df
object, which has classes
c("draws_df", "draws", class(tibble::tibble()))
.
Other formats:
draws
,
draws_array()
,
draws_list()
,
draws_matrix()
,
draws_rvars()
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)
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)
draws_list
formatThe 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.
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)
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)
x 
An object to convert to a 
... 
For 
.nchains 
(positive integer) The number of chains. The default is 
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.
A draws_list
object, which has classes
c("draws_list", "draws", "list")
.
Other formats:
draws
,
draws_array()
,
draws_df()
,
draws_matrix()
,
draws_rvars()
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)
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)
draws_matrix
formatThe 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.
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)
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)
x 
An object to convert to a 
... 
For 
.nchains 
(positive integer) The number of chains. The default is 
Objects of class "draws_matrix"
are matrices (2D arrays) with
dimensions "draw"
and "variable"
. See Examples.
A draws_matrix
object, which has classes
c("draws_matrix", "draws", "matrix")
.
Other formats:
draws
,
draws_array()
,
draws_df()
,
draws_list()
,
draws_rvars()
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)
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)
Gets/sets the arrayrepresentation that backs an rvar
. Should be used rarely.
draws_of(x, with_chains = FALSE) draws_of(x, with_chains = FALSE) < value
draws_of(x, with_chains = FALSE) draws_of(x, with_chains = FALSE) < value
x 
(rvar) An 
with_chains 
(logical) Should the array of draws include a dimension for chains?
If 
value 
(array) An array of values to use as the backing array of 
While rvar
s 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.
rvar
s 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.
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))
.
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
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
draws_rvars
formatThe 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.
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)
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)
x 
An object to convert to a 
... 
For 
.nchains 
(positive integer) The number of chains. The default is 
Objects of class "draws_rvars"
are lists of rvar
objects.
See Examples.
A draws_rvars
object, which has classes
c("draws_rvars", "draws", "list")
.
Other formats:
draws
,
draws_array()
,
draws_df()
,
draws_list()
,
draws_matrix()
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)
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)
draws
objectsThe 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.
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()
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()
.x , object

(draws) A 
... 
Namevalue pairs of summary or diagnostic functions. The provided names will be used as the names of the columns in the result unless the function returns a named vector, in which case the latter names are used. The functions can be specified in any format supported by as_function(). See Examples. 
.args 
(named list) Optional arguments passed to the summary functions. 
.num_args 
(named list) Optional arguments passed to
num() for pretty printing of summaries. Can be controlled
globally via the 
.cores 
(positive integer) The number of cores to use for computing
summaries for different variables in parallel. Coerced to integer if
possible, otherwise errors. The default is 
The default summary functions used are the ones specified by
default_summary_measures()
and default_convergence_measures()
:
default_summary_measures()
default_convergence_measures()
The var()
function should not be used to compute variances due
to its inconsistent behavior with matrices. Instead, please use
distributional::variance()
.
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.
diagnostics
for a list of available diagnostics and links to
their individual help pages.
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"))
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"))
draws
objectsIndex iterations, chains, and draws of draws
objects.
iteration_ids(x) chain_ids(x) draw_ids(x) niterations(x) nchains(x) ndraws(x)
iteration_ids(x) chain_ids(x) draw_ids(x) niterations(x) nchains(x) ndraws(x)
x 
(draws) A 
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.
For iteration_ids()
, chain_ids()
, and draw_ids()
, an integer vector.
For niterations()
, nchains()
, and ndraws()
, a scalar integer.
x < example_draws() iteration_ids(x) niterations(x) chain_ids(x) nchains(x) draw_ids(x) ndraws(x)
x < example_draws() iteration_ids(x) niterations(x) chain_ids(x) nchains(x) draw_ids(x) ndraws(x)
Delete the dimensions of an rvar
which are of size one. See base::drop()
## S4 method for signature 'rvar' drop(x)
## S4 method for signature 'rvar' drop(x)
x 
(rvar) an 
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 rvar
s, unlike numeric
s,
never have NULL
dimensions).
# 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]
# 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, for measuring dispersion in draws from categorical distributions.
entropy(x) ## Default S3 method: entropy(x) ## S3 method for class 'rvar' entropy(x)
entropy(x) ## Default S3 method: entropy(x) ## S3 method for class 'rvar' entropy(x)
x 
(multiple options) A vector to be interpreted as draws from a categorical distribution, such as:

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).
If x
is a factor or numeric, returns a length1 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
.
Allen R. Wilcox (1967). Indices of Qualitative Variation (No. ORNLTM1919). Oak Ridge National Lab., Tenn.
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)
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)
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 indepth
comparison of different effective sample size estimators.
ess_basic(x, ...) ## Default S3 method: ess_basic(x, split = TRUE, ...) ## S3 method for class 'rvar' ess_basic(x, split = TRUE, ...)
ess_basic(x, ...) ## Default S3 method: ess_basic(x, split = TRUE, ...) ## S3 method for class 'rvar' ess_basic(x, split = TRUE, ...)
x 
(multiple options) One of:

... 
Arguments passed to individual methods (if applicable). 
split 
(logical) Should the estimate be computed on split chains? The
default is 
If the input is an array, returns a single numeric value. If any of the draws
is nonfinite, 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.
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 PaulChristian Bürkner (2021). Ranknormalization, folding, and localization: An improved Rhat for assessing convergence of MCMC (with discussion). Bayesian Analysis. 16(2), 667–718. doi:10.1214/20BA1221
Aki Vehtari (2021). Comparison of MCMC effective sample size estimators. Retrieved from https://avehtari.github.io/rhat_ess/ess_comparison.html
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()
mu < extract_variable_matrix(example_draws(), "mu") ess_basic(mu) d < as_draws_rvars(example_draws("multi_normal")) ess_basic(d$Sigma)
mu < extract_variable_matrix(example_draws(), "mu") ess_basic(mu) d < as_draws_rvars(example_draws("multi_normal")) ess_basic(d$Sigma)
Compute a bulk effective sample size estimate (bulkESS) for a single
variable. BulkESS 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 indepth
comparison of different effective sample size estimators.
ess_bulk(x, ...) ## Default S3 method: ess_bulk(x, ...) ## S3 method for class 'rvar' ess_bulk(x, ...)
ess_bulk(x, ...) ## Default S3 method: ess_bulk(x, ...) ## S3 method for class 'rvar' ess_bulk(x, ...)
x 
(multiple options) One of:

... 
Arguments passed to individual methods (if applicable). 
If the input is an array, returns a single numeric value. If any of the draws
is nonfinite, 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.
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and PaulChristian Bürkner (2021). Ranknormalization, folding, and localization: An improved Rhat for assessing convergence of MCMC (with discussion). Bayesian Analysis. 16(2), 667–718. doi:10.1214/20BA1221
Aki Vehtari (2021). Comparison of MCMC effective sample size estimators. Retrieved from https://avehtari.github.io/rhat_ess/ess_comparison.html
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()
mu < extract_variable_matrix(example_draws(), "mu") ess_bulk(mu) d < as_draws_rvars(example_draws("multi_normal")) ess_bulk(d$Sigma)
mu < extract_variable_matrix(example_draws(), "mu") ess_bulk(mu) d < as_draws_rvars(example_draws("multi_normal")) ess_bulk(d$Sigma)
Compute an effective sample size estimate for a mean (expectation) estimate of a single variable.
ess_mean(x, ...) ## S3 method for class 'rvar' ess_mean(x, ...)
ess_mean(x, ...) ## S3 method for class 'rvar' ess_mean(x, ...)
x 
(multiple options) One of:

... 
Arguments passed to individual methods (if applicable). 
If the input is an array, returns a single numeric value. If any of the draws
is nonfinite, 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.
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.
mu < extract_variable_matrix(example_draws(), "mu") ess_mean(mu) d < as_draws_rvars(example_draws("multi_normal")) ess_mean(d$Sigma)
mu < extract_variable_matrix(example_draws(), "mu") ess_mean(mu) d < as_draws_rvars(example_draws("multi_normal")) ess_mean(d$Sigma)
Compute effective sample size estimates for quantile estimates of a single variable.
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, ...) ## Default S3 method: ess_mean(x, ...)
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, ...) ## Default S3 method: ess_mean(x, ...)
x 
(multiple options) One of:

probs 
(numeric vector) Probabilities in 
... 
Arguments passed to individual methods (if applicable). 
names 
(logical) Should the result have a 
If the input is an array,
returns a numeric vector with one element per quantile. If any of the draws is
nonfinite, 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))
.
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and PaulChristian Bürkner (2021). Ranknormalization, folding, and localization: An improved Rhat for assessing convergence of MCMC (with discussion). Bayesian Analysis. 16(2), 667–718. doi:10.1214/20BA1221
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()
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))
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))
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.
ess_sd(x, ...) ## Default S3 method: ess_sd(x, ...) ## S3 method for class 'rvar' ess_sd(x, ...)
ess_sd(x, ...) ## Default S3 method: ess_sd(x, ...) ## S3 method for class 'rvar' ess_sd(x, ...)
x 
(multiple options) One of:

... 
Arguments passed to individual methods (if applicable). 
If the input is an array, returns a single numeric value. If any of the draws
is nonfinite, 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.
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and PaulChristian Bürkner (2021). Ranknormalization, folding, and localization: An improved Rhat for assessing convergence of MCMC (with discussion). Bayesian Analysis. 16(2), 667–718. doi:10.1214/20BA1221
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()
mu < extract_variable_matrix(example_draws(), "mu") ess_sd(mu) d < as_draws_rvars(example_draws("multi_normal")) ess_sd(d$Sigma)
mu < extract_variable_matrix(example_draws(), "mu") ess_sd(mu) d < as_draws_rvars(example_draws("multi_normal")) ess_sd(d$Sigma)
Compute a tail effective sample size estimate (tailESS) for a single
variable. TailESS 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 indepth
comparison of different effective sample size estimators.
ess_tail(x, ...) ## Default S3 method: ess_tail(x, ...) ## S3 method for class 'rvar' ess_tail(x, ...)
ess_tail(x, ...) ## Default S3 method: ess_tail(x, ...) ## S3 method for class 'rvar' ess_tail(x, ...)
x 
(multiple options) One of:

... 
Arguments passed to individual methods (if applicable). 
If the input is an array, returns a single numeric value. If any of the draws
is nonfinite, 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.
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and PaulChristian Bürkner (2021). Ranknormalization, folding, and localization: An improved Rhat for assessing convergence of MCMC (with discussion). Bayesian Analysis. 16(2), 667–718. doi:10.1214/20BA1221
Aki Vehtari (2021). Comparison of MCMC effective sample size estimators. Retrieved from https://avehtari.github.io/rhat_ess/ess_comparison.html
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()
mu < extract_variable_matrix(example_draws(), "mu") ess_tail(mu) d < as_draws_rvars(example_draws("multi_normal")) ess_tail(d$Sigma)
mu < extract_variable_matrix(example_draws(), "mu") ess_tail(mu) d < as_draws_rvars(example_draws("multi_normal")) ess_tail(d$Sigma)
draws
objectsObjects for use in examples, vignettes, and tests.
example_draws(example = "eight_schools")
example_draws(example = "eight_schools")
example 
(string) The name of the example 
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 3dimensional 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
A draws
object.
These objects are only intended to be used in demonstrations and tests. They contain fewer iterations and chains than recommended for performing actual inference.
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.
draws_eight_schools < example_draws("eight_schools") summarise_draws(draws_eight_schools) draws_multi_normal < example_draws("multi_normal") summarise_draws(draws_multi_normal)
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 a vector of draws of a single variable.
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, ...)
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, ...)
x 
(draws) A 
variable 
(string) The name of the variable to extract. Must include
indices for array variables (e.g. 
... 
Arguments passed to individual methods (if applicable). 
A vector of length equal to the number of draws.
Other variable extraction methods:
extract_variable_array()
,
extract_variable_matrix()
x < example_draws() mu < extract_variable(x, variable = "mu") str(mu)
x < example_draws() mu < extract_variable(x, variable = "mu") str(mu)
Extract an array of draws of a single variable, including any dimensions of variables with indices.
extract_variable_array(x, variable, ...) ## Default S3 method: extract_variable_array(x, variable, ...) ## S3 method for class 'draws' extract_variable_array(x, variable, ...)
extract_variable_array(x, variable, ...) ## Default S3 method: extract_variable_array(x, variable, ...) ## S3 method for class 'draws' extract_variable_array(x, variable, ...)
x 
(draws) A 
variable 
(string) The name of the variable to extract. To extract all
dimensions from variables with indices (e.g. 
... 
Arguments passed to individual methods (if applicable). 
An array
with dimension niterations(x)
x nchains(x)
x any remaining
dimensions determined by the indices of the variable x
.
Other variable extraction methods:
extract_variable()
,
extract_variable_matrix()
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)
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 an iterations x chains matrix of draws of a single variable.
This is primarily used for convergence diagnostic functions such as rhat()
.
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, ...)
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, ...)
x 
(draws) A 
variable 
(string) The name of the variable to extract. Must include
indices for array variables (e.g. 
... 
Arguments passed to individual methods (if applicable). 
A matrix
with dimension iterations x chains.
Other variable extraction methods:
extract_variable()
,
extract_variable_array()
x < example_draws() mu < extract_variable_matrix(x, variable = "mu") dim(mu) rhat(mu)
x < example_draws() mu < extract_variable_matrix(x, variable = "mu") dim(mu) rhat(mu)
Executes an expression once for every draw in a draws
object. Used
primarily for its side effects and returns the input x
invisibly.
for_each_draw(x, expr)
for_each_draw(x, expr)
x 
(draws) A 
expr 
(expression) A bare expression that can contain references to
variables in 
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 <<
.)
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.
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) )
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) )
x
a random variable?Test if x
is an rvar
.
is_rvar(x)
is_rvar(x)
x 
(any object) An object to test. 
TRUE
if x
is an rvar
, FALSE
otherwise.
as_rvar()
to convert objects to rvar
s.
x
a factor random variable?Test if x
is an rvar_factor
or rvar_ordered
.
is_rvar_factor(x) is_rvar_ordered(x)
is_rvar_factor(x) is_rvar_ordered(x)
x 
(any object) An object to test. 
TRUE
if x
is an rvar_factor
or rvar_ordered
, FALSE
otherwise.
as_rvar_factor()
and as_rvar_ordered()
to convert objects to
rvar_factor
s and rvar_ordered
s.
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.
match(x, table, ...) ## Default S3 method: match(x, ...) ## S3 method for class 'rvar' match(x, ...) x %in% table
match(x, table, ...) ## Default S3 method: match(x, ...) ## S3 method for class 'rvar' match(x, ...) x %in% table
x 
(multiple options) the values to be matched. Can be:

table 
(vector) the values to be matched against. 
... 
Arguments passed on to

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 nonbase vectors (such
as rvars) are supported.
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
.
x < rvar(c("a","b","b","c","d")) x %in% c("b","d") # for additional examples, see base::match()
x < rvar(c("a","b","b","c","d")) x %in% c("b","d") # for additional examples, see base::match()
Compute the Monte Carlo standard error for the mean (expectation) of a single variable.
mcse_mean(x, ...) ## Default S3 method: mcse_mean(x, ...) ## S3 method for class 'rvar' mcse_mean(x, ...)
mcse_mean(x, ...) ## Default S3 method: mcse_mean(x, ...) ## S3 method for class 'rvar' mcse_mean(x, ...)
x 
(multiple options) One of:

... 
Arguments passed to individual methods (if applicable). 
If the input is an array, returns a single numeric value. If any of the draws
is nonfinite, 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.
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.
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()
mu < extract_variable_matrix(example_draws(), "mu") mcse_mean(mu) d < as_draws_rvars(example_draws("multi_normal")) mcse_mean(d$Sigma)
mu < extract_variable_matrix(example_draws(), "mu") mcse_mean(mu) d < as_draws_rvars(example_draws("multi_normal")) mcse_mean(d$Sigma)
Compute Monte Carlo standard errors for quantile estimates of a single variable.
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, ...)
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, ...)
x 
(multiple options) One of:

probs 
(numeric vector) Probabilities in 
... 
Arguments passed to individual methods (if applicable). 
names 
(logical) Should the result have a 
If the input is an array,
returns a numeric vector with one element per quantile. If any of the draws is
nonfinite, 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))
.
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and PaulChristian Bürkner (2021). Ranknormalization, folding, and localization: An improved Rhat for assessing convergence of MCMC (with discussion). Bayesian Analysis. 16(2), 667–718. doi:10.1214/20BA1221
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()
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)
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)
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).
mcse_sd(x, ...) ## Default S3 method: mcse_sd(x, ...) ## S3 method for class 'rvar' mcse_sd(x, ...)
mcse_sd(x, ...) ## Default S3 method: mcse_sd(x, ...) ## S3 method for class 'rvar' mcse_sd(x, ...)
x 
(multiple options) One of:

... 
Arguments passed to individual methods (if applicable). 
If the input is an array, returns a single numeric value. If any of the draws
is nonfinite, 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.
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and PaulChristian Bürkner (2021). Ranknormalization, folding, and localization: An improved Rhat for assessing convergence of MCMC (with discussion). Bayesian Analysis. 16(2), 667–718. doi:10.1214/20BA1221
J. F. Kenney & E. S. Keeping (1951). Mathematics of Statistics, Vol. II.
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()
mu < extract_variable_matrix(example_draws(), "mu") mcse_sd(mu) d < as_draws_rvars(example_draws("multi_normal")) mcse_sd(d$Sigma)
mu < extract_variable_matrix(example_draws(), "mu") mcse_sd(mu) d < as_draws_rvars(example_draws("multi_normal")) mcse_sd(d$Sigma)
draws
objectsMerge 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)
.
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, ...)
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, ...)
x 
(draws) A 
... 
Arguments passed to individual methods (if applicable). 
A draws
object of the same class as x
.
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))
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 of a vector.
modal_category(x) ## Default S3 method: modal_category(x) ## S3 method for class 'rvar' modal_category(x)
modal_category(x) ## Default S3 method: modal_category(x) ## S3 method for class 'rvar' modal_category(x)
x 
(multiple options) A vector to be interpreted as draws from a categorical distribution, such as:

Finds the modal category (i.e., most frequent value) in x
. In the case of
ties, returns the first tie.
If x
is a factor or numeric, returns a length1 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
.
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)
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)
draws
objectsMutate variables in a draws
object.
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, ...)
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, ...)
.x 
(draws) A 
... 
Namevalue pairs of expressions, each with either length 1 or the
same length as in the entire input (i.e., number of iterations or draws).
The name of each argument will be the name of a new variable, and the value
will be its corresponding value. Use a 
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()
.
Returns a draws
object of the same format as .x
, with variables mutated
according to the expressions provided in ...
.
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)
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)
draws
objectsOrder 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.
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, ...)
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, ...)
x 
(draws) A 
... 
Arguments passed to individual methods (if applicable). 
A draws
object of the same class as x
.
x < as_draws_array(example_draws()) dimnames(x[10:5, 4:3, ]) dimnames(order_draws(x[10:5, 4:3, ]))
x < as_draws_array(example_draws()) dimnames(x[10:5, 4:3, ]) dimnames(order_draws(x[10:5, 4:3, ]))
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).
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, ...)
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, ...)
x 
(multiple options) One of:

... 
Arguments passed to individual methods (if applicable). 
tail 
(string) The tail to diagnose/smooth:
The default is 
r_eff 
(numeric) relative effective sample size estimate. If

ndraws_tail 
(numeric) number of draws for the tail. If

verbose 
(logical) Should diagnostic messages be printed? If

are_log_weights 
(logical) Are the draws log weights? Default is

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 khat values result in
reliable Pareto smoothed estimates. The threshold is lower for
smaller effective sample sizes. If khat 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).
List of Pareto smoothing diagnostics:
khat
: estimated Pareto k shape parameter,
min_ss
: minimum sample size for reliable Pareto smoothed estimate,
khat_threshold
: khatthreshold for reliable Pareto smoothed estimate,
convergence_rate
: Pareto smoothed estimate RMSE convergence rate.
Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao and Jonah Gabry (2024). Pareto Smoothed Importance Sampling. Journal of Machine Learning Research, 25(72):158. PDF
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()
mu < extract_variable_matrix(example_draws(), "mu") pareto_diags(mu) d < as_draws_rvars(example_draws("multi_normal")) pareto_diags(d$Sigma)
mu < extract_variable_matrix(example_draws(), "mu") pareto_diags(mu) d < as_draws_rvars(example_draws("multi_normal")) pareto_diags(d$Sigma)
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).
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, ...)
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, ...)
x 
(multiple options) One of:

... 
Arguments passed to individual methods (if applicable). 
tail 
(string) The tail to diagnose/smooth:
The default is 
r_eff 
(numeric) relative effective sample size estimate. If

ndraws_tail 
(numeric) number of draws for the tail. If

verbose 
(logical) Should diagnostic messages be printed? If

are_log_weights 
(logical) Are the draws log weights? Default is

If the input is an array, returns a single numeric value. If any of the draws
is nonfinite, 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.
Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao and Jonah Gabry (2024). Pareto Smoothed Importance Sampling. Journal of Machine Learning Research, 25(72):158. PDF
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()
mu < extract_variable_matrix(example_draws(), "mu") pareto_khat(mu) d < as_draws_rvars(example_draws("multi_normal")) pareto_khat(d$Sigma)
mu < extract_variable_matrix(example_draws(), "mu") pareto_khat(mu) d < as_draws_rvars(example_draws("multi_normal")) pareto_khat(d$Sigma)
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).
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, ... )
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, ... )
x 
(multiple options) One of:

... 
Arguments passed to individual methods (if applicable). 
return_k 
(logical) Should the Pareto khat be included in
output? If 
extra_diags 
(logical) Should extra Pareto khat diagnostics
be included in output? If 
tail 
(string) The tail to diagnose/smooth:
The default is 
r_eff 
(numeric) relative effective sample size estimate. If

ndraws_tail 
(numeric) number of draws for the tail. If

verbose 
(logical) Should diagnostic messages be printed? If

are_log_weights 
(logical) Are the draws log weights? Default is

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 nonfinite, 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).
Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao and Jonah Gabry (2024). Pareto Smoothed Importance Sampling. Journal of Machine Learning Research, 25(72):158. PDF
pareto_khat
for only calculating khat, and
pareto_diags
for additional diagnostics.
mu < extract_variable_matrix(example_draws(), "mu") pareto_smooth(mu) d < as_draws_rvars(example_draws("multi_normal")) pareto_smooth(d$Sigma)
mu < extract_variable_matrix(example_draws(), "mu") pareto_smooth(mu) d < as_draws_rvars(example_draws("multi_normal")) pareto_smooth(d$Sigma)
draws_array
objectsPretty printing for draws_array
objects.
## 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, ... )
## 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, ... )
x 
(draws) A 
digits 
(nonnegative integer) The minimum number of significant digits
to print. If 
max_iterations 
(positive integer) The maximum number of iterations to
print. Can be controlled globally via the 
max_chains 
(positive integer) The maximum number of chains to print.
Can be controlled globally via the 
max_variables 
(positive integer) The maximum number of variables to
print. Can be controlled globally via the 
reserved 
(logical) Should reserved variables be included in the
output? Defaults to 
... 
Further arguments passed to the underlying 
A draws
object of the same class as x
.
x < as_draws_array(example_draws()) print(x)
x < as_draws_array(example_draws()) print(x)
draws_df
objectsPretty printing for draws_df
objects.
## 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, ... )
## 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, ... )
x 
(draws) A 
digits 
(nonnegative integer) The minimum number of significant digits
to print. If 
max_draws 
(positive integer) The maximum number of draws to print. Can
be controlled globally via the 
max_variables 
(positive integer) The maximum number of variables to
print. Can be controlled globally via the 
reserved 
(logical) Should reserved variables be included in the
output? Defaults to 
... 
Further arguments passed to the underlying 
A draws
object of the same class as x
.
x < as_draws_df(example_draws()) print(x)
x < as_draws_df(example_draws()) print(x)
draws_list
objectsPretty printing for draws_list
objects.
## 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, ... )
## 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, ... )
x 
(draws) A 
digits 
(nonnegative integer) The minimum number of significant digits
to print. If 
max_iterations 
(positive integer) The maximum number of iterations to
print. Can be controlled globally via the 
max_chains 
(positive integer) The maximum number of chains to print.
Can be controlled globally via the 
max_variables 
(positive integer) The maximum number of variables to
print. Can be controlled globally via the 
reserved 
(logical) Should reserved variables be included in the
output? Defaults to 
... 
Further arguments passed to the underlying 
A draws
object of the same class as x
.
x < as_draws_list(example_draws()) print(x)
x < as_draws_list(example_draws()) print(x)
draws_matrix
objectsPretty printing for draws_matrix
objects.
## 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, ... )
## 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, ... )
x 
(draws) A 
digits 
(nonnegative integer) The minimum number of significant digits
to print. If 
max_draws 
(positive integer) The maximum number of draws to print. Can
be controlled globally via the 
max_variables 
(positive integer) The maximum number of variables to
print. Can be controlled globally via the 
reserved 
(logical) Should reserved variables be included in the
output? Defaults to 
... 
Further arguments passed to the underlying 
A draws
object of the same class as x
.
x < as_draws_matrix(example_draws()) print(x)
x < as_draws_matrix(example_draws()) print(x)
draws_rvars
objectsPretty printing for draws_rvars
objects.
## 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, ... )
## 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, ... )
x 
(draws) A 
digits 
(nonnegative integer) The minimum number of significant digits
to print. If 
max_variables 
(positive integer) The maximum number of variables to
print. Can be controlled globally via the 
summary 
(string) The style of summary to display:

reserved 
(logical) Should reserved variables be included in the
output? Defaults to 
... 
Further arguments passed to the underlying 
A draws
object of the same class as x
.
x < as_draws_rvars(example_draws()) print(x)
x < as_draws_rvars(example_draws()) print(x)
draws
objectsPrint output from summarise_draws()
.
## S3 method for class 'draws_summary' print(x, ..., num_args = NULL)
## S3 method for class 'draws_summary' print(x, ..., num_args = NULL)
x 
(draws_summary) A 
... 
Additional arguments passed to 
num_args 
(named list) Optional arguments passed to
num() for pretty printing of summaries. If 
An invisible version of the input object.
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))
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))
Printing and formatting methods for rvar
s.
## 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 )
## 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 )
x , object

(rvar) The 
... 
Further arguments passed to the underlying 
summary 
(string) The style of summary to display:

digits 
(nonnegative integer) The minimum number of significant digits
to print. If 
color 
(logical) Whether or not to use color when formatting the
output. If 
width 
The maxmimum width used to print out lists of factor levels
for 
vec.len 
(nonnegative integer) How many 'first few' elements are
displayed of each vector. If 
indent.str 
(string) The indentation string to use. 
nest.lev 
(nonnegative integer) Current nesting level in the recursive
calls to 
give.attr 
(logical) If 
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.
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
.
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.
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)
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)
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 userfacing diagnostic functions, see
pareto_convergence_rate
and pareto_diags
.
ps_convergence_rate(k, ndraws, ...)
ps_convergence_rate(k, ndraws, ...)
k 
paretok values 
ndraws 
number of draws 
... 
unused 
convergence rate
Other helperfunctions:
ps_khat_threshold()
,
ps_min_ss()
,
ps_tail_length()
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
userfacing diagnostic functions, see pareto_khat_threshold
and
pareto_diags
.
ps_khat_threshold(ndraws, ...)
ps_khat_threshold(ndraws, ...)
ndraws 
number of draws 
... 
unused 
threshold
Other helperfunctions:
ps_convergence_rate()
,
ps_min_ss()
,
ps_tail_length()
Given Paretok 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 userfacing diagnostic functions, see
pareto_min_ss
and pareto_diags
.
ps_min_ss(k, ...)
ps_min_ss(k, ...)
k 
pareto k value 
... 
unused 
minimum sample size
Other helperfunctions:
ps_convergence_rate()
,
ps_khat_threshold()
,
ps_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.
ps_tail_length(ndraws, r_eff, ...)
ps_tail_length(ndraws, r_eff, ...)
ndraws 
number of draws 
r_eff 
relative efficiency 
... 
unused 
tail length
Other helperfunctions:
ps_convergence_rate()
,
ps_khat_threshold()
,
ps_min_ss()
Compute quantiles of a sample and return them in a format consistent with other summary functions in the posterior package.
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, ...)
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, ...)
x 
(multiple options) One of:

probs 
(numeric vector) Probabilities in 
na.rm 
(logical) Should 
... 
Arguments passed to individual methods (if applicable) and then
on to 
names 
(logical) Should the result have a 
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
.
mu < extract_variable_matrix(example_draws(), "mu") quantile2(mu)
mu < extract_variable_matrix(example_draws(), "mu") quantile2(mu)
Execute (nearly) arbitrary R expressions that may include rvar
s,
producing a new rvar
.
rdo(expr, dim = NULL, ndraws = NULL)
rdo(expr, dim = NULL, ndraws = NULL)
expr 
(expression) A bare expression that can (optionally) contain

dim 
(integer vector) One or more integers giving the maximal indices
in each dimension to override the dimensions of the 
ndraws 
(positive integer) The number of draws used to construct new
random variables if no 
This function evaluates expr
possibly multiple times, once for each draw of
the rvar
s it contains, then returns a new rvar
representing the output of those
expressions. To identify rvar
s, rdo()
searches the calling environment for any variables
named in expr
for which is_rvar()
evaluates to TRUE
. If expr
contains no rvar
s,
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 rvar
s
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 rvar
s (which should be fast),
using rvar_rng()
, and/or using operations directly on the arrays that back the rvar
s
(via draws_of()
).
An rvar
.
Other rfun:
rfun()
,
rvar_rng()
mu < rdo(rnorm(10, mean = 1:10, sd = 1)) sigma < rdo(rgamma(1, shape = 1, rate = 1)) x < rdo(rnorm(10, mu, sigma)) x
mu < rdo(rnorm(10, mean = 1:10, sd = 1)) sigma < rdo(rgamma(1, shape = 1, rate = 1)) x < rdo(rnorm(10, mu, sigma)) x
draws
objectsRename variables in a draws
object.
rename_variables(.x, ...) ## S3 method for class 'draws' rename_variables(.x, ...)
rename_variables(.x, ...) ## S3 method for class 'draws' rename_variables(.x, ...)
.x 
(draws) A 
... 
One or more expressions, separated by commas, indicating the
variables to rename. The variable names can be unquoted
( 
Returns a draws
object of the same format as .x
, with variables renamed
according to the expressions provided in ...
.
variables
, variables<
, mutate_variables
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)
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)
draws
objectsRepair indices of draws
objects so that iterations, chains, and draws
are continuously and consistently numbered.
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, ...)
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, ...)
x 
(draws) A 
order 
(logical) Should draws be ordered (via 
... 
Arguments passed to individual methods (if applicable). 
A draws
object of the same class as x
.
x < as_draws_array(example_draws()) (x < x[10:5, 3:4, ]) repair_draws(x)
x < as_draws_array(example_draws()) (x < x[10:5, 3:4, ]) repair_draws(x)
draws
objectsResample draws
objects according to provided weights, for example weights
obtained through importance sampling.
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, ...)
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, ...)
x 
(draws) A 
... 
Arguments passed to individual methods (if applicable). 
weights 
(numeric vector) A vector of positive weights of length

method 
(string) The resampling method to use:
Currently, 
ndraws 
(positive integer) The number of draws to be returned. By
default 
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.
A draws
object of the same class as x
.
Kitagawa, G., Monte Carlo Filter and Smoother for NonGaussian Nonlinear ' State Space Models, Journal of Computational and Graphical Statistics, 5(1):125, 1996.
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())
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())
Get names of reserved variables from objects in the posterior package.
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, ...)
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, ...)
x 
(draws) A 
... 
Arguments passed to individual methods (if applicable). 
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.
A character vector of reserved variables used in x
.
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)
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)
Function that create functions that can accept and/or produce rvar
s.
rfun(.f, rvar_args = NULL, rvar_dots = TRUE, ndraws = NULL)
rfun(.f, rvar_args = NULL, rvar_dots = TRUE, ndraws = NULL)
.f 
(multiple options) A function to turn into a function that accepts and/or produces random variables:

rvar_args 
(character vector) The names of the arguments of 
rvar_dots 
(logical) Should dots ( 
ndraws 
(positive integer). The number of draws used to construct new
random variables if no 
This function wraps an existing function (.f
) such that it returns rvar
s 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 rvar
s passed to it, then returns a new
rvar
representing the output of those function evaluations. If the arguments contain no rvar
s,
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 rvar
s
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 rvar
s (which should be fast),
using rvar_rng()
, and/or using operations directly on the arrays that back the rvar
s
(via draws_of()
).
A function with the same argument specification as .f
, but which can accept and return
rvar
s.
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
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
Compute the Rhat convergence diagnostic for a single variable as the maximum of rank normalized splitRhat and rank normalized foldedsplitRhat as proposed in Vehtari et al. (2021).
rhat(x, ...) ## Default S3 method: rhat(x, ...) ## S3 method for class 'rvar' rhat(x, ...)
rhat(x, ...) ## Default S3 method: rhat(x, ...) ## S3 method for class 'rvar' rhat(x, ...)
x 
(multiple options) One of:

... 
Arguments passed to individual methods (if applicable). 
If the input is an array, returns a single numeric value. If any of the draws
is nonfinite, 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.
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and PaulChristian Bürkner (2021). Ranknormalization, folding, and localization: An improved Rhat for assessing convergence of MCMC (with discussion). Bayesian Analysis. 16(2), 667–718. doi:10.1214/20BA1221
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()
mu < extract_variable_matrix(example_draws(), "mu") rhat(mu) d < as_draws_rvars(example_draws("multi_normal")) rhat(d$Sigma)
mu < extract_variable_matrix(example_draws(), "mu") rhat(mu) d < as_draws_rvars(example_draws("multi_normal")) rhat(d$Sigma)
Compute the basic Rhat convergence diagnostic 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
Rhat convergence diagnostic implemented in rhat()
.
rhat_basic(x, ...) ## Default S3 method: rhat_basic(x, split = TRUE, ...) ## S3 method for class 'rvar' rhat_basic(x, split = TRUE, ...)
rhat_basic(x, ...) ## Default S3 method: rhat_basic(x, split = TRUE, ...) ## S3 method for class 'rvar' rhat_basic(x, split = TRUE, ...)
x 
(multiple options) One of:

... 
Arguments passed to individual methods (if applicable). 
split 
(logical) Should the estimate be computed on split chains? The
default is 
If the input is an array, returns a single numeric value. If any of the draws
is nonfinite, 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.
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 PaulChristian Bürkner (2021). Ranknormalization, folding, and localization: An improved Rhat for assessing convergence of MCMC (with discussion). Bayesian Analysis. 16(2), 667–718. doi:10.1214/20BA1221
Other diagnostics:
ess_basic()
,
ess_bulk()
,
ess_quantile()
,
ess_sd()
,
ess_tail()
,
mcse_mean()
,
mcse_quantile()
,
mcse_sd()
,
pareto_diags()
,
pareto_khat()
,
rhat()
,
rhat_nested()
,
rstar()
mu < extract_variable_matrix(example_draws(), "mu") rhat_basic(mu) d < as_draws_rvars(example_draws("multi_normal")) rhat_basic(d$Sigma)
mu < extract_variable_matrix(example_draws(), "mu") rhat_basic(mu) d < as_draws_rvars(example_draws("multi_normal")) rhat_basic(d$Sigma)
Compute the nested Rhat convergence diagnostic for a single variable as proposed in Margossian et al. (2024).
rhat_nested(x, ...) ## Default S3 method: rhat_nested(x, superchain_ids, ...) ## S3 method for class 'rvar' rhat_nested(x, superchain_ids, ...)
rhat_nested(x, ...) ## Default S3 method: rhat_nested(x, superchain_ids, ...) ## S3 method for class 'rvar' rhat_nested(x, superchain_ids, ...)
x 
(multiple options) One of:

... 
Arguments passed to individual methods (if applicable). 
superchain_ids 
(numeric) Vector of length nchains specifying which superchain each chain belongs to. There should be equal numbers of chains in each superchain. All chains within the same superchain are assumed to have been initialized at the same point. 
Nested Rhat is a convergence diagnostic useful when running many short chains. It is calculated on superchains, which are groups of chains that have been initialized at the same 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 3 in Margossian et al. (2024)).
If the input is an array, returns a single numeric value. If any of the draws
is nonfinite, 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.
Charles C. Margossian, Matthew D. Hoffman, Pavel Sountsov, Lionel RiouDurand, Aki Vehtari and Andrew Gelman (2024). Nested Rhat: Assessing the convergence of Markov chain Monte Carlo when running many short chains. Bayesian Analysis. doi:10.1214/24BA1453
Other diagnostics:
ess_basic()
,
ess_bulk()
,
ess_quantile()
,
ess_sd()
,
ess_tail()
,
mcse_mean()
,
mcse_quantile()
,
mcse_sd()
,
pareto_diags()
,
pareto_khat()
,
rhat()
,
rhat_basic()
,
rstar()
mu < extract_variable_matrix(example_draws(), "mu") rhat_nested(mu, superchain_ids = c(1, 1, 2, 2)) d < as_draws_rvars(example_draws("multi_normal")) rhat_nested(d$Sigma, superchain_ids = c(1, 1, 2, 2))
mu < extract_variable_matrix(example_draws(), "mu") rhat_nested(mu, superchain_ids = c(1, 1, 2, 2)) d < as_draws_rvars(example_draws("multi_normal")) rhat_nested(d$Sigma, superchain_ids = c(1, 1, 2, 2))
The rstar()
function generates a measure of convergence for MCMC draws
based on whether it is possible to determine the Markov chain that generated
a draw with probability greater than chance. To do so, it fits a machine
learning classifier to a training set of MCMC draws and evaluates its
predictive accuracy on a testing set: giving the ratio of accuracy to
predicting a chain uniformly at random.
rstar( x, split = TRUE, uncertainty = FALSE, method = "rf", hyperparameters = NULL, training_proportion = 0.7, nsimulations = 1000, ... )
rstar( x, split = TRUE, uncertainty = FALSE, method = "rf", hyperparameters = NULL, training_proportion = 0.7, nsimulations = 1000, ... )
x 
(draws) A 
split 
(logical) Should the estimate be computed on split chains? The
default is 
uncertainty 
(logical). Indicates whether to provide a vector of R*
values representing uncertainty in the calculated value (if 
method 
(string) The machine learning classifier to use (must be
available in the caret package). The default is 
hyperparameters 
(named list) Hyperparameter settings passed to the classifier.
The default for the random forest classifier ( 
training_proportion 
(positive real) The proportion (in 
nsimulations 
(positive integer) The number of R* values in the
returned vector if 
... 
Other arguments passed to 
The rstar()
function provides a measure of MCMC convergence based
on whether it is possible to determine the chain that generated a
particular draw with a probability greater than chance. To do so, it fits a
machine learning classifier to a subset of the original MCMC draws (the
training set) and evaluates its predictive accuracy on the remaining draws
(the testing set). If predictive accuracy exceeds chance (i.e. predicting
the chain that generated a draw uniformly at random), the diagnostic
measure R* will be above 1, indicating that convergence has yet to occur.
This statistic is recently developed, and it is currently unclear what is a
reasonable threshold for diagnosing convergence.
The statistic, R*, is stochastic, meaning that each time the test is run,
unless the random seed is fixed, it will generally produce a different
result. To minimize the implications of this stochasticity, it is
recommended to repeatedly run this function to calculate a distribution of
R*; alternatively, an approximation to this distribution can be obtained by
setting uncertainty = TRUE
, although this approximation of uncertainty
will generally have a lower mean.
By default, a random forest classifier is used (method = "rf"
), which tends
to perform best for target distributions of around 4 dimensions and above.
For lower dimensional targets, gradient boosted models (called via
method = "gbm"
) tend to have a higher classification accuracy. On a given
MCMC sample, it is recommended to try both of these classifiers.
A numeric vector of length 1 (by default) or length nsimulations
(if uncertainty = TRUE
).
Ben Lambert, Aki Vehtari (2020) R*: A robust MCMC convergence
diagnostic with uncertainty using gradientboosted machines.
arXiv preprint arXiv:2003.07900
.
Other diagnostics:
ess_basic()
,
ess_bulk()
,
ess_quantile()
,
ess_sd()
,
ess_tail()
,
mcse_mean()
,
mcse_quantile()
,
mcse_sd()
,
pareto_diags()
,
pareto_khat()
,
rhat()
,
rhat_basic()
,
rhat_nested()
if (require("caret", quietly = TRUE) && require("randomForest", quietly = TRUE)) { x < example_draws("eight_schools") print(rstar(x)) print(rstar(x, split = FALSE)) print(rstar(x, method = "gbm")) # can pass additional arguments to methods print(rstar(x, method = "gbm", verbose = FALSE)) # with uncertainty, returns a vector of R* values hist(rstar(x, uncertainty = TRUE)) hist(rstar(x, uncertainty = TRUE, nsimulations = 100)) # can use other classification methods in caret library print(rstar(x, method = "knn")) }
if (require("caret", quietly = TRUE) && require("randomForest", quietly = TRUE)) { x < example_draws("eight_schools") print(rstar(x)) print(rstar(x, split = FALSE)) print(rstar(x, method = "gbm")) # can pass additional arguments to methods print(rstar(x, method = "gbm", verbose = FALSE)) # with uncertainty, returns a vector of R* values hist(rstar(x, uncertainty = TRUE)) hist(rstar(x, uncertainty = TRUE, nsimulations = 100)) # can use other classification methods in caret library print(rstar(x, method = "knn")) }
Random variables backed by arrays of arbitrary dimension
rvar( x = double(), dim = NULL, dimnames = NULL, nchains = NULL, with_chains = FALSE )
rvar( x = double(), dim = NULL, dimnames = NULL, nchains = NULL, with_chains = FALSE )
x 
(multiple options) The object to convert to an

dim 
(integer vector) One or more integers giving the maximal indices
in each dimension to override the dimensions of the 
dimnames 
(list) Character vectors giving the names in each dimension
to override the names of the dimensions of the 
nchains 
(positive integer) The number of chains. The if 
with_chains 
(logical) Does 
The "rvar"
class internally represents random variables as arrays of arbitrary
dimension, where the first dimension is used to index draws from the distribution.
Most mathematical operators and functions are supported, including efficient matrix
multiplication and vector and arraystyle indexing. The intent is that an rvar
works as closely as possible to how a base vector/matrix/array does, with a few
differences:
The default behavior when subsetting is not to drop extra dimensions (i.e.
the default drop
argument for [
is FALSE
, not TRUE
).
Rather than base Rstyle recycling, rvar
s use a limited form of broadcasting:
if an operation is being performed on two vectors with different size of the same
dimension, the smaller vector will be recycled up to the size of the larger one
along that dimension so long as it has size 1.
For functions that expect base numeric arrays and for which rvar
s cannot be
used directly as arguments, you can use rfun()
or rdo()
to translate your
code into code that executes across draws from one or more random variables
and returns a random variable as output. Typically rdo()
offers the most
straightforward translation.
As rfun()
and rdo()
incur some performance cost, you can also operate directly
on the underlying array using the draws_of()
function. To reuse existing
random number generator functions to efficiently create rvar
s, use rvar_rng()
.
An object of class "rvar"
representing a random variable.
as_rvar()
to convert objects to rvar
s. See rdo()
, rfun()
, and
rvar_rng()
for higherlevel interfaces for creating rvar
s.
set.seed(1234) # To create a "scalar" `rvar`, pass a onedimensional array or a vector # whose length (here `4000`) is the desired number of draws: x < rvar(rnorm(4000, mean = 1, sd = 1)) x # Create random vectors by adding an additional dimension: n < 4 # length of output vector x < rvar(array(rnorm(4000 * n, mean = rep(1:n, each = 4000), sd = 1), dim = c(4000, n))) x # Create a random matrix: rows < 4 cols < 3 x < rvar(array(rnorm(4000 * rows * cols, mean = 1, sd = 1), dim = c(4000, rows, cols))) x # If the input sample comes from multiple chains, we can indicate that using the # nchains argument (here, 1000 draws each from 4 chains): x < rvar(rnorm(4000, mean = 1, sd = 1), nchains = 4) x # Or if the input sample has chain information as its second dimension, we can # use with_chains to create the rvar x < rvar(array(rnorm(4000, mean = 1, sd = 1), dim = c(1000, 4)), with_chains = TRUE) x
set.seed(1234) # To create a "scalar" `rvar`, pass a onedimensional array or a vector # whose length (here `4000`) is the desired number of draws: x < rvar(rnorm(4000, mean = 1, sd = 1)) x # Create random vectors by adding an additional dimension: n < 4 # length of output vector x < rvar(array(rnorm(4000 * n, mean = rep(1:n, each = 4000), sd = 1), dim = c(4000, n))) x # Create a random matrix: rows < 4 cols < 3 x < rvar(array(rnorm(4000 * rows * cols, mean = 1, sd = 1), dim = c(4000, rows, cols))) x # If the input sample comes from multiple chains, we can indicate that using the # nchains argument (here, 1000 draws each from 4 chains): x < rvar(rnorm(4000, mean = 1, sd = 1), nchains = 4) x # Or if the input sample has chain information as its second dimension, we can # use with_chains to create the rvar x < rvar(array(rnorm(4000, mean = 1, sd = 1), dim = c(1000, 4)), with_chains = TRUE) x
Returns an rvar
obtained by applying a function to margins of an array or rvar
.
Acts like apply()
, except that the function supplied (.f
) should return an rvar
,
and the final result is always an rvar
.
rvar_apply(.x, .margin, .f, ...)
rvar_apply(.x, .margin, .f, ...)
.x 
An array or an 
.margin 
(multiple options) The subscripts which the function will be applied over:

.f 
(function) The function to be applied. The function 
... 
Optional arguments passed to 
This function acts much like apply()
, except that the function passed to it (.f
)
must return rvar
s, and the result is simplified into an rvar
. Unlike
apply()
, it also keeps the dimensions of the returned values along each margin,
rather than simplifying each margin to a vector, and if the results of .f
do
not all have the same dimensions, it applies the rvar
broadcasting rules to
bind results together rather than using vector recycling.
If you wish to apply functions over rvar
s where the result is not intended to
be simplified into an rvar
, you can use the standard apply()
, lapply()
,
sapply()
, or vapply()
functions.
An rvar
.
If the result of each call to .f
returns an rvar
of dimension d
after
being broadcast to a common shape, then rvar_apply()
returns an rvar
of
dimension c(d, dim(.x)[.margin])
. If the last dimension of the result would
be 1
, it is dropped (other dimensions equal to 1
are retained). If d
is
0
, the result has length 0
but not necessarily the 'correct' dimension.
as_rvar()
to convert objects to rvar
s. See rdo()
, rfun()
, and
rvar_rng()
for higherlevel interfaces for creating rvar
s.
set.seed(3456) x < rvar_rng(rnorm, 24, mean = 1:24) dim(x) < c(2,3,4) # we can find the distributions of marginal means of the above array # using rvar_mean along with rvar_apply rvar_apply(x, 1, rvar_mean) rvar_apply(x, 2:3, rvar_mean)
set.seed(3456) x < rvar_rng(rnorm, 24, mean = 1:24) dim(x) < c(2,3,4) # we can find the distributions of marginal means of the above array # using rvar_mean along with rvar_apply rvar_apply(x, 1, rvar_mean) rvar_apply(x, 2:3, rvar_mean)
Random variables backed by factorlike arrays of arbitrary dimension.
rvar_factor( x = factor(), dim = NULL, dimnames = NULL, nchains = NULL, with_chains = FALSE, ... ) rvar_ordered( x = ordered(NULL), dim = NULL, dimnames = NULL, nchains = NULL, with_chains = FALSE, ... )
rvar_factor( x = factor(), dim = NULL, dimnames = NULL, nchains = NULL, with_chains = FALSE, ... ) rvar_ordered( x = ordered(NULL), dim = NULL, dimnames = NULL, nchains = NULL, with_chains = FALSE, ... )
x 
(multiple options) The object to convert to an

dim 
(integer vector) One or more integers giving the maximal indices
in each dimension to override the dimensions of the 
dimnames 
(list) Character vectors giving the names in each dimension
to override the names of the dimensions of the 
nchains 
(positive integer) The number of chains. The if 
with_chains 
(logical) Does 
... 
Arguments passed on to

A subtype of rvar()
that represents a (possibly multidimensional) sample of
a factor or an ordered factor. It is otherwise very similar to the basic rvar()
:
it is backed by a multidimensional array with draws as the first dimension.
The primary difference is that the backing array has class "factor"
(for rvar_factor()
)
or c("ordered", "factor")
(for rvar_ordered()
). If you
pass a factor or ordered factor to rvar()
it will automatically return
an object with the classes "rvar_factor"
or c("rvar_ordered", "rvar_factor")
.
See rvar()
for more details on the internals of the random variable datatype.
An object of class "rvar_factor"
representing a factor
like random variable.
as_rvar_factor()
to convert objects to rvar_factor
s. See rdo()
, rfun()
, and
rvar_rng()
for higherlevel interfaces for creating rvar
s.
set.seed(1234) # To create a "scalar" `rvar_factor`, pass a onedimensional array or a vector # whose length (here `4000`) is the desired number of draws: x < rvar(sample(c("a","a","a","b","c"), 4000, replace = TRUE)) x # Create random vectors by adding an additional dimension: x_array < array(c( sample(c("a","a","a","b","c"), 4000, replace = TRUE), sample(c("a","a","b","c","c"), 4000, replace = TRUE), sample(c("b","b","b","b","c"), 4000, replace = TRUE), sample(c("d","d","b","b","c"), 4000, replace = TRUE) ), dim = c(4000, 4)) rvar_factor(x_array) # You can also create ordered factors rvar_ordered(x_array) # arguments of factor() and ordered() are passed down by the constructor # e.g. we can reorder levels of an ordered factor: rvar_ordered(x_array, levels = c("d","c","b","a")) # Unlike base factors, rvar factors can be matrices or arrays: rvar_factor(x_array, dim = c(2, 2)) # If the input to rvar_factor() is an array with a `"levels"` attribute, it # will use those as the levels of the factor y_array < t(array(rbinom(3000, 1, c(0.1, 0.5, 0.9)) + 1, dim = c(3, 1000))) rvar(y_array) # with levels attr(y_array, "levels") = c("a", "b") rvar_factor(y_array)
set.seed(1234) # To create a "scalar" `rvar_factor`, pass a onedimensional array or a vector # whose length (here `4000`) is the desired number of draws: x < rvar(sample(c("a","a","a","b","c"), 4000, replace = TRUE)) x # Create random vectors by adding an additional dimension: x_array < array(c( sample(c("a","a","a","b","c"), 4000, replace = TRUE), sample(c("a","a","b","c","c"), 4000, replace = TRUE), sample(c("b","b","b","b","c"), 4000, replace = TRUE), sample(c("d","d","b","b","c"), 4000, replace = TRUE) ), dim = c(4000, 4)) rvar_factor(x_array) # You can also create ordered factors rvar_ordered(x_array) # arguments of factor() and ordered() are passed down by the constructor # e.g. we can reorder levels of an ordered factor: rvar_ordered(x_array, levels = c("d","c","b","a")) # Unlike base factors, rvar factors can be matrices or arrays: rvar_factor(x_array, dim = c(2, 2)) # If the input to rvar_factor() is an array with a `"levels"` attribute, it # will use those as the levels of the factor y_array < t(array(rbinom(3000, 1, c(0.1, 0.5, 0.9)) + 1, dim = c(3, 1000))) rvar(y_array) # with levels attr(y_array, "levels") = c("a", "b") rvar_factor(y_array)
A version of ifelse()
that returns an rvar
.
rvar_ifelse(test, yes, no)
rvar_ifelse(test, yes, no)
test 
(logical 
yes 
( 
no 
( 
An rvar
with the common type of yes
and no
(as determined by
vctrs::vec_cast_common()
) and a shape determined by broadcasting test
,
yes
, and no
to a common shape (see the section on broadcasting rules in
vignette("rvar")
). For every element of draws_of(test)
, the corresponding
element of draws_of(yes)
or draws_of(no)
is placed into the result,
depending on whether the element of test
is TRUE
or FALSE
.
x < rvar(1:4) y < rvar(5:8) i < rvar(c(TRUE,FALSE,TRUE,FALSE)) z < rvar_ifelse(i, x, y) z draws_of(z)
x < rvar(1:4) y < rvar(5:8) i < rvar(c(TRUE,FALSE,TRUE,FALSE)) z < rvar_ifelse(i, x, y) z draws_of(z)
Compute special value predicates (checking for finite / infinite values, NaN
, and NA
)
on all draws within a random variable, returning a random variable.
rvar_is_finite(x) rvar_is_infinite(x) rvar_is_nan(x) rvar_is_na(x)
rvar_is_finite(x) rvar_is_infinite(x) rvar_is_nan(x) rvar_is_na(x)
x 
(rvar) An 
These functions return a new rvar
that is the result of applying
is.finite()
, is.infinite()
, is.nan()
, or is.na()
to every draw
in the input random variable.
A logical rvar
of the same length as the input.
rvarsummariesoverdraws for summary functions across draws, including
implementations of is.finite()
, is.infinite()
, is.nan()
, and is.na()
for
rvar
s.
Other rvarsummaries:
rvarsummariesoverdraws
,
rvarsummarieswithindraws
x < rvar(c(1, Inf, Inf, NaN, NA)) x rvar_is_finite(x) rvar_is_infinite(x) rvar_is_nan(x) rvar_is_na(x)
x < rvar(c(1, Inf, Inf, NaN, NA)) x rvar_is_finite(x) rvar_is_infinite(x) rvar_is_nan(x) rvar_is_na(x)
Specialized alternative to rdo()
or rfun()
for creating rvar
s from
existing randomnumber generator functions (such as rnorm()
, rbinom()
, etc).
rvar_rng(.f, n, ..., ndraws = NULL)
rvar_rng(.f, n, ..., ndraws = NULL)
.f 
(function) A function (or string naming a function) representing a
randomnumber generating function that follows the pattern of base random
number generators (like

n 
(positive integer) The length of the output 
... 
Arguments passed to 
ndraws 
(positive integer) The number of draws used to construct the
returned random variable if no 
This function unwraps the arrays underlying the input rvar
s in
...
and then passes them to .f
, relying on the vectorization of .f
to evaluate it across draws from the input rvar
s. This is why the arguments
of .f
must be vectorized. It asks for n
times the number of draws
in the input rvar
s (or ndraws
if none are given) draws from the
random number generator .f
, then reshapes the output from .f
into an
rvar
with length n
.
rvar_rng()
is a fast alternative to rdo()
or rfun()
, but you must
ensure that .f
satisfies the preconditions described above for the result
to be correct. Most base random number generators satisfy these conditions.
It is advisable to test against rdo()
or rfun()
(which should be correct,
but slower) if you are uncertain.
A singledimensional rvar
of length n
.
mu < rvar_rng(rnorm, 10, mean = 1:10, sd = 1) sigma < rvar_rng(rgamma, 1, shape = 1, rate = 1) x < rvar_rng(rnorm, 10, mu, sigma) x
mu < rvar_rng(rnorm, 10, mean = 1:10, sd = 1) sigma < rvar_rng(rgamma, 1, shape = 1, rate = 1) x < rvar_rng(rnorm, 10, mu, sigma) x
The probability density function (density()
), cumulative distribution
function (cdf()
), and quantile function / inverse CDF (quantile()
) of
an rvar
.
## S3 method for class 'rvar' density(x, at, ...) ## S3 method for class 'rvar_factor' density(x, at, ...) ## S3 method for class 'rvar' cdf(x, q, ...) ## S3 method for class 'rvar_factor' cdf(x, q, ...) ## S3 method for class 'rvar_ordered' cdf(x, q, ...) ## S3 method for class 'rvar' quantile(x, probs, ...) ## S3 method for class 'rvar_factor' quantile(x, probs, ...) ## S3 method for class 'rvar_ordered' quantile(x, probs, ...)
## S3 method for class 'rvar' density(x, at, ...) ## S3 method for class 'rvar_factor' density(x, at, ...) ## S3 method for class 'rvar' cdf(x, q, ...) ## S3 method for class 'rvar_factor' cdf(x, q, ...) ## S3 method for class 'rvar_ordered' cdf(x, q, ...) ## S3 method for class 'rvar' quantile(x, probs, ...) ## S3 method for class 'rvar_factor' quantile(x, probs, ...) ## S3 method for class 'rvar_ordered' quantile(x, probs, ...)
x 
(rvar) An 
... 
Additional arguments passed onto underlying methods:

q , at

(numeric vector) One or more quantiles. 
probs 
(numeric vector) One or more probabilities in 
If x
is a scalar rvar
, returns a vector of the same length as the input
(q
, at
, or probs
) containing values from the corresponding function
of the given rvar
.
If x
has length greater than 1, returns an array with dimensions
c(length(y), dim(x))
where y
is q
, at
, or probs
, where each
result[i,...]
is the value of the corresponding function,f(y[i])
, for
the corresponding cell in the input array, x[...]
.
set.seed(1234) x = rvar(rnorm(100)) density(x, seq(2, 2, length.out = 10)) cdf(x, seq(2, 2, length.out = 10)) quantile(x, ppoints(10)) x2 = c(rvar(rnorm(100, mean = 0.5)), rvar(rnorm(100, mean = 0.5))) density(x2, seq(2, 2, length.out = 10)) cdf(x2, seq(2, 2, length.out = 10)) quantile(x2, ppoints(10))
set.seed(1234) x = rvar(rnorm(100)) density(x, seq(2, 2, length.out = 10)) cdf(x, seq(2, 2, length.out = 10)) quantile(x, ppoints(10)) x2 = c(rvar(rnorm(100, mean = 0.5)), rvar(rnorm(100, mean = 0.5))) density(x2, seq(2, 2, length.out = 10)) cdf(x2, seq(2, 2, length.out = 10)) quantile(x2, ppoints(10))
Matrix multiplication of random variables.
x %**% y ## S3 method for class 'rvar' matrixOps(x, y)
x %**% y ## S3 method for class 'rvar' matrixOps(x, y)
x 
(multiple options) The object to be postmultiplied by If a vector is used, it is treated as a row vector. 
y 
(multiple options) The object to be premultiplied by If a vector is used, it is treated as a column vector. 
If x
or y
are vectors, they are converted into matrices prior to multiplication, with x
converted to a row vector and y
to a column vector. Numerics and logicals can be multiplied
by rvar
s 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)
.
In R >= 4.3, you can also use %*%
in place of %**%
for matrix multiplication
of rvar
s. In R < 4.3, S3 classes cannot properly override %*%
, so
you must use %**%
for matrix multiplication of rvar
s.
An rvar
representing the matrix product of x
and y
.
# d has mu (mean vector of length 3) and Sigma (3x3 covariance matrix) d < as_draws_rvars(example_draws("multi_normal")) d$Sigma # trivial example: multiplication by a nonrandom matrix d$Sigma %**% diag(1:3) # Decompose Sigma into R s.t. R'R = Sigma ... R < chol(d$Sigma) # ... and recreate Sigma using matrix multiplication t(R) %**% R
# d has mu (mean vector of length 3) and Sigma (3x3 covariance matrix) d < as_draws_rvars(example_draws("multi_normal")) d$Sigma # trivial example: multiplication by a nonrandom matrix d$Sigma %**% diag(1:3) # Decompose Sigma into R s.t. R'R = Sigma ... R < chol(d$Sigma) # ... and recreate Sigma using matrix multiplication t(R) %**% R
Operations for slicing rvar
s and replacing parts of rvar
s.
## S3 method for class 'rvar' x[[i, ...]] ## S3 replacement method for class 'rvar' x[[i, ...]] < value ## S3 method for class 'rvar' x[..., drop = FALSE] ## S3 replacement method for class 'rvar' x[i, ...] < value
## S3 method for class 'rvar' x[[i, ...]] ## S3 replacement method for class 'rvar' x[[i, ...]] < value ## S3 method for class 'rvar' x[..., drop = FALSE] ## S3 replacement method for class 'rvar' x[i, ...] < value
x 
an 
i , ...

indices; see Details. 
value 
( 
drop 
(logical) Should singular dimensions be dropped when slicing
array 
The rvar
slicing operators ([
and [[
) attempt to implement the same
semantics as the base array slicing operators. There are some
exceptions; most notably, rvar
slicing defaults to drop = FALSE
instead
of drop = TRUE
.
[[
The [[
operator extracts (or replaces) single elements. It always
returns (or replaces) a scalar (length1) rvar
.
The x[[i,...]]
operator can be used as follows:
x[[<numeric>]]
for scalar numeric i
: gives the i
th element of x
. If x
is
multidimensional (i.e. length(dim(x)) > 1
), extra dimensions are ignored
when indexing. For example, if x
is a $6 \times 2$
rvar
array, the
7th element, x[[7]]
, will be the first element of the second column, x[1,2]
.
x[[<numeric rvar>]]
for scalar numeric rvar
i
: a generalization of indexing when
i
is a scalar numeric. Within each draw of x
, selects the element
corresponding to the value of i
within that same draw.
x[[<character>]]
for scalar character i
: gives the element of x
with name
equal to i
. Unlike with base arrays, does not work with
multidimensional rvar
s.
x[[i_1,i_2,...,i_n]]
for scalar numeric or character i_1
, i_2
, etc.
Must provide exactly the same number of indices as dimensions in x
. Selects
the element at the corresponding position in the rvar
by number and/or
dimname (as a string).
[
The [
operator extracts (or replaces) multiple elements. It always returns
(or replaces) a possiblymultidimensional rvar
.
The x[i,...]
operator can be used as follows:
x[<logical>]
for vector logical i
: i
is recycled to the same length as x
,
ignoring multiple dimensions in x
, then an rvar
vector is returned
containing the elements in x
where i
is TRUE
.
x[<logical rvar>]
for scalar logical rvar
i
: returns an rvar
the same shape
as x
containing only those draws where i
is TRUE
.
x[<numeric>]
for vector numeric i
: an rvar
vector is returned
containing the i
th elements of x
, ignoring dimensions.
x[<matrix>]
for numeric matrix i
, where ncol(i) == length(dim(x))
: each row
of i
should give the multidimensional index for a single element in x
. The
result is an rvar
vector of length nrow(i)
containing elements of x
selected by each row of i
.
x[i_1,i_2,...,i_n]
for vector numeric, character, or logical i_1
,
i_2
, etc. Returns a slice of x
containing all elements from the dimensions
specified in i_1
, i_2
, etc. If an argument is left empty, all elements
from that dimension are included. Unlike base arrays, trailing dimensions
can be omitted entirely and will still be selected; for example, if x
has
three dimensions, both x[1,,]
and x[1,]
can be used to create a
slice that includes all elements from the last two dimensions. Unlike base
arrays, [
defaults to drop = FALSE
, so results retain the same number of
dimensions as x
.
x < rvar(array(1:24, dim = c(4,2,3))) dimnames(x) < list(c("a","b"), c("d","e","f")) x ## Slicing single elements # x[[<numeric>]] x[[2]] # x[[<numeric rvar>]] # notice the draws of x[1:4]... draws_of(x[1:4]) x[[rvar(c(1,3,4,4))]] # ... x[[rvar(c(1,3,4,4))]] creates a mixures of those draws draws_of(x[[rvar(c(1,3,4,4))]]) # x[[i_1,i_2,...]] x[[2,"e"]] ## Slicing multiple elements # x[<logical>] x[c(TRUE,TRUE,FALSE)] # x[<logical rvar>] # select every other draw x[rvar(c(TRUE,FALSE,TRUE,FALSE))] # x[<numeric>] x[1:3] # x[<matrix>] x[rbind( c(1,2), c(1,3), c(2,2) )] # x[i_1,i_2,...,i_n] x[1,] x[1,2:3] x[,2:3]
x < rvar(array(1:24, dim = c(4,2,3))) dimnames(x) < list(c("a","b"), c("d","e","f")) x ## Slicing single elements # x[[<numeric>]] x[[2]] # x[[<numeric rvar>]] # notice the draws of x[1:4]... draws_of(x[1:4]) x[[rvar(c(1,3,4,4))]] # ... x[[rvar(c(1,3,4,4))]] creates a mixures of those draws draws_of(x[[rvar(c(1,3,4,4))]]) # x[[i_1,i_2,...]] x[[2,"e"]] ## Slicing multiple elements # x[<logical>] x[c(TRUE,TRUE,FALSE)] # x[<logical rvar>] # select every other draw x[rvar(c(TRUE,FALSE,TRUE,FALSE))] # x[<numeric>] x[1:3] # x[<matrix>] x[rbind( c(1,2), c(1,3), c(2,2) )] # x[i_1,i_2,...,i_n] x[1,] x[1,2:3] x[,2:3]
Compute summaries within elements of an rvar
and over draws of each element,
producing an array of the same shape as the input random variable (except in
the case of range()
, see Details).
E(x, ...) ## S3 method for class 'rvar' mean(x, ...) Pr(x, ...) ## Default S3 method: Pr(x, ...) ## S3 method for class 'logical' Pr(x, ...) ## S3 method for class 'rvar' Pr(x, ...) ## S3 method for class 'rvar' median(x, ...) ## S3 method for class 'rvar' min(x, ...) ## S3 method for class 'rvar' max(x, ...) ## S3 method for class 'rvar' sum(x, ...) ## S3 method for class 'rvar' prod(x, ...) ## S3 method for class 'rvar' all(x, ...) ## S3 method for class 'rvar' any(x, ...) ## S3 method for class 'rvar' Summary(...) ## S3 method for class 'rvar' variance(x, ...) var(x, ...) ## Default S3 method: var(x, ...) ## S3 method for class 'rvar' var(x, ...) sd(x, ...) ## Default S3 method: sd(x, ...) ## S3 method for class 'rvar' sd(x, ...) mad(x, ...) ## Default S3 method: mad(x, ...) ## S3 method for class 'rvar' mad(x, ...) ## S3 method for class 'rvar_ordered' mad(x, ...) ## S3 method for class 'rvar' range(x, ...) ## S3 method for class 'rvar' is.finite(x) ## S3 method for class 'rvar' is.infinite(x) ## S3 method for class 'rvar' is.nan(x) ## S3 method for class 'rvar' is.na(x)
E(x, ...) ## S3 method for class 'rvar' mean(x, ...) Pr(x, ...) ## Default S3 method: Pr(x, ...) ## S3 method for class 'logical' Pr(x, ...) ## S3 method for class 'rvar' Pr(x, ...) ## S3 method for class 'rvar' median(x, ...) ## S3 method for class 'rvar' min(x, ...) ## S3 method for class 'rvar' max(x, ...) ## S3 method for class 'rvar' sum(x, ...) ## S3 method for class 'rvar' prod(x, ...) ## S3 method for class 'rvar' all(x, ...) ## S3 method for class 'rvar' any(x, ...) ## S3 method for class 'rvar' Summary(...) ## S3 method for class 'rvar' variance(x, ...) var(x, ...) ## Default S3 method: var(x, ...) ## S3 method for class 'rvar' var(x, ...) sd(x, ...) ## Default S3 method: sd(x, ...) ## S3 method for class 'rvar' sd(x, ...) mad(x, ...) ## Default S3 method: mad(x, ...) ## S3 method for class 'rvar' mad(x, ...) ## S3 method for class 'rvar_ordered' mad(x, ...) ## S3 method for class 'rvar' range(x, ...) ## S3 method for class 'rvar' is.finite(x) ## S3 method for class 'rvar' is.infinite(x) ## S3 method for class 'rvar' is.nan(x) ## S3 method for class 'rvar' is.na(x)
x 
(rvar) An 
... 
Further arguments passed to underlying functions (e.g.,

Summaries include expectations (E()
or mean()
), probabilities (Pr()
),
medians (median()
), spread (var()
, variance()
, sd()
, mad()
), sums and
products (sum()
, prod()
), extrema and ranges (min()
, max()
, range()
),
logical summaries (all()
, any()
), and special value predicates (is.finite()
,
is.infinite()
, is.nan()
, is.na()
).
Unless otherwise stated, these functions return a numeric array with the same shape
(same dimensions) as the input rvar
, x
.
range(x)
returns an array with dimensions c(2, dim(x))
, where the last
dimension contains the minimum and maximum values.
is.infinite(x)
, is.nan(x)
, and is.na(x)
return logical arrays, where each
element is TRUE
if any draws in its corresponding element in x
match
the predicate. Each elements in the result of is.finite(x)
is TRUE
if
all draws in the corresponding element in x
are finite.
Both E()
, mean()
, and Pr()
return the means of each element in the input.
Pr()
additionally checks that the provided rvar
is a logical variable (hence, taking its expectation results in a probability).
For consistency, E()
and Pr()
are also defined for base arrays so that
they can be used as summary functions in summarise_draws()
.
A numeric or logical vector with the same dimensions as the given random variable, where
each entry in the vector is the mean, median, or variance of the corresponding entry in x
.
rvarsummarieswithindraws for summary functions within draws. rvardist for density, CDF, and quantile functions of random variables.
Other rvarsummaries:
rvarsummarieswithindraws
,
rvar_is_finite()
set.seed(5678) x = rvar_rng(rnorm, 4, mean = 1:4, sd = 2) # These should all be ~= c(1, 2, 3, 4) E(x) mean(x) median(x) # This ... Pr(x < 1.5) # ... should be about the same as this: pnorm(1.5, mean = 1:4, sd = 2)
set.seed(5678) x = rvar_rng(rnorm, 4, mean = 1:4, sd = 2) # These should all be ~= c(1, 2, 3, 4) E(x) mean(x) median(x) # This ... Pr(x < 1.5) # ... should be about the same as this: pnorm(1.5, mean = 1:4, sd = 2)
Compute summaries of random variables over array elements and within draws,
producing a new random variable of length 1 (except in the case of
rvar_range()
, see Details).
rvar_mean(..., na.rm = FALSE) rvar_median(..., na.rm = FALSE) rvar_sum(..., na.rm = FALSE) rvar_prod(..., na.rm = FALSE) rvar_min(..., na.rm = FALSE) rvar_max(..., na.rm = FALSE) rvar_sd(..., na.rm = FALSE) rvar_var(..., na.rm = FALSE) rvar_mad(..., constant = 1.4826, na.rm = FALSE) rvar_range(..., na.rm = FALSE) rvar_quantile(..., probs, names = FALSE, na.rm = FALSE) rvar_all(..., na.rm = FALSE) rvar_any(..., na.rm = FALSE)
rvar_mean(..., na.rm = FALSE) rvar_median(..., na.rm = FALSE) rvar_sum(..., na.rm = FALSE) rvar_prod(..., na.rm = FALSE) rvar_min(..., na.rm = FALSE) rvar_max(..., na.rm = FALSE) rvar_sd(..., na.rm = FALSE) rvar_var(..., na.rm = FALSE) rvar_mad(..., constant = 1.4826, na.rm = FALSE) rvar_range(..., na.rm = FALSE) rvar_quantile(..., probs, names = FALSE, na.rm = FALSE) rvar_all(..., na.rm = FALSE) rvar_any(..., na.rm = FALSE)
... 
(rvar) One or more 
na.rm 
(logical) Should 
constant 
(scalar real) For 
probs 
(numeric vector) For 
names 
(logical) For 
These functions compute statistics within each draw of the random variable. For summaries over draws (such as expectations), see rvarsummariesoverdraws.
Each function defined here corresponds to the base function of the same name
without the rvar_
prefix (e.g., rvar_mean()
calls mean()
under the hood, etc).
An rvar
of length 1 (for range()
, length 2; for quantile()
, length
equal to length(probs)
) with the same number
of draws as the input rvar(s) containing the summary statistic computed within
each draw of the input rvar(s).
rvarsummariesoverdraws for summary functions across draws (e.g. expectations). rvardist for density, CDF, and quantile functions of random variables.
Other rvarsummaries:
rvarsummariesoverdraws
,
rvar_is_finite()
set.seed(5678) x = rvar_rng(rnorm, 4, mean = 1:4, sd = 2) # These will give similar results to mean(1:4), # median(1:4), sum(1:4), prod(1:4), etc rvar_mean(x) rvar_median(x) rvar_sum(x) rvar_prod(x) rvar_range(x) rvar_quantile(x, probs = c(0.25, 0.5, 0.75), names = TRUE)
set.seed(5678) x = rvar_rng(rnorm, 4, mean = 1:4, sd = 2) # These will give similar results to mean(1:4), # median(1:4), sum(1:4), prod(1:4), etc rvar_mean(x) rvar_median(x) rvar_sum(x) rvar_prod(x) rvar_range(x) rvar_quantile(x, probs = c(0.25, 0.5, 0.75), names = TRUE)
Split chains by halving the number of iterations per chain and doubling the number of chains.
split_chains(x, ...)
split_chains(x, ...)
x 
(draws) A 
... 
Arguments passed to individual methods (if applicable). 
A draws
object of the same class as x
.
x < example_draws() niterations(x) nchains(x) x < split_chains(x) niterations(x) nchains(x)
x < example_draws() niterations(x) nchains(x) x < split_chains(x) niterations(x) nchains(x)
draws
objectsSubset draws
objects by variables, iterations, chains, and draws indices.
subset_draws(x, ...) ## S3 method for class 'draws_matrix' subset_draws( x, variable = NULL, iteration = NULL, chain = NULL, draw = NULL, regex = FALSE, unique = TRUE, exclude = FALSE, scalar = FALSE, ... ) ## S3 method for class 'draws_array' subset_draws( x, variable = NULL, iteration = NULL, chain = NULL, draw = NULL, regex = FALSE, unique = TRUE, exclude = FALSE, scalar = FALSE, ... ) ## S3 method for class 'draws_df' subset_draws( x, variable = NULL, iteration = NULL, chain = NULL, draw = NULL, regex = FALSE, unique = TRUE, exclude = FALSE, scalar = FALSE, ... ) ## S3 method for class 'draws_list' subset_draws( x, variable = NULL, iteration = NULL, chain = NULL, draw = NULL, regex = FALSE, unique = TRUE, exclude = FALSE, scalar = FALSE, ... ) ## S3 method for class 'draws_rvars' subset_draws( x, variable = NULL, iteration = NULL, chain = NULL, draw = NULL, regex = FALSE, unique = TRUE, exclude = FALSE, scalar = FALSE, ... ) ## S3 method for class 'rvar' subset_draws(x, variable = NULL, ...) ## S3 method for class 'draws' subset(x, ...)
subset_draws(x, ...) ## S3 method for class 'draws_matrix' subset_draws( x, variable = NULL, iteration = NULL, chain = NULL, draw = NULL, regex = FALSE, unique = TRUE, exclude = FALSE, scalar = FALSE, ... ) ## S3 method for class 'draws_array' subset_draws( x, variable = NULL, iteration = NULL, chain = NULL, draw = NULL, regex = FALSE, unique = TRUE, exclude = FALSE, scalar = FALSE, ... ) ## S3 method for class 'draws_df' subset_draws( x, variable = NULL, iteration = NULL, chain = NULL, draw = NULL, regex = FALSE, unique = TRUE, exclude = FALSE, scalar = FALSE, ... ) ## S3 method for class 'draws_list' subset_draws( x, variable = NULL, iteration = NULL, chain = NULL, draw = NULL, regex = FALSE, unique = TRUE, exclude = FALSE, scalar = FALSE, ... ) ## S3 method for class 'draws_rvars' subset_draws( x, variable = NULL, iteration = NULL, chain = NULL, draw = NULL, regex = FALSE, unique = TRUE, exclude = FALSE, scalar = FALSE, ... ) ## S3 method for class 'rvar' subset_draws(x, variable = NULL, ...) ## S3 method for class 'draws' subset(x, ...)
x 
(draws) A 
... 
Arguments passed to individual methods (if applicable). 
variable 
(character vector) The variables to select. All elements of nonscalar variables can be selected at once. 
iteration 
(integer vector) The iteration indices to select. 
chain 
(integer vector) The chain indices to select. 
draw 
(integer vector) The draw indices to be
select. Subsetting draw indices will lead to an automatic merging
of chains via 
regex 
(logical) Should 
unique 
(logical) Should duplicated selection of chains,
iterations, or draws be allowed? If 
exclude 
(logical) Should the selected subset be excluded?
If 
scalar 
(logical) Should only scalar variables be selected?
If 
To ensure that multiple consecutive subsetting operations work correctly,
subset()
repairs the draws
object before and after
subsetting.
A draws
object of the same class as x
.
x < example_draws() subset_draws(x, variable = c("mu", "tau")) subset_draws(x, chain = 2) subset_draws(x, iteration = 5:10, chain = 3:4) # extract the first chain twice 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)
x < example_draws() subset_draws(x, variable = c("mu", "tau")) subset_draws(x, chain = 2) subset_draws(x, iteration = 5:10, chain = 3:4) # extract the first chain twice 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)
draws
objectsThin draws
objects to reduce their size and autocorrelation in
the chains.
thin_draws(x, thin = NULL, ...) ## S3 method for class 'draws' thin_draws(x, thin = NULL, ...) ## S3 method for class 'rvar' thin_draws(x, thin = NULL, ...)
thin_draws(x, thin = NULL, ...) ## S3 method for class 'draws' thin_draws(x, thin = NULL, ...) ## S3 method for class 'rvar' thin_draws(x, thin = NULL, ...)
x 
(draws) A 
thin 
(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

... 
Arguments passed to individual methods (if applicable). 
A draws
object of the same class as x
.
Teemu Säilynoja, PaulChristian Bürkner, and Aki Vehtari (2022). Graphical test for discrete uniformity and its applications in goodnessoffit evaluation and multiple sample comparison. Statistics and Computing. 32, 32. doi:10.1007/s11222022100906
x < example_draws() niterations(x) x < thin_draws(x, thin = 5) niterations(x)
x < example_draws() niterations(x) x < thin_draws(x, thin = 5) niterations(x)
draws
objectsGet variable names from draws
objects.
variables(x, ...) ## S3 method for class 'draws_matrix' variables(x, reserved = FALSE, with_indices = TRUE, ...) ## S3 method for class 'draws_array' variables(x, reserved = FALSE, with_indices = TRUE, ...) ## S3 method for class 'draws_df' variables(x, reserved = FALSE, with_indices = TRUE, ...) ## S3 method for class 'draws_list' variables(x, reserved = FALSE, with_indices = TRUE, ...) ## S3 method for class 'draws_rvars' variables(x, reserved = FALSE, with_indices = FALSE, ...) nvariables(x, ...)
variables(x, ...) ## S3 method for class 'draws_matrix' variables(x, reserved = FALSE, with_indices = TRUE, ...) ## S3 method for class 'draws_array' variables(x, reserved = FALSE, with_indices = TRUE, ...) ## S3 method for class 'draws_df' variables(x, reserved = FALSE, with_indices = TRUE, ...) ## S3 method for class 'draws_list' variables(x, reserved = FALSE, with_indices = TRUE, ...) ## S3 method for class 'draws_rvars' variables(x, reserved = FALSE, with_indices = FALSE, ...) nvariables(x, ...)
x 
(draws) A 
... 
Arguments passed to individual methods (if applicable). 
reserved 
(logical) Should reserved variables be included in the
output? Defaults to 
with_indices 
(logical) Should indices be included in variable
names? For example, if the object includes variables named 
variables()
returns a vector of all variable names, and nvariables()
returns the number of variables.
For variables()
, a character vector.
For nvariables()
, a scalar integer.
variables<
, rename_variables
, drawsindex
x < example_draws() variables(x) nvariables(x) variables(x) < letters[1:nvariables(x)]
x < example_draws() variables(x) nvariables(x) variables(x) < letters[1:nvariables(x)]
draws
objectsSet variable names for all variables in a draws
object. The
set_variables()
form is useful when using pipe operators.
variables(x, ...) < value ## S3 replacement method for class 'draws_matrix' variables(x, with_indices = TRUE, ...) < value ## S3 replacement method for class 'draws_array' variables(x, with_indices = TRUE, ...) < value ## S3 replacement method for class 'draws_df' variables(x, with_indices = TRUE, ...) < value ## S3 replacement method for class 'draws_list' variables(x, with_indices = TRUE, ...) < value ## S3 replacement method for class 'draws_rvars' variables(x, with_indices = FALSE, ...) < value set_variables(x, variables, ...)
variables(x, ...) < value ## S3 replacement method for class 'draws_matrix' variables(x, with_indices = TRUE, ...) < value ## S3 replacement method for class 'draws_array' variables(x, with_indices = TRUE, ...) < value ## S3 replacement method for class 'draws_df' variables(x, with_indices = TRUE, ...) < value ## S3 replacement method for class 'draws_list' variables(x, with_indices = TRUE, ...) < value ## S3 replacement method for class 'draws_rvars' variables(x, with_indices = FALSE, ...) < value set_variables(x, variables, ...)
x 
(draws) A 
... 
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 
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.
Returns a draws
object of the same format as x
, with
variables named as specified.
variables
, rename_variables
, drawsindex
x < example_draws() variables(x) nvariables(x) variables(x) < letters[1:nvariables(x)] # or equivalently... x < set_variables(x, letters[1:nvariables(x)])
x < example_draws() variables(x) nvariables(x) variables(x) < letters[1:nvariables(x)] # or equivalently... x < set_variables(x, letters[1:nvariables(x)])
draws
objectsAdd weights to draws
objects, with one weight per draw, for use in
subsequent weighting operations. For reasons of numerical accuracy, weights
are stored in the form of unnormalized logweights (in a variable called
.log_weight
). See weights.draws()
for details how to extract weights from
draws
objects.
weight_draws(x, weights, ...) ## S3 method for class 'draws_matrix' weight_draws(x, weights, log = FALSE, pareto_smooth = FALSE, ...) ## S3 method for class 'draws_array' weight_draws(x, weights, log = FALSE, pareto_smooth = FALSE, ...) ## S3 method for class 'draws_df' weight_draws(x, weights, log = FALSE, pareto_smooth = FALSE, ...) ## S3 method for class 'draws_list' weight_draws(x, weights, log = FALSE, pareto_smooth = FALSE, ...) ## S3 method for class 'draws_rvars' weight_draws(x, weights, log = FALSE, pareto_smooth = FALSE, ...)
weight_draws(x, weights, ...) ## S3 method for class 'draws_matrix' weight_draws(x, weights, log = FALSE, pareto_smooth = FALSE, ...) ## S3 method for class 'draws_array' weight_draws(x, weights, log = FALSE, pareto_smooth = FALSE, ...) ## S3 method for class 'draws_df' weight_draws(x, weights, log = FALSE, pareto_smooth = FALSE, ...) ## S3 method for class 'draws_list' weight_draws(x, weights, log = FALSE, pareto_smooth = FALSE, ...) ## S3 method for class 'draws_rvars' weight_draws(x, weights, log = FALSE, pareto_smooth = FALSE, ...)
x 
(draws) A 
weights 
(numeric vector) A vector of weights of length 
... 
Arguments passed to individual methods (if applicable). 
log 
(logical) Are the weights passed already on the log scale? The
default is 
pareto_smooth 
(logical) Should the weights be Paretosmoothed?
The default is 
A draws
object of the same class as x
.
weights.draws()
, resample_draws()
x < example_draws() # sample some random weights for illustration wts < rexp(ndraws(x)) head(wts) # add weights x < weight_draws(x, weights = wts) # extract weights head(weights(x)) # defaults to normalized weights head(weights(x, normalize=FALSE)) # recover original weights head(weights(x, log=TRUE)) # get normalized logweights # add weights which are already on the log scale log_wts < log(wts) head(log_wts) x < weight_draws(x, weights = log_wts, log = TRUE) # extract weights head(weights(x)) head(weights(x, log=TRUE, normalize = FALSE)) # recover original log_wts # add weights on log scale and Pareto smooth them x < weight_draws(x, weights = log_wts, log = TRUE, pareto_smooth = TRUE)
x < example_draws() # sample some random weights for illustration wts < rexp(ndraws(x)) head(wts) # add weights x < weight_draws(x, weights = wts) # extract weights head(weights(x)) # defaults to normalized weights head(weights(x, normalize=FALSE)) # recover original weights head(weights(x, log=TRUE)) # get normalized logweights # add weights which are already on the log scale log_wts < log(wts) head(log_wts) x < weight_draws(x, weights = log_wts, log = TRUE) # extract weights head(weights(x)) head(weights(x, log=TRUE, normalize = FALSE)) # recover original log_wts # add weights on log scale and Pareto smooth them x < weight_draws(x, weights = log_wts, log = TRUE, pareto_smooth = TRUE)
Extract weights from draws
objects, with one weight per draw.
See weight_draws
for details how to add weights to draws
objects.
## S3 method for class 'draws' weights(object, log = FALSE, normalize = TRUE, ...)
## S3 method for class 'draws' weights(object, log = FALSE, normalize = TRUE, ...)
object 
(draws) A 
log 
(logical) Should the weights be returned on the log scale?
Defaults to 
normalize 
(logical) Should the weights be normalized to sum to 1 on
the standard scale? Defaults to 
... 
Arguments passed to individual methods (if applicable). 
A vector of weights, with one weight per draw.
x < example_draws() # sample some random weights for illustration wts < rexp(ndraws(x)) head(wts) # add weights x < weight_draws(x, weights = wts) # extract weights head(weights(x)) # defaults to normalized weights head(weights(x, normalize=FALSE)) # recover original weights head(weights(x, log=TRUE)) # get normalized logweights # add weights which are already on the log scale log_wts < log(wts) head(log_wts) x < weight_draws(x, weights = log_wts, log = TRUE) # extract weights head(weights(x)) head(weights(x, log=TRUE, normalize = FALSE)) # recover original log_wts # add weights on log scale and Pareto smooth them x < weight_draws(x, weights = log_wts, log = TRUE, pareto_smooth = TRUE)
x < example_draws() # sample some random weights for illustration wts < rexp(ndraws(x)) head(wts) # add weights x < weight_draws(x, weights = wts) # extract weights head(weights(x)) # defaults to normalized weights head(weights(x, normalize=FALSE)) # recover original weights head(weights(x, log=TRUE)) # get normalized logweights # add weights which are already on the log scale log_wts < log(wts) head(log_wts) x < weight_draws(x, weights = log_wts, log = TRUE) # extract weights head(weights(x)) head(weights(x, log=TRUE, normalize = FALSE)) # recover original log_wts # add weights on log scale and Pareto smooth them x < weight_draws(x, weights = log_wts, log = TRUE, pareto_smooth = TRUE)