apply self-calib to precondition prior

.Rprofile

@@ -0,0 +1,3 @@
+source("~/.Rprofile")
+# Allows to change how all vignettes are run at once (especially to test rstan)
+options("SBC.vignettes_cmdstanr" = TRUE)

.gitignore

@@ -7,3 +7,5 @@ Meta
 *.exe
 
 *_SBC_cache
+vignettes/stan/*.rds
+vignettes/small_model_workflow/*.rds

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: SBC
 Title: Simulation Based Calibration for rstan/cmdstanr models
-Version: 0.0.1.9000
+Version: 0.1.1.9000
 Authors@R: 
     c(person(given = "Shinyoung",
            family = "Kim",
@@ -16,10 +16,18 @@ Authors@R:
            family = "Modrák",
            role = c("aut"),
            email = "martin.modrak@biomed.cas.cz",
-           comment = c()))
-Description: SBC helps perform Simulation Based Calibration on models using Stan or brms.
-  It contains various classes and functions to extract samples from models, calculate SBC metrics and visualize utility plots.
+           comment = c(ORCID= "0000-0002-8886-7797")),
+    person(given = "Teemu",
+           family = "Säilynoja",
+           role = "aut")
+           )
+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.
 License: MIT + file LICENSE
+URL: https://hyunjimoon.github.io/SBC/, https://github.com/hyunjimoon/SBC/
+BugReports: https://github.com/hyunjimoon/SBC/issues
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
@@ -42,5 +50,6 @@ Suggests:
     rstan,
     knitr,
     rmarkdown,
-    brms
+    brms,
+    mclust
 VignetteBuilder: knitr

NAMESPACE

@@ -6,22 +6,32 @@ S3method(SBC_backend_default_thin_ranks,default)
 S3method(SBC_backend_hash_for_cache,SBC_backend_brms)
 S3method(SBC_backend_hash_for_cache,SBC_backend_cmdstan_sample)
 S3method(SBC_backend_hash_for_cache,SBC_backend_cmdstan_variational)
+S3method(SBC_backend_hash_for_cache,SBC_backend_rjags)
+S3method(SBC_backend_hash_for_cache,SBC_backend_rstan_optimizing)
 S3method(SBC_backend_hash_for_cache,SBC_backend_rstan_sample)
 S3method(SBC_backend_hash_for_cache,default)
-S3method(SBC_backend_iid_samples,SBC_backend_cmdstan_variational)
-S3method(SBC_backend_iid_samples,default)
+S3method(SBC_backend_iid_draws,SBC_backend_cmdstan_variational)
+S3method(SBC_backend_iid_draws,SBC_backend_mock_rng)
+S3method(SBC_backend_iid_draws,SBC_backend_rstan_optimizing)
+S3method(SBC_backend_iid_draws,default)
 S3method(SBC_fit,SBC_backend_brms)
 S3method(SBC_fit,SBC_backend_cmdstan_sample)
 S3method(SBC_fit,SBC_backend_cmdstan_variational)
 S3method(SBC_fit,SBC_backend_mock)
+S3method(SBC_fit,SBC_backend_mock_rng)
+S3method(SBC_fit,SBC_backend_rjags)
+S3method(SBC_fit,SBC_backend_rstan_optimizing)
 S3method(SBC_fit,SBC_backend_rstan_sample)
 S3method(SBC_fit_to_diagnostics,CmdStanMCMC)
 S3method(SBC_fit_to_diagnostics,CmdStanVB)
+S3method(SBC_fit_to_diagnostics,RStanOptimizingFit)
 S3method(SBC_fit_to_diagnostics,brmsfit)
 S3method(SBC_fit_to_diagnostics,default)
 S3method(SBC_fit_to_diagnostics,stanfit)
 S3method(SBC_fit_to_draws_matrix,CmdStanMCMC)
 S3method(SBC_fit_to_draws_matrix,CmdStanVB)
+S3method(SBC_fit_to_draws_matrix,RStanOptimizingFit)
+S3method(SBC_fit_to_draws_matrix,SBC_rjags_fit)
 S3method(SBC_fit_to_draws_matrix,brmsfit)
 S3method(SBC_fit_to_draws_matrix,default)
 S3method(check_all_SBC_diagnostics,SBC_results)
@@ -31,6 +41,8 @@ S3method(generate_datasets,SBC_generator_custom)
 S3method(generate_datasets,SBC_generator_function)
 S3method(get_diagnostic_messages,SBC_ADVI_diagnostics)
 S3method(get_diagnostic_messages,SBC_ADVI_diagnostics_summary)
+S3method(get_diagnostic_messages,SBC_RStanOptimizing_diagnostics)
+S3method(get_diagnostic_messages,SBC_RStanOptimizing_diagnostics_summary)
 S3method(get_diagnostic_messages,SBC_results)
 S3method(get_diagnostic_messages,SBC_results_summary)
 S3method(length,SBC_datasets)
@@ -39,13 +51,16 @@ S3method(plot_contraction,SBC_results)
 S3method(plot_contraction,data.frame)
 S3method(plot_coverage,SBC_results)
 S3method(plot_coverage,data.frame)
+S3method(plot_coverage_diff,SBC_results)
+S3method(plot_coverage_diff,data.frame)
 S3method(plot_rank_hist,SBC_results)
 S3method(plot_rank_hist,data.frame)
 S3method(plot_sim_estimated,SBC_results)
 S3method(plot_sim_estimated,data.frame)
 S3method(print,SBC_nuts_diagnostics_summary)
 S3method(print,SBC_results_summary)
 S3method(summary,SBC_ADVI_diagnostics)
+S3method(summary,SBC_RStanOptimizing_diagnostics)
 S3method(summary,SBC_nuts_diagnostics)
 S3method(summary,SBC_results)
 export(SBC_backend_brms)
@@ -54,25 +69,34 @@ export(SBC_backend_cmdstan_sample)
 export(SBC_backend_cmdstan_variational)
 export(SBC_backend_default_thin_ranks)
 export(SBC_backend_hash_for_cache)
-export(SBC_backend_iid_samples)
+export(SBC_backend_iid_draws)
 export(SBC_backend_mock)
+export(SBC_backend_mock_rng)
+export(SBC_backend_rjags)
+export(SBC_backend_rstan_optimizing)
 export(SBC_backend_rstan_sample)
 export(SBC_datasets)
 export(SBC_diagnostic_messages)
+export(SBC_example_backend)
+export(SBC_example_generator)
+export(SBC_example_results)
 export(SBC_fit)
 export(SBC_fit_to_diagnostics)
 export(SBC_fit_to_draws_matrix)
 export(SBC_generator_brms)
 export(SBC_generator_custom)
 export(SBC_generator_function)
+export(SBC_print_example_model)
 export(SBC_results)
+export(SBC_statistics_from_single_fit)
 export(bind_datasets)
 export(bind_results)
 export(calculate_prior_sd)
 export(calculate_ranks_draws_matrix)
 export(calculate_sds_draws_matrix)
 export(check_all_SBC_diagnostics)
 export(cjs_dist)
+export(compute_SBC)
 export(compute_gen_quants)
 export(compute_results)
 export(data_for_ecdf_plots)
@@ -84,17 +108,26 @@ export(empirical_coverage)
 export(generate_datasets)
 export(generated_quantities)
 export(get_diagnostic_messages)
+export(invtf_param_vec)
 export(max_diff)
 export(plot_contraction)
 export(plot_coverage)
+export(plot_coverage_diff)
 export(plot_ecdf)
 export(plot_ecdf_diff)
 export(plot_rank_hist)
 export(plot_sim_estimated)
 export(rank2unif)
+export(recompute_SBC_statistics)
 export(recompute_statistics)
+export(rnorm_max_coupling)
+export(self_calib_adaptive)
+export(self_calib_gaussian)
+export(self_calib_gmm)
 export(set2set)
-export(statistics_from_single_fit)
+export(tf_param)
+export(tf_param_vec)
+export(update_quantile_approximation)
 export(validate_SBC_datasets)
 export(validate_SBC_results)
 export(validate_generated_quantities)

R/SBC-deprecated.R

@@ -0,0 +1,8 @@
+#' @title Deprecated functions in package \pkg{SBC}.
+#' @description The functions listed below are deprecated and will be defunct in
+#'   the near future. When possible, alternative functions with similar
+#'   functionality are also mentioned. Help pages for deprecated functions are
+#'   available at \code{help("<function>-deprecated")}.
+#' @name SBC-deprecated
+#' @keywords internal
+NULL

R/backend-mock.R

@@ -38,3 +38,42 @@ SBC_fit.SBC_backend_mock <- function(backend, generated, cores) {
   backend$result
 }
 
+#' @export
+SBC_backend_mock_rng <- function(..., n_draws = 1000) {
+  var_to_rng <- list(...)
+  if(is.null(names(var_to_rng)) ||
+     length(unique(names(var_to_rng))) != length(var_to_rng) ||
+     any(names(var_to_rng) == "")
+  ) {
+    stop("All arguments must have a unique name")
+  }
+  var_to_rng <- purrr::map(var_to_rng, purrr::as_mapper)
+
+  purrr::iwalk(var_to_rng, function(rng, name) {
+    tryCatch({
+      res <- rng(13)
+    }, error = function(e) {
+      message(e)
+      stop("Test invocation for argument '", name, "' failed.\n",
+           "All arguments must be convertible to a function that takes the number of draws as input and returns a 1D array of draws.")
+    })
+    if(!is.numeric(res) || length(res) != 13) {
+      stop("Test invocation for argument '", name, "' returned unexpected result.\n",
+           "All arguments must be convertible to a function that takes the number of draws as input and returns a 1D array of draws.")
+    }
+  })
+
+  structure(list(var_to_rng = var_to_rng, n_draws = n_draws), class = "SBC_backend_mock_rng")
+}
+
+#' @export
+SBC_fit.SBC_backend_mock_rng <- function(backend, generated, cores) {
+  draws_list <- purrr::map(backend$var_to_rng, ~ .x(backend$n_draws))
+
+  do.call(posterior::draws_matrix, draws_list)
+}
+
+#' @export
+SBC_backend_iid_draws.SBC_backend_mock_rng <- function(backend) {
+  TRUE
+}

R/backend-rjags.R

@@ -0,0 +1,62 @@
+#' Create a JAGS backend using `rjags`
+#'
+#' @param file model file or connection to model code (passed to [rjags::jags.model()])
+#' @param n.iter number of iterations for sampling (passed to [rjags::coda.samples())
+#' @param n.burnin number of iterations used for burnin
+#' @param variable.names names of variables to monitor (passed to [rjags::coda.samples()])
+#' @param thin thinning (passed to [rjags::coda.samples()])
+#' @param na.rm whether to omit variables containing NA (passed to [rjags::coda.samples()])
+#' @param ... additional optional arguments passed to [rjags::jags.model()]
+#'   - most notably `n.chains`, `n.adapt` and `inits`.
+#' @export
+SBC_backend_rjags <- function(file, n.iter, n.burnin, variable.names, thin = 1, na.rm = TRUE, ...) {
+  args = list(...)
+  if(any(names(args) == "data")) {
+    stop(paste0("Argument 'data' cannot be provided when defining a backend",
+                " as it needs to be set by the SBC package"))
+  }
+
+  structure(list(file = file,
+                 n.iter = n.iter,
+                 variable.names = variable.names,
+                 n.burnin = n.burnin,
+                 thin = thin,
+                 na.rm = na.rm,
+                 args = args), class = "SBC_backend_rjags")
+}
+
+
+#' @export
+SBC_fit.SBC_backend_rjags <- function(backend, generated, cores) {
+  args_all <- c(list(file = backend$file, data = generated), backend$args)
+
+  model <- do.call(rjags::jags.model, args_all)
+  if(backend$n.burnin > 0) {
+    stats::update(model, n.iter = backend$n.burnin, progress.bar = "none")
+  }
+  samples <- rjags::coda.samples(model,
+                          variable.names = backend$variable.names,
+                          n.iter = backend$n.iter,
+                          thin = backend$thin,
+                          na.rm = backend$na.rm,
+                          progress.bar = "none")
+
+  structure(list(model = model, samples = samples),
+            class = "SBC_rjags_fit")
+}
+
+
+#' @export
+SBC_fit_to_draws_matrix.SBC_rjags_fit <- function(fit) {
+  posterior::as_draws_matrix(fit$samples)
+}
+
+#' @export
+SBC_backend_hash_for_cache.SBC_backend_rjags <- function(backend) {
+  model_code <- readLines(backend$file)
+
+  backend_for_cache <- backend
+  backend_for_cache$file <- NULL
+  backend_for_cache$model_code <- model_code
+  rlang::hash(backend_for_cache)
+}