The loo
methods for arrays, matrices, and functions compute PSIS-LOO
CV, efficient approximate leave-one-out (LOO) cross-validation for Bayesian
models using Pareto smoothed importance sampling (PSIS). This is an
implementation of the methods described in Vehtari, Gelman, and Gabry (2017a,
2017b).
The loo_i
function enables testing log-likelihood
functions for use with the loo.function
method.
loo(x, ...)# S3 method for array
loo(x, ..., r_eff = NULL, save_psis = FALSE,
cores = getOption("mc.cores", 1))
# S3 method for matrix
loo(x, ..., r_eff = NULL, save_psis = FALSE,
cores = getOption("mc.cores", 1))
# S3 method for function
loo(x, ..., data = NULL, draws = NULL, r_eff = NULL,
save_psis = FALSE, cores = getOption("mc.cores", 1))
loo_i(i, llfun, ..., data = NULL, draws = NULL, r_eff = NULL)
A log-likelihood array, matrix, or function. See the Methods (by class) section below for a detailed description of how to specify the inputs for each method.
Vector of relative effective sample size estimates for the
likelihood (exp(log_lik)
) of each observation. This is related to
the relative efficiency of estimating the normalizing term in
self-normalizing importance sampling. If r_eff
is not provided then
the reported PSIS effective sample sizes and Monte Carlo error estimates
will be over-optimistic. See the relative_eff
helper function
for computing r_eff
.
Should the "psis"
object created internally by
loo
be saved in the returned object? The loo
function calls
psis
internally but by default discards the (potentially
large) "psis"
object after using it to compute the LOO-CV summaries.
Setting save_psis
to TRUE
will add a psis_object
component to the list returned by loo
. Currently this is only needed
if you plan to use the E_loo
function to compute weighted
expectations after running loo
.
The number of cores to use for parallelization. This defaults to
the option mc.cores
which can be set for an entire R session by
options(mc.cores = NUMBER)
. The old option loo.cores
is now
deprecated but will be given precedence over mc.cores
until
loo.cores
is removed in a future release. As of version
2.0.0 the default is now 1 core if mc.cores
is not set, but we
recommend using as many (or close to as many) cores as possible.
For the loo
function method and the loo_i
function, the data, posterior draws, and other arguments to pass to the
log-likelihood function. See the Methods (by class) section below
for details on how to specify these arguments.
For loo_i
, an integer in 1:N
.
For loo_i
, the same as x
for the
loo.function
method. A log-likelihood function as described in the
Methods (by class) section.
The loo
methods return a named list with class
c("psis_loo", "loo")
and components:
estimates
A matrix with two columns (Estimate
, SE
) and four rows
(elpd_loo
, mcse_elpd_loo
, p_loo
, looic
). This
contains point estimates and standard errors of the expected log pointwise
predictive density (elpd_loo
), the Monte Carlo standard error of
elpd_loo
(mcse_elpd_loo
), the effective number of parameters
(p_loo
) and the LOO information criterion looic
(which is
just -2 * elpd_loo
, i.e., converted to deviance scale).
pointwise
A matrix with four columns (and number of rows equal to the number of
observations) containing the pointwise contributions of each of the above
measures (elpd_loo
, mcse_elpd_loo
, p_loo
,
looic
).
diagnostics
A named list containing two vectors:
pareto_k
: Estimates of the shape parameter \(k\) of the
generalized Pareto fit to the importance ratios for each leave-one-out
distribution. See the pareto-k-diagnostic
page for details.
n_eff
: PSIS effective sample size estimates.
psis_object
This component will be NULL
unless the save_psis
argument is
set to TRUE
when calling loo
. In that case psis_object
will be the object of class "psis"
that is created when the
loo
function calls psis
internally to do the PSIS
procedure.
The loo_i
function returns a named list with components
pointwise
and diagnostics
. These components have the same
structure as the pointwise
and diagnostics
components of the
object returned by loo
except they contain results for only a single
observation.
array
: An \(I\) by \(C\) by \(N\) array, where \(I\)
is the number of MCMC iterations per chain, \(C\) is the number of
chains, and \(N\) is the number of data points.
matrix
: An \(S\) by \(N\) matrix, where \(S\) is the size
of the posterior sample (with all chains merged) and \(N\) is the number
of data points.
function
: A function f
that takes arguments data_i
and draws
and
returns a vector containing the log-likelihood for a single observation
i
evaluated at each posterior draw. The function should be written
such that, for each observation i
in 1:N
, evaluating
f(data_i = data[i,, drop=FALSE], draws = draws)
results in a vector
of length S
(size of posterior sample). The log-likelihood function
can also have additional arguments but data_i
and draws
are
required.
If using the function method then the arguments data
and draws
must also be specified in the call to loo
:
data
: A data frame or matrix containing the data (e.g.
observed outcome and predictors) needed to compute the pointwise
log-likelihood. For each observation i
, the i
th row of
data
will be passed to the data_i
argument of the
log-likelihood function.
draws
: An object containing the posterior draws for any
parameters needed to compute the pointwise log-likelihood. Unlike
data
, which is indexed by observation, for each observation the
entire object draws
will be passed to the draws
argument of
the log-likelihood function.
The ...
can be used to pass additional arguments to your
log-likelihood function. These arguments are used like the draws
argument in that they are recycled for each observation.
Package developers can
define loo
methods for fitted models objects. See the example
loo.stanfit
method in the Examples section below for an
example of defining a method that calls loo.array
. The
loo.stanreg
method in rstanarm is an example of defining a
method that calls loo.function
.
The loo
function is an S3 generic and methods are provided
for computing LOO from 3-D pointwise log-likelihood arrays, pointwise
log-likelihood matrices, and log-likelihood functions. The array and matrix
methods are most convenient, but for models fit to very large datasets the
loo.function
method is more memory efficient and may be preferable.
Vehtari, A., Gelman, A., and Gabry, J. (2017a). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing. 27(5), 1413--1432. doi:10.1007/s11222-016-9696-4. (published version, arXiv preprint).
Vehtari, A., Gelman, A., and Gabry, J. (2017b). Pareto smoothed importance sampling. arXiv preprint: http://arxiv.org/abs/1507.02646/
The loo package vignettes for demonstrations.
psis
for the underlying Pareto Smoothed Importance
Sampling (PSIS) procedure used in the LOO-CV approximation.
pareto-k-diagnostic for convenience functions for looking at diagnostics.
compare
for model comparison.
# NOT RUN { ### Array and matrix methods (using example objects included with loo package) # Array method LLarr <- example_loglik_array() rel_n_eff <- relative_eff(exp(LLarr)) loo(LLarr, r_eff = rel_n_eff, cores = 2) # Matrix method LLmat <- example_loglik_matrix() rel_n_eff <- relative_eff(exp(LLmat), chain_id = rep(1:2, each = 500)) loo(LLmat, r_eff = rel_n_eff, cores = 2) # } # NOT RUN { ### Usage with stanfit objects # see ?extract_log_lik log_lik1 <- extract_log_lik(stanfit1, merge_chains = FALSE) rel_n_eff <- relative_eff(exp(log_lik1)) loo(log_lik1, r_eff = rel_n_eff, cores = 2) # } # NOT RUN { ### Using log-likelihood function instead of array or matrix set.seed(124) # Simulate data and draw from posterior N <- 50; K <- 10; S <- 100; a0 <- 3; b0 <- 2 p <- rbeta(1, a0, b0) y <- rbinom(N, size = K, prob = p) a <- a0 + sum(y); b <- b0 + N * K - sum(y) fake_posterior <- as.matrix(rbeta(S, a, b)) dim(fake_posterior) # S x 1 fake_data <- data.frame(y,K) dim(fake_data) # N x 2 llfun <- function(data_i, draws) { # each time called internally within loo the arguments will be equal to: # data_i: ith row of fake_data (fake_data[i,, drop=FALSE]) # draws: entire fake_posterior matrix dbinom(data_i$y, size = data_i$K, prob = draws, log = TRUE) } # Use the loo_i function to check that llfun works on a single observation # before running on all obs. For example, using the 3rd obs in the data: loo_3 <- loo_i(i = 3, llfun = llfun, data = fake_data, draws = fake_posterior) print(loo_3$pointwise[, "elpd_loo"]) # Use loo.function method loo_with_fn <- loo(llfun, draws = fake_posterior, data = fake_data) # If we look at the elpd_loo contribution from the 3rd obs it should be the # same as what we got above with the loo_i function and i=3: print(loo_with_fn$pointwise[3, "elpd_loo"]) print(loo_3$pointwise[, "elpd_loo"]) # Check that the loo.matrix method gives same answer as loo.function method log_lik_matrix <- sapply(1:N, function(i) { llfun(data_i = fake_data[i,, drop=FALSE], draws = fake_posterior) }) loo_with_mat <- loo(log_lik_matrix) all.equal(loo_with_mat$estimates, loo_with_fn$estimates) # should be TRUE! # } # NOT RUN { ### For package developers: defining loo methods # An example of a possible loo method for 'stanfit' objects (rstan package). # A similar method is planned for a future release of rstan (or is already # released, depending on when you are reading this). In order for users # to be able to call loo(stanfit) instead of loo.stanfit(stanfit) the # NAMESPACE needs to be handled appropriately (roxygen2 and devtools packages # are good for that). # loo.stanfit <- function(x, pars = "log_lik", ..., save_psis = FALSE, cores = getOption("mc.cores", 1)) { stopifnot(length(pars) == 1L) LLarray <- loo::extract_log_lik(stanfit = x, parameter_name = pars, merge_chains = FALSE) r_eff <- loo::relative_eff(x = exp(LLarray), cores = cores) loo::loo.array(LLarray, r_eff = r_eff, cores = cores, save_psis = save_psis) } # } # NOT RUN { # }
Run the code above in your browser using DataCamp Workspace