Working with Posteriors

Summary statistics

We can easily customize the summary statistics reported by $summary() and $print().

fit <- cmdstanr::cmdstanr_example("schools", method = "sample")
fit$summary()

By default, all variables are summarized with the following functions:

posterior::default_summary_measures()

To change the variables summarized, use the variables argument:

fit$summary(variables = c("mu", "tau"))

We can also change which functions are used:

fit$summary(variables = c("mu", "tau"), mean, sd)

To summarize all variables with non-default functions, it is necessary to explicitly set the variables argument, either to NULL or the full vector of variable names.

fit$summary(variables = NULL, "mean", "median")

Summary functions can be specified by character string, function, or using a formula (or anything else supported by rlang::as_function()). If these arguments are named, those names will be used in the tibble output. If the summary results are named they will take precedence.

my_sd <- function(x) c(My_SD = sd(x))
fit$summary(
  c("mu", "tau"), 
  MEAN = mean, 
  "median",
  my_sd,
  ~quantile(.x, probs = c(0.1, 0.9)),
  Minimum = function(x) min(x)
)

Arguments to all summary functions can also be specified with .args.

fit$summary(c("mu", "tau"), quantile, .args = list(probs = c(0.025, .05, .95, .975)))

Each summary function is applied separately to each variable and receives a matrix whose rows are saved iterations and whose columns are chains.

fit$summary(variables = NULL, dim, colMeans)

For this reason users may have unexpected results if they use stats::var() directly, as it will return a covariance matrix. An alternative is the distributional::variance() function, which can also be accessed via posterior::variance().

fit$summary(c("mu", "tau"), posterior::variance, ~var(as.vector(.x)))

Summary functions need not return numeric values when used with $summary(). The $print() method requires numeric summary columns because it rounds them to the requested number of digits.

strict_pos <- function(x) if (all(x > 0)) "yes" else "no"
fit$summary(variables = NULL, "Strictly Positive" = strict_pos)
# fit$print(variables = NULL, "Strictly Positive" = strict_pos)

For more information, see posterior::summarise_draws(), which is called by $summary().

Extracting posterior draws/samples

The $draws() method extracts draws in formats provided by the posterior package. The Getting started with CmdStanR vignette introduces the most commonly used formats and how to convert between them.

fit$draws("mu")
fit$draws("theta")
fit$draws(c("mu", "theta[1]"), format = "df")

For MCMC fits, inc_warmup = TRUE includes warmup draws, but only if save_warmup = TRUE was specified when fitting the model.

For more ways to manipulate draws, see the posterior package vignettes and documentation.

Structured draws similar to rstan::extract()

The posterior package provides two useful ways to work with variables while preserving their original dimensions.

posterior::extract_list_of_variable_arrays() returns a named list containing one array per variable. Setting with_chains = FALSE combines the chains, giving the same general structure as the list returned by rstan::extract():

draw_arrays <- posterior::extract_list_of_variable_arrays(
  fit$draws(),
  variables = c("mu", "theta"),
  with_chains = FALSE
)
str(draw_arrays)

The first dimension of each array indexes draws, and any remaining dimensions match the dimensions of the corresponding Stan variable.

Alternatively, the posterior package’s rvar format represents each variable as a multidimensional random variable, with its posterior draws handled behind the scenes:

draws_rvars <- posterior::as_draws_rvars(
  fit$draws(c("mu", "theta"))
)
theta_rvar <- draws_rvars$theta

# Compute the difference for every draw using natural vector indexing
# The posterior draws are handled automatically
theta_difference <- theta_rvar[1] - theta_rvar[2]
theta_difference

hist(
  posterior::draws_of(theta_difference),
  main = "Difference between theta[1] and theta[2]",
  xlab = "theta[1] - theta[2]"
)

# Direct access to the underlying draws is also available
theta_array <- posterior::draws_of(theta_rvar)
dim(theta_array)

The object theta_rvar behaves like the vector declared in the Stan program. theta_array provides direct access to its underlying draws, with the first dimension indexing draws. See the rvar vignette for details.