,
Teemu Säilynoja [aut],
Luna Fazio [ctb]| Title: | Simulation Based Calibration for rstan/cmdstanr models |
|---|---|
| Description: | SBC helps perform Simulation Based Calibration on Bayesian models. SBC lets you check for bugs in your model code and/or algorithm that fits the model. SBC focuses on models built with 'Stan' <https://mc-stan.org>, but can support other modelling languages as well. |
| Authors: | Shinyoung Kim [aut],
Angie.H Moon [aut, cre],
Martin Modrák [aut] ,
Teemu Säilynoja [aut],
Luna Fazio [ctb] |
| Maintainer: | Angie.H Moon <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.3.0.9000 |
| Built: | 2024-06-22 09:46:06 UTC |
| Source: | https://github.com/hyunjimoon/SBC |
SBC_datasets object.Subset an SBC_datasets object.
## S3 method for class 'SBC_datasets' x[indices]## S3 method for class 'SBC_datasets' x[indices]
Subset the results.
## S3 method for class 'SBC_results' x[indices]## S3 method for class 'SBC_results' x[indices]
indices |
integer indices or a binary vector of the same length as the number fits, selecting which fits to keep. |
Combine multiple datasets together.
bind_datasets(...)bind_datasets(...)
... |
datasets to bind |
Combine two lists of derived quantities
bind_derived_quantities(dq1, dq2)bind_derived_quantities(dq1, dq2)
Combine two sets globals for use in derived quantities or backend
bind_globals(globals1, globals2)bind_globals(globals1, globals2)
compute_SBC(), derived_quantities()
Primarily useful for iteratively adding more simulations to your SBC check.
bind_results(...)bind_results(...)
... |
objects of type |
An example usage can be found in the small_model_workflow vignette.
Calculate prior standard deviation of a dataset
calculate_prior_sd(datasets)calculate_prior_sd(datasets)
a named vector of prior SDs
When there are ties (e.g. for discrete variables), the rank is currently drawn stochastically among the ties. NA is assumed to be potentially equal to any value (Issue #78)
calculate_ranks_draws_matrix(variables, dm, params = NULL)calculate_ranks_draws_matrix(variables, dm, params = NULL)
variables |
a vector of values to check |
dm |
draws_matrix of the fit (assumed to be already thinned if that was necessary) |
params |
DEPRECATED. Use |
Check diagnostics and issue warnings when those fail.
check_all_SBC_diagnostics(x) ## Default S3 method: check_all_SBC_diagnostics(x) ## S3 method for class 'SBC_results' check_all_SBC_diagnostics(x)check_all_SBC_diagnostics(x) ## Default S3 method: check_all_SBC_diagnostics(x) ## S3 method for class 'SBC_results' check_all_SBC_diagnostics(x)
TRUE if all checks are OK, FALSE otherwise.
Computes the cumulative Jensen-Shannon distance between two samples.
cjs_dist( x, y, x_weights = rep(1/length(x), length(x)), y_weights = rep(1/length(y), length(y)), ... )cjs_dist( x, y, x_weights = rep(1/length(x), length(x)), y_weights = rep(1/length(y), length(y)), ... )
x |
numeric vector of draws from first distribution |
y |
numeric vector of draws from second distribution |
x_weights |
numeric vector of weights of first distribution |
y_weights |
numeric vector of weights of second distribution |
... |
unused |
The Cumulative Jensen-Shannon distance is a symmetric metric based on the cumulative Jensen-Shannon divergence. The divergence CJS(P || Q) between two cumulative distribution functions P and Q is defined as:
The symmetric metric is defined as:
This has an upper bound of \sqrt \sum (P(x) + Q(x))
distance value based on CJS computation.
Nguyen H-V., Vreeken J. (2015). Non-parametric
Jensen-Shannon Divergence. In: Appice A., Rodrigues P., Santos
Costa V., Gama J., Jorge A., Soares C. (eds) Machine Learning
and Knowledge Discovery in Databases. ECML PKDD 2015. Lecture
Notes in Computer Science, vol 9285. Springer, Cham.
doi:10.1007/978-3-319-23525-7_11
combine_all_variables will merge all variables in a single plot, while
combine_array_elements will merge all elements of any array into a single
panel of the plot
combine_all_variables(x) combine_array_elements(x)combine_all_variables(x) combine_array_elements(x)
x |
parameter names |
Combine two named lists and overwrite elements with the same name using the value from args2
combine_args(args1, args2)combine_args(args1, args2)
Compute derived quantities based on given data and posterior draws.
compute_dquants(draws, generated, dquants, gen_quants = NULL)compute_dquants(draws, generated, dquants, gen_quants = NULL)
gen_quants |
Deprecated, use |
Performs the main SBC routine given datasets and a backend.
compute_SBC( datasets, backend, cores_per_fit = default_cores_per_fit(length(datasets)), keep_fits = TRUE, thin_ranks = SBC_backend_default_thin_ranks(backend), ensure_num_ranks_divisor = 2, chunk_size = default_chunk_size(length(datasets)), dquants = NULL, cache_mode = "none", cache_location = NULL, globals = list(), gen_quants = NULL )compute_SBC( datasets, backend, cores_per_fit = default_cores_per_fit(length(datasets)), keep_fits = TRUE, thin_ranks = SBC_backend_default_thin_ranks(backend), ensure_num_ranks_divisor = 2, chunk_size = default_chunk_size(length(datasets)), dquants = NULL, cache_mode = "none", cache_location = NULL, globals = list(), gen_quants = NULL )
datasets |
an object of class |
backend |
the model + sampling algorithm. The built-in backends can be constructed
using |
cores_per_fit |
how many cores should the backend be allowed to use for a single fit?
Defaults to the maximum number that does not produce more parallel chains
than you have cores. See |
keep_fits |
boolean, when |
thin_ranks |
how much thinning should be applied to posterior draws before computing ranks for SBC. Should be large enough to avoid any noticeable autocorrelation of the thinned draws See details below. |
ensure_num_ranks_divisor |
Potentially drop some posterior samples to ensure that this number divides the total number of SBC ranks (see Details). |
chunk_size |
How many simulations within the |
dquants |
Derived quantities to include in SBC. Use |
cache_mode |
Type of caching of results, currently the only supported modes are
|
cache_location |
The filesystem location of cache. For |
globals |
A list of names of objects that are defined
in the global environment and need to present for the backend to work (
if they are not already available in package).
It is added to the |
gen_quants |
Deprecated, use dquants instead |
An object of class SBC_results().
Parallel processing is supported via the future package, for most uses, it is most sensible
to just call plan(multisession) once in your R session and all
cores your computer will be used. For more details refer to the documentation
of the future package.
When using backends based on MCMC, there are two possible moments when
draws may need to be thinned. They can be thinned directly within the backend
and they may be thinned only to compute the ranks for SBC as specified by the
thin_ranks argument. The main reason those are separate is that computing the
ranks requires no or negligible autocorrelation while some autocorrelation
may be easily tolerated for summarising the fit results or assessing convergence.
In fact, thinning too aggressively in the backend may lead to overly noisy
estimates of posterior means, quantiles and the posterior::rhat() and
posterior::ess_tail() diagnostics. So for well-adapted Hamiltonian Monte-Carlo
chains (e.g. Stan-based backends), we recommend no thinning in the backend and
even value of thin_ranks between 6 and 10 is usually sufficient to remove
the residual autocorrelation. For a backend based on Metropolis-Hastings,
it might be sensible to thin quite aggressively already in the backend and
then have some additional thinning via thin_ranks.
Backends that don't require thining should implement SBC_backend_iid_draws()
or SBC_backend_default_thin_ranks() to avoid thinning by default.
Some of the visualizations and post processing steps
we use in the SBC package (e.g. plot_rank_hist(), empirical_coverage())
work best if the total number of possible SBC ranks is a "nice" number
(lots of divisors).
However, the number of ranks is one plus the number of posterior samples
after thinning - therefore as long as the number of samples is a "nice"
number, the number of ranks usually will not be. To remedy this, you can
specify ensure_num_ranks_divisor - the method will drop at most
ensure_num_ranks_divisor - 1 samples to make the number of ranks divisible
by ensure_num_ranks_divisor. The default 2 prevents the most annoying
pathologies while discarding at most a single sample.
Maybe not export in the end? Useful for debugging
data_for_ecdf_plots(x, ..., prob = 0.95, gamma = NULL, K = NULL)data_for_ecdf_plots(x, ..., prob = 0.95, gamma = NULL, K = NULL)
By default will make every worker process a single chunk.
You can set the options(SBC.min_chunk_size = value) to enforce a minimum
chunk size globally (chunk size can still be larger if you have substantially more
fits to run than workers.
default_chunk_size(n_fits, n_workers = future::nbrOfWorkers())default_chunk_size(n_fits, n_workers = future::nbrOfWorkers())
When parallel processing is disabled, this just returns the number of available
cores. Otherwise, it chooses the largest integer that keeps
cores_per_fit * (n_fits / chunk_size) <= total_cores, i.e. it avoids
running so many chains in parallel that there will be more chains than cores.
default_cores_per_fit( n_fits, total_cores = future::availableCores(), chunk_size = default_chunk_size(n_fits) )default_cores_per_fit( n_fits, total_cores = future::availableCores(), chunk_size = default_chunk_size(n_fits) )
When the expression contains non-library functions/objects, and parallel processing
is enabled, those must be
named in the .globals parameter (hopefully we'll be able to detect those
automatically in the future). Note that recompute_SBC_statistics() currently
does not use parallel processing, so .globals don't need to be set.
derived_quantities(..., .globals = list())derived_quantities(..., .globals = list())
... |
named expressions representing the quantitites |
.globals |
A list of names of objects that are defined
in the global environment and need to present for the gen. quants. to evaluate.
It is added to the |
# Derived quantity computing the total log likelihood of a normal distribution # with known sd = 1 normal_lpdf <- function(y, mu, sigma) { sum(dnorm(y, mean = mu, sd = sigma, log = TRUE)) } # Note the use of .globals to make the normal_lpdf function available # within the expression log_lik_dq <- derived_quantities(log_lik = normal_lpdf(y, mu, 1), .globals = "normal_lpdf" )# Derived quantity computing the total log likelihood of a normal distribution # with known sd = 1 normal_lpdf <- function(y, mu, sigma) { sum(dnorm(y, mean = mu, sd = sigma, log = TRUE)) } # Note the use of .globals to make the normal_lpdf function available # within the expression log_lik_dq <- derived_quantities(log_lik = normal_lpdf(y, mu, 1), .globals = "normal_lpdf" )
Uses ranks to compute coverage and surrounding uncertainty of posterior credible intervals.
The uncertainty is only approximate (treating coverage for each interval as a set of independent
Bernoulli trials, while in fact they are not independent), so for making claims on presence/
absence of detectable discrepancies we strongly recommend using plot_ecdf() or plot_ecdf_diff().
The uncertainty about the coverage can however be useful for guiding decisions on whether
more SBC steps should be performed (i.e. whether we can rule out that the coverage of
the given backend differs too much for our purposes from the optimal value).
empirical_coverage(stats, width, prob = 0.95, interval_type = "central")empirical_coverage(stats, width, prob = 0.95, interval_type = "central")
stats |
a data.frame of rank statistics (e.g. as returned in the |
width |
a vector of values between 0 and 1 representing widths of credible intervals for which we compute coverage. |
prob |
determines width of the uncertainty interval around the observed coverage |
inteval_type |
|
Note that while coverage of central posterior intervals (with the default type = "central")
is often of the biggest practical interest, perfect calibration of central intervals
still leaves space for substantial problems with the model (e.g. if the posterior 25% - 50% intervals
contain 50% of the true values and the posterior 50% - 75% interval never contains the true value,
the central 50% interval still has the ideal 50% coverage), so investigating central
intervals should always be accompanied by checks with plot_ecdf() or plot_ecdf_diff()
or by using type = "leftmost", because if all leftmost credible intervals are well calibrated,
then all intervals are well calibrated.
A data.frame with columns variable, width (width of the interval as given
in the width parameter), width_represented the closest width that can be represented by
the ranks in the input (any discrepancy needs to be judged against this rather than width),
estimate - observed coverage for the interval, ci_low, ci_high the uncertainty
interval around estimate (width of the interval is given by the prob argument).
Generate datasets.
generate_datasets(generator, n_sims, n_datasets = NULL)generate_datasets(generator, n_sims, n_datasets = NULL)
generator |
a generator object - build e.g. via |
n_sims |
the number of simulated datasets to use |
n_datasets |
DEPRECATED, use |
object of class SBC_datasets
TODO: seed
SBC_results or other objects.Get diagnostic messages for SBC_results or other objects.
get_diagnostic_messages(x)get_diagnostic_messages(x)
An object of class SBC_diagnostic_messages, inheriting a data.frame.
plot_rank_hist().Guess the number of bins for plot_rank_hist().
guess_rank_hist_bins(max_rank, N)guess_rank_hist_bins(max_rank, N)
max_rank |
the maximum rank observed |
N |
the number of ranks observed |
Max difference between binned samples with the same length
max_diff(x, y)max_diff(x, y)
x |
numeric vector of density from first distribution |
y |
numeric vector of density from second distribution |
... |
unused |
distance value based on max difference
The rationale for this plot and its interpretation is explained in Mike Betancourt's Towards A Principled Bayesian Workflow.
plot_contraction( x, prior_sd, variables = NULL, scale = "sd", alpha = 0.8, parameters = NULL )plot_contraction( x, prior_sd, variables = NULL, scale = "sd", alpha = 0.8, parameters = NULL )
x |
object containing results (a data.frame or |
prior_sd |
a named vector of prior standard deviations for your variables.
Either pass in analytically obtained values or use |
variables |
variables to show in the plot or |
scale |
which scale of variability you want to see - either |
alpha |
the alpha for the points |
parameters |
DEPRECATED, use |
a ggplot2 plot object
plot_coverage will plot the observed coverage,
while plot_coverage_diff will show the difference between observed
and expected coverage.
Please refer to empirical_coverage() for details on computation
and limitations of this plot as well as details on the arguments.
See vignette("rank_visualizations") for
more details.
plot_coverage( x, variables = NULL, prob = 0.95, interval_type = "central", parameters = NULL, max_points = NULL ) plot_coverage_diff( x, variables = NULL, prob = 0.95, interval_type = "central", parameters = NULL, max_points = NULL )plot_coverage( x, variables = NULL, prob = 0.95, interval_type = "central", parameters = NULL, max_points = NULL ) plot_coverage_diff( x, variables = NULL, prob = 0.95, interval_type = "central", parameters = NULL, max_points = NULL )
x |
object containing results (a data.frame or |
variables |
variables to show in the plot or |
prob |
the with of the uncertainty interval to be shown |
parameters |
DEPRECATED. Use |
max_points |
maximum number of points where to evaluate the coverage.
If set to |
a ggplot2 plot object
empirical_coverage
See vignette("rank_visualizations") for
more details.
See the methods for data_for_ecdf_plots() for available data formats.
plot_ecdf( x, variables = NULL, K = NULL, gamma = NULL, prob = 0.95, size = 1, alpha = 0.33, combine_variables = NULL, ecdf_alpha = NULL, ..., parameters = NULL ) plot_ecdf_diff( x, variables = NULL, K = NULL, gamma = NULL, prob = 0.95, size = 1, alpha = 0.33, combine_variables = NULL, ecdf_alpha = NULL, ..., parameters = NULL )plot_ecdf( x, variables = NULL, K = NULL, gamma = NULL, prob = 0.95, size = 1, alpha = 0.33, combine_variables = NULL, ecdf_alpha = NULL, ..., parameters = NULL ) plot_ecdf_diff( x, variables = NULL, K = NULL, gamma = NULL, prob = 0.95, size = 1, alpha = 0.33, combine_variables = NULL, ecdf_alpha = NULL, ..., parameters = NULL )
x |
object supporting the |
variables |
optional subset of variables to show in the plot |
K |
number of uniformly spaced evaluation points for the ECDF or ECDFs. Affects the granularity of the plot and can significantly speed up the computation of the simultaneous confidence bands. Defaults to the smaller of number of ranks per variable and the maximum rank. |
gamma |
TODO |
prob |
the width of the plotted confidence interval for the ECDF. |
size |
size passed to |
alpha |
alpha level of the confidence band |
combine_variables |
optionally specify a named list where each entry is a character vectors which specifies a group of variables that will be displayed in the same panel. Panel title will be the name of the list element. A function that takes a character vector as an input and produces a list can also be specified (see combine-functions). |
ecdf_alpha |
the alpha level of the empirical CDF. Can be either a single number or
a function taking the number of variables that were combined (when |
... |
additional arguments passed to |
parameters |
DEPRECATED, use |
arxiv::1903.08008 by A. Vehtari et al.
The expected uniform distribution and an approximate confidence interval
is also shown. The confidence interval cannot be taken too seriously
as it is derived assuming the bins are independent (which they are not).
The plot_ecdf() and plot_ecdf_diff() plots provide better confidence interval
but are somewhat less interpretable. See vignette("rank_visualizations") for
more details.
plot_rank_hist( x, variables = NULL, bins = NULL, prob = 0.95, ..., parameters = NULL )plot_rank_hist( x, variables = NULL, bins = NULL, prob = 0.95, ..., parameters = NULL )
x |
Object supporting the plotting method. |
variables |
Names of variables to show |
bins |
number of bins to be used in the histogram, if left unspecified,
it is determined by |
prob |
The width of the approximate confidence interval shown. |
By default the support is for SBC_results objects and data frames in the same
format as the $stats element of SBC_results.
Plot the simulated "true" values versus posterior estimates
plot_sim_estimated( x, variables = NULL, estimate = "mean", uncertainty = c("q5", "q95"), alpha = NULL, parameters = NULL )plot_sim_estimated( x, variables = NULL, estimate = "mean", uncertainty = c("q5", "q95"), alpha = NULL, parameters = NULL )
x |
object containing results (a data.frame or |
variables |
variables to show in the plot or |
estimate |
which estimate to use for the central tendency,
must correspond a field already computed in the results (most likely |
uncertainty |
which estimates to use for uncertainty (a character vector of length 2)
must correspond a field already computed in the results. Pass |
alpha |
the alpha for the points and uncertainty intervals |
parameters |
DEPRECATED, use |
a ggplot2 plot object
Distance between binned draws (rank for SBC) and discrete uniform
rank2unif(results, par, bins = 20)rank2unif(results, par, bins = 20)
par |
names of parameter to plot |
bins |
number of bins to use for summary |
ranks |
array of dimension (n_iter, n_pars) where n_iter=number of posterior draw iterations, n_pars the number of parameters of interest |
thin |
integer in which thinning was applied |
Useful for example to recompute SBC ranks with a different choice of thin_ranks
or added derived quantities.
recompute_SBC_statistics( old_results, datasets, backend, thin_ranks = SBC_backend_default_thin_ranks(backend), ensure_num_ranks_divisor = 2, dquants = NULL, gen_quants = NULL )recompute_SBC_statistics( old_results, datasets, backend, thin_ranks = SBC_backend_default_thin_ranks(backend), ensure_num_ranks_divisor = 2, dquants = NULL, gen_quants = NULL )
datasets |
an object of class |
backend |
backend used to fit the results. Used to pull various defaults and other setting influencing the computation of statistics. |
thin_ranks |
how much thinning should be applied to posterior draws before computing ranks for SBC. Should be large enough to avoid any noticeable autocorrelation of the thinned draws See details below. |
ensure_num_ranks_divisor |
Potentially drop some posterior samples to ensure that this number divides the total number of SBC ranks (see Details). |
dquants |
Derived quantities to include in SBC. Use |
gen_quants |
Deprecated, use dquants instead |
An S3 object of class SBC_results with updated $stats and $default_diagnostics fields.
brms package.Build a backend based on the brms package.
SBC_backend_brms( ..., template_data, out_stan_file = NULL, template_dataset = NULL )SBC_backend_brms( ..., template_data, out_stan_file = NULL, template_dataset = NULL )
... |
arguments passed to |
template_data |
a representative value for the |
out_stan_file |
A filename for the generated Stan code. Useful for debugging and for avoiding unnecessary recompilations. |
template_dataset |
DEPRECATED. Use |
SBC_generator_brms
object.Build a brms backend, reusing the compiled model from a previously created SBC_generator_brms
object.
SBC_backend_brms_from_generator(generator, ...)SBC_backend_brms_from_generator(generator, ...)
cmdstanr.Backend based on sampling via cmdstanr.
SBC_backend_cmdstan_sample(model, ..., init_factory = NULL)SBC_backend_cmdstan_sample(model, ..., init_factory = NULL)
model |
an object of class |
... |
other arguments passed to the |
init_factory |
an optional function that takes in a dataset and returns a value that
can be passed to the |
cmdstanr.Backend based on variational approximation via cmdstanr.
SBC_backend_cmdstan_variational(model, ..., n_retries_init = 1)SBC_backend_cmdstan_variational(model, ..., n_retries_init = 1)
model |
an object of class |
... |
other arguments passed to the |
n_retries_init |
number of times to retry the variational fit if the algorithm has trouble initializing (e.g. too many dropped evaluations (see https://discourse.mc-stan.org/t/advi-too-many-dropped-evaluations-even-for-well-behaved-models/24338), or "cannot compute ELBO using the initial variational distribution") |
The default implementation plays it relatively safe and returns 10, unless
SBC_backend_iid_draws() returns TRUE in which case it returns 1.
SBC_backend_default_thin_ranks(backend) ## Default S3 method: SBC_backend_default_thin_ranks(backend)SBC_backend_default_thin_ranks(backend) ## Default S3 method: SBC_backend_default_thin_ranks(backend)
S3 generic that allows backends to override how a hash is computed. By default rlang::hash()
is used.
SBC_backend_hash_for_cache(backend)SBC_backend_hash_for_cache(backend)
Most backends (e.g. those based on variatns of MCMC) don't produce
independent draws and thus diagnostics like Rhat and ESS are important
and draws may need thinning. Backends that already produce independent
draws (e.g. ADVI/optimizing) can implement this method to return TRUE
to signal this is the case. If this method returns TRUE, ESS and Rhat will
always attain their best possible values and SBC_backend_default_thin_ranks()
will return 1.
The default implementation returns FALSE.
SBC_backend_iid_draws(backend) ## Default S3 method: SBC_backend_iid_draws(backend)SBC_backend_iid_draws(backend) ## Default S3 method: SBC_backend_iid_draws(backend)
backend |
to check |
Mock backend is useful for testing the package.
It will ignore all data passed to it and always
provide result as the draws generated by the backend.
SBC_backend_mock( result = posterior::draws_matrix(a = rnorm(100)), output = NULL, message = NULL, warning = NULL, error = NULL )SBC_backend_mock( result = posterior::draws_matrix(a = rnorm(100)), output = NULL, message = NULL, warning = NULL, error = NULL )
result |
a |
rjags
Create a JAGS backend using rjags
SBC_backend_rjags( file, n.iter, n.burnin, variable.names, thin = 1, na.rm = TRUE, ... )SBC_backend_rjags( file, n.iter, n.burnin, variable.names, thin = 1, na.rm = TRUE, ... )
file |
model file or connection to model code (passed to |
n.iter |
number of iterations for sampling (passed to [rjags::coda.samples()) |
n.burnin |
number of iterations used for burnin |
variable.names |
names of variables to monitor (passed to |
thin |
thinning (passed to |
na.rm |
whether to omit variables containing NA (passed to |
... |
additional optional arguments passed to
|
optimizing method from rstan.SBC backend using the optimizing method from rstan.
SBC_backend_rstan_optimizing(model, ..., n_retries_hessian = 1)SBC_backend_rstan_optimizing(model, ..., n_retries_hessian = 1)
model |
a |
... |
other arguments passed to |
n_retries_hessian |
the number of times the backend is allow to retry optimization (with different seeed) to produce a usable Hessian that can produce draws. In some cases, the Hessian may be numerically unstable and not be positive definite. |
sampling method from rstan.SBC backend using the sampling method from rstan.
SBC_backend_rstan_sample(model, ...)SBC_backend_rstan_sample(model, ...)
model |
a |
... |
other arguments passed to |
SBC_datasets object.In most cases, you may want to use generate_datasets to build the object, but
for full control, you can also create datasets directly via this function.
SBC_datasets(variables, generated, parameters = NULL)SBC_datasets(variables, generated, parameters = NULL)
variables |
draws of "true" values of unobserved parameters or other derived variables.
An object of class |
generated |
a list of objects that can be passed as data to the backend you plan to use.
(e.g. list of values for Stan-based backends, a data frame for |
parameters |
DEPRECATED. Use variables instead. |
Note that this will involve compiling a Stan model and may take a while.
SBC_example_backend( example = c("normal_sd", "normal_bad"), interface = c("rstan", "cmdstanr", "rjags") )SBC_example_backend( example = c("normal_sd", "normal_bad"), interface = c("rstan", "cmdstanr", "rjags") )
example |
name of the example model. |
interface |
name of the interface to be used to fit the model |
Construct a generator used in the examples.
SBC_example_generator(example = c("normal"), N = 100)SBC_example_generator(example = c("normal"), N = 100)
example |
name of example |
N |
size of the dataset the generator should simulate |
an object that can be passed to generate_datasets()
Except for example = "visualizations", all examples will actually
compile and fit Stan models and thus may take a while to complete.
SBC_example_results( example = c("normal_ok", "normal_bad", "visualizations"), interface = c("rstan", "cmdstanr", "rjags"), N = 100, n_sims = 50 )SBC_example_results( example = c("normal_ok", "normal_bad", "visualizations"), interface = c("rstan", "cmdstanr", "rjags"), N = 100, n_sims = 50 )
example |
|
interface |
name of the interface to be used for the backend |
N |
number of datapoints to simulate from the generator for each simulation |
n_sims |
number of simulations to perform |
Needs to be implemented by all backends.
All implementations have to return an object for which you can safely
call SBC_fit_to_draws_matrix() and get some draws.
If that's not possible an error should be raised.
SBC_fit(backend, generated, cores)SBC_fit(backend, generated, cores)
The diagnostics object has to be a data.frame but may
inherit additional classes - in particular it may be useful
for the returning object to implement get_diagnostic_messages().
SBC_fit_to_diagnostics(fit, fit_output, fit_messages, fit_warnings)SBC_fit_to_diagnostics(fit, fit_output, fit_messages, fit_warnings)
fit |
The fit returned by |
fit_output |
a character string capturing what the backend wrote to stdout |
fit_messages |
a character vector of messages the backend raised |
fit_warnings |
a character vector of warnings the backend raised |
an single row data.frame that includes diagnostics or NULL, if no diagnostics available.
draws_matrix object.Needs to be implemented for all types of objects the backend can
return from SBC_fit(). Default implementation just calls,
posterior::as_draws_matrix(), so if the fit object already supports
this, it will work out of the box.
SBC_fit_to_draws_matrix(fit) ## Default S3 method: SBC_fit_to_draws_matrix(fit)SBC_fit_to_draws_matrix(fit) ## Default S3 method: SBC_fit_to_draws_matrix(fit)
Brms generator uses a brms model with sample_prior = "only" to generate
new datasets.
SBC_generator_brms( formula, data, ..., generate_lp = TRUE, generate_lp_chunksize = 5000/nrow(data), out_stan_file = NULL )SBC_generator_brms( formula, data, ..., generate_lp = TRUE, generate_lp_chunksize = 5000/nrow(data), out_stan_file = NULL )
formula, data, ...
|
arguments passed to |
generate_lp |
whether to compute the overall log-likelihood of the model as an additional variable. This can be somewhat computationally expensive, but improves sensitivity of the SBC process. |
generate_lp_chunksize |
Will determine whether
the chunk size used with |
out_stan_file |
A filename for the generated Stan code. Useful for debugging and for avoiding unnecessary recompilations. |
This creates a very thin wrapper around the function and can store additional arguments, but does not do anything more..
SBC_generator_custom(f, ...)SBC_generator_custom(f, ...)
f |
function accepting at least an |
... |
Additional arguments passed to |
Running:
gen <- SBC_generator_custom(f, <<some other args>>) datasets <- generate_datasets(gen, n_sims = my_n_sims)
is equivalent to just running
datasets <- f(<<some other args>>, n_sims = my_n_sims)
So whenever you control the code calling generate_datasets,
it usually makes more sense to just create an SBC_datasets
object directly and avoid using SBC_generator_custom and generate_datasets at all.
SBC_generator_custom can however be useful, when a code you
do not control calls generate_datasets for you and the
built-in generators do not provide you with enough flexibility.
Generate datasets via a function that creates a single dataset.
SBC_generator_function(f, ...)SBC_generator_function(f, ...)
f |
function returning a list with elements |
... |
Additional arguments passed to |
Print the Stan code of a model used in the examples.
SBC_print_example_model( example = c("normal_sd", "normal_bad"), interface = c("rstan", "cmdstanr", "rjags") )SBC_print_example_model( example = c("normal_sd", "normal_bad"), interface = c("rstan", "cmdstanr", "rjags") )
example |
name of the example model. |
SBC_results objectThis will build and validate an SBC_results object from its constituents.
SBC_results( stats, fits, backend_diagnostics, default_diagnostics, outputs, messages, warnings, errors )SBC_results( stats, fits, backend_diagnostics, default_diagnostics, outputs, messages, warnings, errors )
The SBC_results contains the following fields:
$stats statistics for all variables and fits (one row per variable-fit combination)
$fits the raw fits (unless keep_fits = FALSE) or NULL if the fit failed
$errors error messages that caused fit failures
$outputs, $messages, $warnings the outputs/messages/warnings written by fits
$default_diagnostics a data frame of default convergence/correctness diagnostics (one row per fit)
$backend_diagnostics a data frame of backend-specific diagnostics (one row per fit)
Potentially useful for doing some advanced stuff, but should not
be used in regular workflow. Use recompute_SBC_statistics() to update
an [SBC_results] objects with different thin_ranks or other settings.
SBC_statistics_from_single_fit( fit, variables, generated, thin_ranks, ensure_num_ranks_divisor, dquants, backend, gen_quants = NULL )SBC_statistics_from_single_fit( fit, variables, generated, thin_ranks, ensure_num_ranks_divisor, dquants, backend, gen_quants = NULL )
thin_ranks |
how much thinning should be applied to posterior draws before computing ranks for SBC. Should be large enough to avoid any noticeable autocorrelation of the thinned draws See details below. |
ensure_num_ranks_divisor |
Potentially drop some posterior samples to ensure that this number divides the total number of SBC ranks (see Details). |
dquants |
Derived quantities to include in SBC. Use |
backend |
the model + sampling algorithm. The built-in backends can be constructed
using |
gen_quants |
Deprecated, use dquants instead |
Summarize relational property of overall prior and posterior draws
set2set(priors, posteriors, par, bins = 20)set2set(priors, posteriors, par, bins = 20)
priors |
A posterior::draws_rvars of dimension(n_iterations=1, n_chains=n_sbc_iterations, n_variables=n_variables) which stores prior draws |
posteriors |
A posterior::draws_Rvars of dimension(n_iterations=n_posterior_draws, n_chains=n_sbc_iterations, n_variables=n_variables), which stores fitted posterior draws |
par |
names of parameter to summarize |
bins |
number of bins for prior and post density |
Validate a definition of derived quantities evaluated in R.
validate_derived_quantities(x)validate_derived_quantities(x)
wasserstein distance between binned samples
wasserstein(x, y)wasserstein(x, y)
x |
numeric vector of density from first distribution |
y |
numeric vector of density from second distribution |
... |
unused |
distance value based on max difference