diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index d0cbcc132..5d7ceecb2 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -3,10 +3,19 @@ on: branches: - main - master + paths: + - '.github/workflows/R-CMD-check.yaml' # Run when the workflow file changes + - 'R/**' # Run when R code changes + - 'tests/**' # Run when test code changes pull_request: branches: - main - master + paths: + - '.github/workflows/R-CMD-check.yaml' # Run when the workflow file changes + - 'R/**' + - 'tests/**' + workflow_dispatch: schedule: - cron: '1 23 * * Sun' diff --git a/.github/workflows/geweke-ahmc.yaml b/.github/workflows/geweke-ahmc.yaml new file mode 100644 index 000000000..1bef2ac2d --- /dev/null +++ b/.github/workflows/geweke-ahmc.yaml @@ -0,0 +1,78 @@ +name: geweke-tests-ahmc + +on: + push: + branches: + - main + - master + paths: + - '.github/workflows/geweke-ahmc.yaml' # Run when the workflow file changes + - 'R/**' # Run when R code changes + - 'inst/greta-geweke-ahmc.R' # Run when the test script changes + pull_request: + branches: + - main + - master + paths: + - '.github/workflows/geweke-ahmc.yaml' + - 'R/**' + - 'inst/greta-geweke-ahmc.R' + workflow_dispatch: + +jobs: + run: + runs-on: ubuntu-latest + permissions: + contents: write + timeout-minutes: 360 + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + steps: + - name: Checkout + uses: actions/checkout@v4 + - uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: | + any::rcmdcheck + any::devtools + any::reticulate + any::ggplot2 + any::bayesplot + any::glue + local::. + cache-version: 4 + upgrade: 'TRUE' + - name: Install greta deps + shell: Rscript {0} + run: | + library(greta) + greta::install_greta_deps(timeout = 50) + - name: Situation Report on greta install + shell: Rscript {0} + run: greta::greta_sitrep() + - name: Run Geweke tests + run: Rscript inst/greta-geweke-ahmc.R + + - name: Upload diagnostic plots + uses: actions/upload-artifact@v4 + with: + name: geweke-ahmc-diagnostic-plots + path: geweke-*.png + retention-days: 90 # How long to keep the artifact + + - name: Session info + shell: Rscript {0} + run: | + options(width = 100) + pkgs <- installed.packages()[, "Package"] + sessioninfo::session_info(pkgs, include_base = TRUE) + + # - name: Upload plots and statistics as artifacts + # uses: actions/upload-artifact@v3 + # with: + # name: geweke-results-${{ matrix.os }} + # path: plots/ + # if-no-files-found: warn diff --git a/.github/workflows/geweke-hmc.yaml b/.github/workflows/geweke-hmc.yaml new file mode 100644 index 000000000..e08df6175 --- /dev/null +++ b/.github/workflows/geweke-hmc.yaml @@ -0,0 +1,71 @@ +name: geweke-tests-hmc + +on: + push: + branches: + - main + - master + paths: + - '.github/workflows/geweke-hmc.yaml' # Run when the workflow file changes + - 'R/**' # Run when R code changes + - 'inst/greta-geweke-hmc.R' # Run when the test script changes + pull_request: + branches: + - main + - master + paths: + - '.github/workflows/geweke-hmc.yaml' + - 'R/**' + - 'inst/greta-geweke-hmc.R' + workflow_dispatch: + +jobs: + run: + runs-on: ubuntu-latest + permissions: + contents: write + timeout-minutes: 360 + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + steps: + - name: Checkout + uses: actions/checkout@v4 + - uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: any::rcmdcheck any::devtools any::reticulate local::. + cache-version: 4 + upgrade: 'TRUE' + - name: Install greta deps + shell: Rscript {0} + run: | + library(greta) + greta::install_greta_deps(timeout = 50) + - name: Situation Report on greta install + shell: Rscript {0} + run: greta::greta_sitrep() + - name: Run Geweke tests + run: Rscript inst/greta-geweke-hmc.R + + - name: Upload plot as artifact + uses: actions/upload-artifact@v4 + with: + name: geweke-hmc-diagnostic-plots + path: geweke-*.png + retention-days: 90 # How long to keep the artifact + + - name: Session info + shell: Rscript {0} + run: | + options(width = 100) + pkgs <- installed.packages()[, "Package"] + sessioninfo::session_info(pkgs, include_base = TRUE) + + # - name: Upload plots and statistics as artifacts + # uses: actions/upload-artifact@v3 + # with: + # name: geweke-results-${{ matrix.os }} + # path: plots/ + # if-no-files-found: warn diff --git a/.github/workflows/geweke-rwmh.yaml b/.github/workflows/geweke-rwmh.yaml new file mode 100644 index 000000000..f6f8d5c00 --- /dev/null +++ b/.github/workflows/geweke-rwmh.yaml @@ -0,0 +1,71 @@ +name: geweke-tests-rwmh + +on: + push: + branches: + - main + - master + paths: + - '.github/workflows/geweke-rwmh.yaml' # Run when the workflow file changes + - 'R/**' # Run when R code changes + - 'inst/greta-geweke-rwmh.R' # Run when the test script changes + pull_request: + branches: + - main + - master + paths: + - '.github/workflows/geweke-rwmh.yaml' + - 'R/**' + - 'inst/greta-geweke-rwmh.R' + workflow_dispatch: + +jobs: + run: + runs-on: ubuntu-latest + permissions: + contents: write + timeout-minutes: 360 + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + steps: + - name: Checkout + uses: actions/checkout@v4 + - uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: any::rcmdcheck any::devtools any::reticulate local::. + cache-version: 4 + upgrade: 'TRUE' + - name: Install greta deps + shell: Rscript {0} + run: | + library(greta) + greta::install_greta_deps(timeout = 50) + - name: Situation Report on greta install + shell: Rscript {0} + run: greta::greta_sitrep() + - name: Run Geweke tests + run: Rscript inst/greta-geweke-rwmh.R + + - name: Upload plot as artifact + uses: actions/upload-artifact@v4 + with: + name: geweke-rwmh-diagnostic-plots + path: geweke-*.png + retention-days: 90 # How long to keep the artifact + + - name: Session info + shell: Rscript {0} + run: | + options(width = 100) + pkgs <- installed.packages()[, "Package"] + sessioninfo::session_info(pkgs, include_base = TRUE) + + # - name: Upload plots and statistics as artifacts + # uses: actions/upload-artifact@v3 + # with: + # name: geweke-results-${{ matrix.os }} + # path: plots/ + # if-no-files-found: warn diff --git a/.github/workflows/geweke-slice.yaml b/.github/workflows/geweke-slice.yaml new file mode 100644 index 000000000..2fd2b99c8 --- /dev/null +++ b/.github/workflows/geweke-slice.yaml @@ -0,0 +1,78 @@ +name: geweke-tests-slice + +on: + push: + branches: + - main + - master + paths: + - '.github/workflows/geweke-slice.yaml' # Run when the workflow file changes + - 'R/**' # Run when R code changes + - 'inst/greta-geweke-slice.R' # Run when the test script changes + pull_request: + branches: + - main + - master + paths: + - '.github/workflows/geweke-slice.yaml' + - 'R/**' + - 'inst/greta-geweke-slice.R' + workflow_dispatch: + +jobs: + run: + runs-on: ubuntu-latest + permissions: + contents: write + timeout-minutes: 360 + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + steps: + - name: Checkout + uses: actions/checkout@v4 + - uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: | + any::rcmdcheck + any::devtools + any::reticulate + any::ggplot2 + any::bayesplot + any::glue + local::. + cache-version: 4 + upgrade: 'TRUE' + - name: Install greta deps + shell: Rscript {0} + run: | + library(greta) + greta::install_greta_deps(timeout = 50) + - name: Situation Report on greta install + shell: Rscript {0} + run: greta::greta_sitrep() + - name: Run Geweke tests + run: Rscript inst/greta-geweke-slice.R + + - name: Upload diagnostic plots + uses: actions/upload-artifact@v4 + with: + name: geweke-slice-diagnostic-plots + path: geweke-*.png + retention-days: 90 # How long to keep the artifact + + - name: Session info + shell: Rscript {0} + run: | + options(width = 100) + pkgs <- installed.packages()[, "Package"] + sessioninfo::session_info(pkgs, include_base = TRUE) + + # - name: Upload plots and statistics as artifacts + # uses: actions/upload-artifact@v3 + # with: + # name: geweke-results-${{ matrix.os }} + # path: plots/ + # if-no-files-found: warn diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml index 2468d333b..1c6253b46 100644 --- a/.github/workflows/pkgdown.yaml +++ b/.github/workflows/pkgdown.yaml @@ -3,6 +3,9 @@ on: branches: - main - master + paths: + - 'R/**' # Run when R code changes + - '.github/workflows/pkgdown.yaml' # Run when the workflow file changes name: pkgdown diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index 558fef72a..3db183f39 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -2,9 +2,21 @@ # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help on: push: - branches: [main, master] + branches: + - main + - master + paths: + - '.github/workflows/test-coverage.yaml' # Run when the workflow file changes + - 'R/**' # Run when R code changes + - 'tests/**' # Run when test code changes pull_request: - branches: [main, master] + branches: + - main + - master + paths: + - '.github/workflows/test-coverage.yaml' + - 'R/**' + - 'tests/**' name: test-coverage.yaml diff --git a/NAMESPACE b/NAMESPACE index ede86f97e..ce9a7e573 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -157,6 +157,7 @@ export(adagrad) export(adagrad_da) export(adam) export(adamax) +export(adaptive_hmc) export(apply) export(are_null) export(as.greta_model) diff --git a/R/greta-sitrep.R b/R/greta-sitrep.R index 46f63abdc..4240cd5bd 100644 --- a/R/greta-sitrep.R +++ b/R/greta-sitrep.R @@ -6,12 +6,12 @@ #' @param verbosity character. How verbose the output of the situation report. #' Possible options: "minimal" (default), "detailed", and "quiet". "Minimal" #' provides just information in python version, tensorflow version, -#' tensorflow proability, and whether greta conda environment is available. +#' tensorflow probability, and whether greta conda environment is available. #' "Quiet" presents no information, but prepares greta to be used. "Detailed" #' gives information on the version and path for R, greta, python, #' tensorflow, tensorflow probability, the greta conda environment, and a #' statement on greta usability. -#' @return Message on greta situation report. See "verbsoity" parameter details +#' @return Message on greta situation report. See "verbosity" parameter details #' above for more information. #' @export #' diff --git a/R/package.R b/R/package.R index 3b6df155d..6ceec328d 100644 --- a/R/package.R +++ b/R/package.R @@ -47,6 +47,7 @@ reticulate::py_version #' @importFrom utils globalVariables utils::globalVariables( c( + "as_tensor_spec", "N", "greta_deps_tf_tfp", "greta_logfile", diff --git a/R/sampler_class.R b/R/sampler_class.R index 988228e34..9b0281552 100644 --- a/R/sampler_class.R +++ b/R/sampler_class.R @@ -142,13 +142,15 @@ sampler <- R6Class( # how big would we like the bursts to be ideal_burst_size <- ifelse(one_by_one, 1L, pb_update) - - self$run_warmup( - n_samples = n_samples, - pb_update = pb_update, - ideal_burst_size = ideal_burst_size, - verbose = verbose - ) + do_warmup <- self$warmup > 0 + if (do_warmup) { + self$run_warmup( + n_samples = n_samples, + pb_update = pb_update, + ideal_burst_size = ideal_burst_size, + verbose = verbose + ) + } self$run_sampling( n_samples = n_samples, @@ -565,14 +567,26 @@ sampler <- R6Class( self$free_state <- free_state } + self$update_acceptance(results = batch_results$trace) + self$update_rejection(results = batch_results$trace) + }, + + update_acceptance = function(results) { if (self$uses_metropolis) { # log acceptance probability - log_accept_stats <- as.array(batch_results$trace$log_accept_ratio) - is_accepted <- as.array(batch_results$trace$is_accepted) + log_accept_stats <- as.array(results$log_accept_ratio) + is_accepted <- as.array(results$is_accepted) + # related to tuning information self$accept_history <- rbind(self$accept_history, is_accepted) accept_stats_batch <- pmin(1, exp(log_accept_stats)) + # also related to tuning information self$mean_accept_stat <- mean(accept_stats_batch, na.rm = TRUE) + } + }, + update_rejection = function(results) { + if (self$uses_metropolis) { + log_accept_stats <- as.array(results$log_accept_ratio) # numerical rejections parameter sets bad <- sum(!is.finite(log_accept_stats)) self$numerical_rejections <- self$numerical_rejections + bad @@ -665,3 +679,14 @@ sampler <- R6Class( } ) ) + + +tune_tf <- R6Class( + "tune_tf", + inherit = sampler +) + +tune_r <- R6Class( + "tune_r", + inherit = sampler +) diff --git a/R/samplers.R b/R/samplers.R index 1ad23219c..ac8c3ce37 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -94,6 +94,45 @@ slice <- function(max_doublings = 5) { obj } +#' @rdname samplers +#' +#' @param epsilon leapfrog stepsize hyperparameter (positive, will be tuned) +#' @param diag_sd estimate of the posterior marginal standard deviations +#' (positive, will be tuned). +#' @param max_leapfrog_steps numeric. Default 1000. Maximum number of leapfrog +#' steps used. The algorithm will determine the optimal number less than this. +#' @param method character length one. Currently can only be "SNAPER" but in +#' the future this may expand to other adaptive samplers. +#' @details For `adaptive_hmc()`. The Lmin and Lmax parameters are learned and +#' so not provided in this. The number of chains cannot be less than 2, due +#' to how adaptive HMC works. `diag_sd` is used to rescale the parameter +#' space to make it more uniform, and make sampling more efficient. +#' @export +adaptive_hmc <- function( + max_leapfrog_steps = 1000, + epsilon = 0.1, + diag_sd = 1, + method = "SNAPER" +) { + method <- rlang::arg_match( + arg = method, + values = "SNAPER" + ) + + # nolint end + obj <- list( + parameters = list( + max_leapfrog_steps = max_leapfrog_steps, + epsilon = epsilon, + diag_sd = diag_sd + ), + class = adaptive_hmc_sampler + ) + class(obj) <- c("adaptive_hmc_sampler", "sampler") + obj +} + + #' @noRd #' @export print.sampler <- function(x, ...) { @@ -120,16 +159,6 @@ print.sampler <- function(x, ...) { cat(msg) } -tune_tf <- R6Class( - "tune_tf", - inherit = sampler -) - -tune_r <- R6Class( - "tune_r", - inherit = sampler -) - hmc_sampler <- R6Class( "hmc_sampler", inherit = tune_r, @@ -148,8 +177,9 @@ hmc_sampler <- R6Class( free_state_size <- length(sampler_param_vec) - 2 - # TF1/2 check - # this will likely get replaced... + # TODO we will wrap these up into adaptive sampling to learn better + # similar to the adaptive_hmc - to put more things into TFP, + # rather than going back and forth between R and TF hmc_l <- sampler_param_vec[0] hmc_epsilon <- sampler_param_vec[1] @@ -169,6 +199,12 @@ hmc_sampler <- R6Class( # build the kernel # nolint start + # TODO HMC with momentum would get adaptation addition + # one of the is dual averaging step size adapation (see adaptive_hmc) + # the other is adapt_momentum + # diagonal mass + # https://www.tensorflow.org/probability/api_docs/python/tfp/experimental/mcmc/DiagonalMassMatrixAdaptation + # sampler_kernel <- tfp$mcmc$HamiltonianMonteCarlo( target_log_prob_fn = dag$tf_log_prob_function_adjusted, @@ -318,3 +354,313 @@ slice_sampler <- R6Class( } ) ) + +adaptive_hmc_sampler <- R6Class( + "adaptive_hmc_sampler", + inherit = tune_tf, + public = list( + parameters = list( + # Lmin = 10, + # Lmax = 20, + max_leapfrog_steps = 1000, + # TODO clean up these parameter usage else where + # epsilon = 0.005, + # diag_sd = 1, + # TODO some kind of validity check of method? Currently this can only be + # "SNAPER". + method = "SNAPER" + ), + accept_target = 0.651, + + # TODO + # Eventually we will adapt the other samplers to have an interface + # more in line with this - where we do the adaptation in tensorflow + sampler_kernel = NULL, + + define_tf_kernel = function(sampler_param_vec) { + dag <- self$model$dag + tfe <- dag$tf_environment + + adaptive_hmc_max_leapfrog_steps <- tf$cast( + x = self$parameters$max_leapfrog_steps, + dtype = tf$int32 + ) + + kernel_base <- tfp$experimental$mcmc$SNAPERHamiltonianMonteCarlo( + target_log_prob_fn = dag$tf_log_prob_function_adjusted, + step_size = 1, + num_adaptation_steps = as.integer(self$warmup), + max_leapfrog_steps = adaptive_hmc_max_leapfrog_steps + ) + + # learns the step size + self$sampler_kernel <- tfp$mcmc$DualAveragingStepSizeAdaptation( + inner_kernel = kernel_base, + num_adaptation_steps = as.integer(self$warmup) + ) + }, + + # current_state = NULL, + + # given MCMC kernel `sampler_kernel` and initial model parameter state + # `free_state`, adapt the kernel tuning parameters whilst simultaneously + # burning-in the model parameter state. Return both finalised kernel + # tuning parameters and the burned-in model parameter state + warm_up_sampler = function(sampler_kernel, free_state) { + # get the predetermined adaptation period of the kernel + + # make the uncompiled function (with curried arguments) + warmup_raw <- function() { + tfp$mcmc$sample_chain( + num_results = as.integer(sampler_kernel$num_adaptation_steps), + current_state = free_state, + kernel = sampler_kernel, + return_final_kernel_results = TRUE, + trace_fn = function(current_state, kernel_results) { + kernel_results$step #kernel_results + } + ) + } + + # compile it into a concrete function + warmup <- tensorflow::tf_function(warmup_raw) + + # execute it + result <- warmup() + + # If no warmup, return the free state to use as initials in sampling, + # Otherwise, use the last chunk of all_states, which has been tuned + no_warmup <- as.integer(sampler_kernel$num_adaptation_steps) <= 0 + if (no_warmup) { + self$free_state <- tensorflow::as_tensor(free_state) + } else { + # when there is warmup, get the warmed up result + self$free_state <- get_last_state(result$all_states) + } + # return the last (burned-in) state of the model parameters and the final + # (tuned) kernel parameters + + self$warm_results <- list( + kernel_results = result$final_kernel_results + ) + }, + + sampler_parameter_values = function() { + # random number of integration steps + max_leapfrog_steps <- self$parameters$max_leapfrog_steps + epsilon <- self$parameters$epsilon + diag_sd <- matrix(self$parameters$diag_sd) + method <- self$parameters$method + + # return named list for replacing tensors + list( + adaptive_hmc_max_leapfrog_steps = max_leapfrog_steps, + method = method + ) + }, + + # given a warmed up sampler object, return a compiled TF function + # that generates a new burst of samples from samples from it + sampler_function = NULL, + make_sampler_function = function() { + # make the uncompiled function (with curried arguments) + sample_raw <- function(current_state, n_samples) { + results <- tfp$mcmc$sample_chain( + # how many iterations + num_results = n_samples, + # where to start from + current_state = current_state, + # kernel + kernel = self$sampler_kernel, + # tuned sampler settings + previous_kernel_results = self$warm_results$kernel_results, + # what to trace (nothing) + trace_fn = function(current_state, kernel_results) { + # could compute badness here to save memory? + # is.finite(kernel_results$inner_results$inner_results$inner_results$log_accept_ratio) + kernel_results + } + ) + # return the parameter states and the kernel results + list( + all_states = results$all_states, + kernel_results = results$trace + ) + } + + # compile it into a concrete function and return + self$sampler_function <- tensorflow::tf_function( + sample_raw, + list( + as_tensorspec(tensorflow::as_tensor(self$free_state)), + tf$TensorSpec(shape = c(), dtype = tf$int32) + ) + ) + }, + + # prepare a slot to put the warmed up results into + warm_results = NULL, + + run_warmup = function( + n_samples, + pb_update, + ideal_burst_size, + verbose + ) { + # skip warmup progress bar if warmup is 0 + nonzero_warmup <- self$warmup > 0 + if (verbose && nonzero_warmup) { + pb_warmup <- create_progress_bar( + phase = "warmup", + iter = c(self$warmup, n_samples), + pb_update = self$warmup, + width = self$pb_width + ) + iterate_progress_bar( + pb = pb_warmup, + it = 0, + rejects = 0, + chains = self$n_chains, + file = self$pb_file + ) + } else { + pb_sampling <- NULL + } + self$define_tf_kernel() + self$warm_up_sampler( + sampler_kernel = self$sampler_kernel, + free_state = self$free_state + ) + + # skip warmup progress bar if warmup is 0 + if (verbose && nonzero_warmup) { + # update the progress bar/percentage log + iterate_progress_bar( + pb = pb_warmup, + it = self$warmup, + # TODO check if we can add rejection info + rejects = 0, + chains = self$n_chains, + file = self$pb_file + ) + + self$write_percentage_log( + total = self$warmup, + completed = self$warmup, + stage = "warmup" + ) + } # close verbose + }, + + run_sampling = function( + n_samples, + pb_update, + ideal_burst_size, + trace_batch_size, + thin, + verbose + ) { + perform_sampling <- n_samples > 0 + if (perform_sampling) { + # on exiting during the main sampling period (even if killed by the + # user) trace the free state values + + on.exit(self$trace_values(trace_batch_size), add = TRUE) + + # main sampling + if (verbose) { + pb_sampling <- create_progress_bar( + phase = "sampling", + iter = c(self$warmup, n_samples), + pb_update = pb_update, + width = self$pb_width + ) + iterate_progress_bar( + pb = pb_sampling, + it = 0, + rejects = 0, + chains = self$n_chains, + file = self$pb_file + ) + } else { + pb_sampling <- NULL + } + + ### Adaptive start + # split up warmup iterations into bursts of sampling + burst_lengths <- as.integer(self$burst_lengths( + n_samples, + ideal_burst_size + )) + completed_iterations <- cumsum(burst_lengths) + + # maybe warm up a sampler + if (is.null(self$warm_results)) { + self$run_warmup( + verbose = verbose, + n_samples = n_samples + ) + } + # maybe compile a sampling function + if (is.null(self$sampler_function)) { + self$make_sampler_function() + } + + for (burst in seq_along(burst_lengths)) { + burst_size <- burst_lengths[burst] + batch_results <- self$sampler_function( + current_state = self$free_state, + n_samples = burst_size + ) + + free_state_draws <- as.array(batch_results$all_states) + + # if there is one sample at a time, and it's rejected, conversion from + # python back to R can drop a dimension, so handle that here. Ugh. + if (n_dim(free_state_draws) != 3) { + dim(free_state_draws) <- c(1, dim(free_state_draws)) + } + + self$last_burst_free_states <- split_chains(free_state_draws) + + n_draws <- nrow(free_state_draws) + if (n_draws > 0) { + free_state <- free_state_draws[n_draws, , , drop = FALSE] + dim(free_state) <- dim(free_state)[-1] + self$free_state <- free_state + } + + self$update_rejection(batch_results) + + self$trace() + + if (verbose) { + # update the progress bar/percentage log + iterate_progress_bar( + pb = pb_sampling, + it = completed_iterations[burst], + rejects = self$numerical_rejections, + chains = self$n_chains, + file = self$pb_file + ) + + self$write_percentage_log( + total = n_samples, + completed = completed_iterations[burst], + stage = "sampling" + ) + } + } + } # end sampling + ### Adaptive end + }, + + update_rejection = function(results) { + if (self$uses_metropolis) { + # accumulate and report on the badness + bad <- sum(bad_steps(results$kernel_results)) + self$numerical_rejections <- self$numerical_rejections + bad + } + } + ) +) diff --git a/R/utils.R b/R/utils.R index 679a479ca..2ab79af00 100644 --- a/R/utils.R +++ b/R/utils.R @@ -1079,3 +1079,44 @@ n_warmup <- function(x) { x_info <- attr(x, "model_info") x_info$warmup } + +as_tensorspec <- function(tensor) { + tf$TensorSpec(shape = tensor$shape, dtype = tensor$dtype) +} + +maybe_make_tensor_shape <- function(x) { + if (tf$is_tensor(x)) { + as_tensor_spec(x) + } else { + x + } +} + +# get the final model parameter state from a chain as returned in the all_states +# object from tfp$mcmc$sample_chain +get_last_state <- function(all_states) { + n_iter <- dim(all_states)[1] + tf$gather(all_states, n_iter - 1L, 0L) +} + +# find out if MCMC steps had non-finite acceptance probabilities +bad_steps <- function(kernel_results) { + log_accept_ratios <- recursive_get_log_accept_ratio(kernel_results) + !is.finite(log_accept_ratios) +} + + +# recursively extract the log accaptance ratio from the MCMC kernel +recursive_get_log_accept_ratio <- function(kernel_results) { + nm <- names(kernel_results) + if ("log_accept_ratio" %in% nm) { + log_accept_ratios <- kernel_results$log_accept_ratio + } else if ("inner_results" %in% nm) { + log_accept_ratios <- recursive_get_log_accept_ratio( + kernel_results$inner_results + ) + } else { + stop("non-standard kernel structure") + } + as.array(log_accept_ratios) +} diff --git a/inst/WORDLIST b/inst/WORDLIST index 731a78645..89d589d8c 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -30,23 +30,22 @@ Nesterov NumPy ORCID OpenBUGS -Optimizers PSAT ProximalAdagradOptimizer ProximalGradientDescentOptimizer README RStudio RStudio's +SNAPER Scalable TF TFP Tensorflow -Vectorised WinBUGS XLA -acyclic al analysed +arg args automagically bayesplot @@ -106,7 +105,6 @@ optimisation optimisations optimiser optimisers -optimizers ouputs parallelisation parallelise @@ -114,19 +112,11 @@ parallelising pkgdown poisson polygamma -posteriori -priori pythonic realisation realisations -recoded regularisation reinstalls -repo -reproducibility -rescale -rescaled -rescaling reticulate scalable schoolers @@ -139,12 +129,10 @@ tf tfp tril tuh -unadjusted uncentered uncoached vectorising visualise -warmup winbuilder wishart zhang diff --git a/inst/greta-geweke-ahmc.R b/inst/greta-geweke-ahmc.R new file mode 100644 index 000000000..a23ea1847 --- /dev/null +++ b/inst/greta-geweke-ahmc.R @@ -0,0 +1,88 @@ +devtools::load_all(".") + +# parameters ---- +n <- 10 +mu1 <- 0 +sd1 <- 2 +sd2 <- 1 + +# prior (n draws) ---- +p_theta <- function(n) { + rnorm(n, mu1, sd1) +} + +# likelihood ---- +p_x_bar_theta <- function(theta) { + rnorm(n, theta, sd2) +} + +# define the greta model ---- +x <- as_data(rep(0, n)) +greta_theta <- normal(mu1, sd1) +distribution(x) <- normal(greta_theta, sd2) +model <- model(greta_theta) + +# mcmc parameters ---- +n_iter <- 250 +n_warmup <- 10000 +n_chains <- 20 +n_thin <- 5 +geweke_sampler <- adaptive_hmc() + +# checking ---- +time_taken <- system.time({ + geweke_ahmc <- check_geweke( + sampler = geweke_sampler, + model = model, + data = x, + p_theta = p_theta, + p_x_bar_theta = p_x_bar_theta, + chains = n_chains, + niter = n_iter, + warmup = n_warmup, + thin = n_thin + ) +}) + +# processing ---- +geweke_thin <- apply_thinning(geweke_ahmc, n_thin) + +geweke_stat <- geweke_ks(geweke_thin) + +geweke_stat + +the_time <- round(time_taken[3]) + +cat("This took", the_time, "seconds") + +# ---- plotting ---- +plot_metadata <- build_plot_metadata( + sampler = geweke_sampler, + time = the_time, + n_warmup = n_warmup, + n_iter = n_iter, + n_chains = n_chains, + n_thin = n_thin, + ks_stat = geweke_stat +) + +png(plot_metadata$qq_filename, width = 1200, height = 800, res = 150) + +geweke_qq(geweke_thin, title = plot_metadata$qq_title) +dev.off() + +png(plot_metadata$coda_filename, width = 1200, height = 800, res = 150) +plot(geweke_ahmc$draws) +dev.off() + +bplot_dens <- bayesplot::mcmc_dens_chains(geweke_ahmc$draws) + +ggplot2::ggsave( + filename = plot_metadata$bayesplot_filename, + plot = bplot_dens, + width = 1200, + height = 800, + dpi = 150, + units = "px", + bg = "white" +) diff --git a/inst/greta-geweke-hmc.R b/inst/greta-geweke-hmc.R new file mode 100644 index 000000000..7a577b8c4 --- /dev/null +++ b/inst/greta-geweke-hmc.R @@ -0,0 +1,88 @@ +devtools::load_all(".") + +# parameters ---- +n <- 10 +mu1 <- 0 +sd1 <- 2 +sd2 <- 1 + +# prior (n draws) ---- +p_theta <- function(n) { + rnorm(n, mu1, sd1) +} + +# likelihood ---- +p_x_bar_theta <- function(theta) { + rnorm(n, theta, sd2) +} + +# define the greta model ---- +x <- as_data(rep(0, n)) +greta_theta <- normal(mu1, sd1) +distribution(x) <- normal(greta_theta, sd2) +model <- model(greta_theta) + +# mcmc parameters ---- +n_iter <- 2000 +n_warmup <- 2000 +n_chains <- 10 +n_thin <- 5 +geweke_sampler <- hmc() + +# checking ---- +time_taken <- system.time({ + geweke_ahmc <- check_geweke( + sampler = geweke_sampler, + model = model, + data = x, + p_theta = p_theta, + p_x_bar_theta = p_x_bar_theta, + chains = n_chains, + niter = n_iter, + warmup = n_warmup, + thin = n_thin + ) +}) + +# processing ---- +geweke_thin <- apply_thinning(geweke_ahmc, n_thin) + +geweke_stat <- geweke_ks(geweke_thin) + +geweke_stat + +the_time <- round(time_taken[3]) + +cat("This took", the_time, "seconds") + +# ---- plotting ---- +plot_metadata <- build_plot_metadata( + sampler = geweke_sampler, + time = the_time, + n_warmup = n_warmup, + n_iter = n_iter, + n_chains = n_chains, + n_thin = n_thin, + ks_stat = geweke_stat +) + +png(plot_metadata$qq_filename, width = 1200, height = 800, res = 150) + +geweke_qq(geweke_thin, title = plot_metadata$qq_title) +dev.off() + +png(plot_metadata$coda_filename, width = 1200, height = 800, res = 150) +plot(geweke_ahmc$draws) +dev.off() + +bplot_dens <- bayesplot::mcmc_dens_chains(geweke_ahmc$draws) + +ggplot2::ggsave( + filename = plot_metadata$bayesplot_filename, + plot = bplot_dens, + width = 1200, + height = 800, + dpi = 150, + units = "px", + bg = "white" +) diff --git a/inst/greta-geweke-rwmh.R b/inst/greta-geweke-rwmh.R new file mode 100644 index 000000000..c4d41beec --- /dev/null +++ b/inst/greta-geweke-rwmh.R @@ -0,0 +1,88 @@ +devtools::load_all(".") + +# parameters ---- +n <- 10 +mu1 <- 0 +sd1 <- 2 +sd2 <- 1 + +# prior (n draws) ---- +p_theta <- function(n) { + rnorm(n, mu1, sd1) +} + +# likelihood ---- +p_x_bar_theta <- function(theta) { + rnorm(n, theta, sd2) +} + +# define the greta model ---- +x <- as_data(rep(0, n)) +greta_theta <- normal(mu1, sd1) +distribution(x) <- normal(greta_theta, sd2) +model <- model(greta_theta) + +# mcmc parameters ---- +n_iter <- 2000 +n_warmup <- 2000 +n_chains <- 2 +n_thin <- 5 +geweke_sampler <- rwmh() + +# checking ---- +time_taken <- system.time({ + geweke_ahmc <- check_geweke( + sampler = geweke_sampler, + model = model, + data = x, + p_theta = p_theta, + p_x_bar_theta = p_x_bar_theta, + chains = n_chains, + niter = n_iter, + warmup = n_warmup, + thin = n_thin + ) +}) + +# processing ---- +geweke_thin <- apply_thinning(geweke_ahmc, n_thin) + +geweke_stat <- geweke_ks(geweke_thin) + +geweke_stat + +the_time <- round(time_taken[3]) + +cat("This took", the_time, "seconds") + +# ---- plotting ---- +plot_metadata <- build_plot_metadata( + sampler = geweke_sampler, + time = the_time, + n_warmup = n_warmup, + n_iter = n_iter, + n_chains = n_chains, + n_thin = n_thin, + ks_stat = geweke_stat +) + +png(plot_metadata$qq_filename, width = 1200, height = 800, res = 150) + +geweke_qq(geweke_thin, title = plot_metadata$qq_title) +dev.off() + +png(plot_metadata$coda_filename, width = 1200, height = 800, res = 150) +plot(geweke_ahmc$draws) +dev.off() + +bplot_dens <- bayesplot::mcmc_dens_chains(geweke_ahmc$draws) + +ggplot2::ggsave( + filename = plot_metadata$bayesplot_filename, + plot = bplot_dens, + width = 1200, + height = 800, + dpi = 150, + units = "px", + bg = "white" +) diff --git a/inst/greta-geweke-slice.R b/inst/greta-geweke-slice.R new file mode 100644 index 000000000..d40caf87c --- /dev/null +++ b/inst/greta-geweke-slice.R @@ -0,0 +1,88 @@ +devtools::load_all(".") + +# parameters ---- +n <- 10 +mu1 <- 0 +sd1 <- 2 +sd2 <- 1 + +# prior (n draws) ---- +p_theta <- function(n) { + rnorm(n, mu1, sd1) +} + +# likelihood ---- +p_x_bar_theta <- function(theta) { + rnorm(n, theta, sd2) +} + +# define the greta model ---- +x <- as_data(rep(0, n)) +greta_theta <- normal(mu1, sd1) +distribution(x) <- normal(greta_theta, sd2) +model <- model(greta_theta) + +# mcmc parameters ---- +n_iter <- 2000 +n_warmup <- 2000 +n_chains <- 2 +n_thin <- 1 +geweke_sampler <- slice() + +# checking ---- +time_taken <- system.time({ + geweke_ahmc <- check_geweke( + sampler = geweke_sampler, + model = model, + data = x, + p_theta = p_theta, + p_x_bar_theta = p_x_bar_theta, + chains = n_chains, + niter = n_iter, + warmup = n_warmup, + thin = n_thin + ) +}) + +# processing ---- +geweke_thin <- apply_thinning(geweke_ahmc, n_thin) + +geweke_stat <- geweke_ks(geweke_thin) + +geweke_stat + +the_time <- round(time_taken[3]) + +cat("This took", the_time, "seconds") + +# ---- plotting ---- +plot_metadata <- build_plot_metadata( + sampler = geweke_sampler, + time = the_time, + n_warmup = n_warmup, + n_iter = n_iter, + n_chains = n_chains, + n_thin = n_thin, + ks_stat = geweke_stat +) + +png(plot_metadata$qq_filename, width = 1200, height = 800, res = 150) + +geweke_qq(geweke_thin, title = plot_metadata$qq_title) +dev.off() + +png(plot_metadata$coda_filename, width = 1200, height = 800, res = 150) +plot(geweke_ahmc$draws) +dev.off() + +bplot_dens <- bayesplot::mcmc_dens_chains(geweke_ahmc$draws) + +ggplot2::ggsave( + filename = plot_metadata$bayesplot_filename, + plot = bplot_dens, + width = 1200, + height = 800, + dpi = 150, + units = "px", + bg = "white" +) diff --git a/man/greta_sitrep.Rd b/man/greta_sitrep.Rd index 8fd90b953..79716fad4 100644 --- a/man/greta_sitrep.Rd +++ b/man/greta_sitrep.Rd @@ -10,14 +10,14 @@ greta_sitrep(verbosity = c("minimal", "detailed", "quiet")) \item{verbosity}{character. How verbose the output of the situation report. Possible options: "minimal" (default), "detailed", and "quiet". "Minimal" provides just information in python version, tensorflow version, -tensorflow proability, and whether greta conda environment is available. +tensorflow probability, and whether greta conda environment is available. "Quiet" presents no information, but prepares greta to be used. "Detailed" gives information on the version and path for R, greta, python, tensorflow, tensorflow probability, the greta conda environment, and a statement on greta usability.} } \value{ -Message on greta situation report. See "verbsoity" parameter details +Message on greta situation report. See "verbosity" parameter details above for more information. } \description{ diff --git a/man/samplers.Rd b/man/samplers.Rd index 076eded3b..9a7e9f640 100644 --- a/man/samplers.Rd +++ b/man/samplers.Rd @@ -5,6 +5,7 @@ \alias{hmc} \alias{rwmh} \alias{slice} +\alias{adaptive_hmc} \title{MCMC samplers} \usage{ hmc(Lmin = 5, Lmax = 10, epsilon = 0.1, diag_sd = 1) @@ -12,6 +13,13 @@ hmc(Lmin = 5, Lmax = 10, epsilon = 0.1, diag_sd = 1) rwmh(proposal = c("normal", "uniform"), epsilon = 0.1, diag_sd = 1) slice(max_doublings = 5) + +adaptive_hmc( + max_leapfrog_steps = 1000, + epsilon = 0.1, + diag_sd = 1, + method = "SNAPER" +) } \arguments{ \item{Lmin}{minimum number of leapfrog steps (positive integer, Lmin > Lmax)} @@ -27,6 +35,12 @@ slice(max_doublings = 5) \item{max_doublings}{the maximum number of iterations of the 'doubling' algorithm used to adapt the size of the slice} + +\item{max_leapfrog_steps}{numeric. Default 1000. Maximum number of leapfrog +steps used. The algorithm will determine the optimal number less than this.} + +\item{method}{character length one. Currently can only be "SNAPER" but in +the future this may expand to other adaptive samplers.} } \value{ a \code{sampler} object that can be passed to \code{\link[=mcmc]{mcmc()}}. @@ -53,4 +67,9 @@ is implemented for uniform and normal proposals. \code{slice()} implements a multivariate slice sampling algorithm. The parameter \code{max_doublings} is not tuned during warmup. + +For \code{adaptive_hmc()}. The Lmin and Lmax parameters are learned and +so not provided in this. The number of chains cannot be less than 2, due +to how adaptive HMC works. \code{diag_sd} is used to rescale the parameter +space to make it more uniform, and make sampling more efficient. } diff --git a/tests/testthat/_snaps/adaptive-hmc.md b/tests/testthat/_snaps/adaptive-hmc.md new file mode 100644 index 000000000..2b3d1601c --- /dev/null +++ b/tests/testthat/_snaps/adaptive-hmc.md @@ -0,0 +1,34 @@ +# adaptive_hmc() errors when given 1 chain + + Code + draws <- mcmc(m, n_samples = 1, warmup = 1, chains = 1, sampler = adaptive_hmc()) + Message + + warmup 0/1 | eta: ?s + Condition + Error in `py_call_impl()`: + ! ValueError: SNAPERHMC requires at least 2 chains. Got: 1 + Run `reticulate::py_last_error()` for details. + +# bad mcmc proposals are rejected + + Code + draws <- mcmc(m, chains = 2, n_samples = 2, warmup = 0, verbose = FALSE, + sampler = adaptive_hmc(), initial_values = initials(z = 1e+120)) + Message + Only one set of initial values was provided, and was used for all chains + Condition + Error in `self$check_valid_parameters()`: + ! The log density could not be evaluated at these initial values + Try using these initials as the `values` argument in `calculate()` to see what values of subsequent s these initial values lead to. + +--- + + Code + mcmc(m, chains = 1, n_samples = 1, warmup = 0, sampler = adaptive_hmc(), + verbose = FALSE) + Condition + Error in `self$check_reasonable_starting_values()`: + ! Could not find reasonable starting values after 20 attempts. + Please specify initial values manually via the `initial_values` argument + diff --git a/tests/testthat/_snaps/greta-sitrep.md b/tests/testthat/_snaps/greta-sitrep.md index a0ac6579a..10031eb79 100644 --- a/tests/testthat/_snaps/greta-sitrep.md +++ b/tests/testthat/_snaps/greta-sitrep.md @@ -246,18 +246,18 @@ Message -- R --------------------------------------------------------------------------- - * version: 4.4.2 + * version: 4.5.0 * path: '/Library/Frameworks/R.framework/Resources' -- greta ----------------------------------------------------------------------- * version: 0.5.0.9000 - * path: '/Users/nick/github/greta-dev/greta' + * path: '/Users/nick_1/github/greta-dev/greta' -- python ---------------------------------------------------------------------- i checking if python available v python (v3.10) available - * path: '/Users/nick/Library/r-miniconda-arm64' + * path: '/Users/nick_1/Library/r-miniconda-arm64' -- greta conda environment ----------------------------------------------------- i checking if greta conda environment available @@ -269,7 +269,7 @@ i checking if TensorFlow available v TensorFlow (v2.15.0) available - * R path: '/Users/nick/Library/R/arm64/4.4/library/tensorflow' + * R path: '/Users/nick_1/Library/R/arm64/4.5/library/tensorflow' * Exists in conda env: TRUE -- TensorFlow Probability ------------------------------------------------------ diff --git a/tests/testthat/_snaps/inference.md b/tests/testthat/_snaps/inference.md index 1fa4c7f77..fcf305e33 100644 --- a/tests/testthat/_snaps/inference.md +++ b/tests/testthat/_snaps/inference.md @@ -218,3 +218,11 @@ hmc sampler object with parameters: Lmin = 1, Lmax = 10, epsilon = 0.1, diag_sd = 1 +--- + + Code + adaptive_hmc(max_leapfrog_steps = 1L) + Output + adaptive_hmc_sampler object with parameters: + max_leapfrog_steps = 1, epsilon = 0.1, diag_sd = 1 + diff --git a/tests/testthat/helpers.R b/tests/testthat/helpers.R index fe589749c..cfc75fb6b 100644 --- a/tests/testthat/helpers.R +++ b/tests/testthat/helpers.R @@ -661,7 +661,8 @@ check_geweke <- function( p_x_bar_theta, niter = 2000, warmup = 1000, - thin = 1 + thin = 1, + chains = 1 ) { # sample independently target_theta <- p_theta(niter) @@ -674,12 +675,14 @@ check_geweke <- function( p_theta = p_theta, p_x_bar_theta = p_x_bar_theta, sampler = sampler, - warmup = warmup + warmup = warmup, + chains = chains ) geweke_checks <- list( - target_theta = do_thinning(target_theta, thin), - greta_theta = do_thinning(greta_theta, thin) + target_theta = target_theta, + greta_theta = greta_theta$theta, + draws = greta_theta$draws ) geweke_checks @@ -688,9 +691,13 @@ check_geweke <- function( geweke_qq <- function(geweke_checks, title) { # visualise correspondence quants <- (1:99) / 100 + # target q1 <- stats::quantile(geweke_checks$target_theta, quants) + # greta sampled q2 <- stats::quantile(geweke_checks$greta_theta, quants) - plot(q2, q1, main = title) + # q1 is target + # q2 is greta + plot(q2, q1, main = title, xlab = "Q2 greta", ylab = "Q1 target") graphics::abline(0, 1) } @@ -714,8 +721,10 @@ p_theta_greta <- function( data, p_theta, p_x_bar_theta, + # this argument gets passed along based on sampler used sampler = hmc(), - warmup = 1000 + warmup = 1000, + chains ) { # set up and initialize trace theta <- rep(NA, niter) @@ -726,7 +735,7 @@ p_theta_greta <- function( model, warmup = warmup, n_samples = 1, - chains = 1, + chains = chains, sampler = sampler, verbose = FALSE ) @@ -752,24 +761,37 @@ p_theta_greta <- function( # samples, now using this value of x (slow, but necessary in eager mode) dag$tf_log_prob_function <- NULL dag$define_tf_log_prob_function() + sampler <- attr(draws, "model_info")$samplers[[1]] - sampler$define_tf_evaluate_sample_batch() - # take anoteher sample + # we need a condition inside the geweke test where if it's adaptive_hmc we + # recreate the log prob function, recreate the default kernel, replace that + # inside the warmed up results, and then recreate the sampler function + if (inherits(sampler, "adaptive_hmc_sampler")) { + sampler$define_tf_kernel() + sampler$make_sampler_function() + } else { + sampler$define_tf_evaluate_sample_batch() + } + + # take another sample draws <- extra_samples( draws, n_samples = 1, verbose = FALSE ) - # trace the sample + # trace the sample - just one chain theta[i] <- tail(as.numeric(draws[[1]]), 1) } # kill the progress_bar cli::cli_progress_done() - theta + list( + theta = theta, + draws = draws + ) } # test mcmc for models with analytic posteriors @@ -915,6 +937,63 @@ do_thinning <- function(x, thinning = 1) { x[idx] } +apply_thinning <- function(geweke_result, n_thin) { + geweke_thin <- list( + target_theta = do_thinning(geweke_result$target_theta, n_thin), + greta_theta = do_thinning(geweke_result$greta_theta, n_thin) + ) + geweke_thin +} + +build_plot_metadata <- function( + sampler, + time, + n_warmup, + n_iter, + n_chains, + n_thin, + ks_stat +) { + sampler_class <- class(sampler) + + sampler <- switch( + sampler_class[1], + "adaptive_hmc_sampler" = "ahmc", + "slice sampler" = "slice", + "rwmh sampler" = "rwmh", + "hmc sampler" = "hmc" + ) + + qq_title <- glue::glue( + "Sampler {sampler} Geweke test (time taken = {time}s)\\ + \nwarmup = {n_warmup}, iter = {n_iter}, chains = {n_chains}, thin = {n_thin}\np-value = {round(ks_stat$p.value, 8)}" + ) + + param_stub <- glue::glue( + "warm-{n_warmup}-it-{n_iter}\\ + -chain-{n_chains}-thin-{n_thin}.png" + ) + + qq_filename <- glue::glue( + "geweke-{sampler}-qq-{param_stub}" + ) + + coda_filename <- glue::glue( + "geweke-{sampler}-coda-{param_stub}" + ) + + bayesplot_filename <- glue::glue( + "geweke-{sampler}-dens-{param_stub}" + ) + + list( + qq_title = qq_title, + qq_filename = qq_filename, + coda_filename = coda_filename, + bayesplot_filename = bayesplot_filename + ) +} + get_distribution_name <- function(x) { x_node <- get_node(x) @@ -934,6 +1013,7 @@ get_distribution_name <- function(x) { check_samples <- function( x, iid_function, + # TODO note that we might want to change this to adaptive_hmc sampler = hmc(), n_effective = 3000, title = NULL, diff --git a/tests/testthat/test-adaptive-hmc.R b/tests/testthat/test-adaptive-hmc.R new file mode 100644 index 000000000..825db0bfa --- /dev/null +++ b/tests/testthat/test-adaptive-hmc.R @@ -0,0 +1,106 @@ +set.seed(2025 - 02 - 13) + +test_that("adaptive_hmc() errors when given 1 chain", { + skip_if_not(check_tf_version()) + x <- normal(0, 1) + m <- model(x) + expect_snapshot( + error = TRUE, + draws <- mcmc( + m, + n_samples = 1, + warmup = 1, + chains = 1, + sampler = adaptive_hmc() + ) + ) +}) + +test_that("bad mcmc proposals are rejected", { + skip_if_not(check_tf_version()) + + # set up for numerical rejection of initial location + x <- rnorm(10000, 1e60, 1) + z <- normal(-1e60, 1e-60) + distribution(x) <- normal(z, 1e-60) + m <- model(z, precision = "single") + + # # catch badness in the progress bar + out <- get_output( + mcmc( + m, + n_samples = 10, + warmup = 0, + pb_update = 10, + sampler = adaptive_hmc() + ) + ) + + expect_match(out, "100% bad") + + expect_snapshot( + error = TRUE, + draws <- mcmc( + m, + chains = 2, + n_samples = 2, + warmup = 0, + verbose = FALSE, + sampler = adaptive_hmc(), + initial_values = initials(z = 1e120) + ) + ) + + # really bad proposals + x <- rnorm(100000, 1e120, 1) + z <- normal(-1e120, 1e-120) + distribution(x) <- normal(z, 1e-120) + m <- model(z, precision = "single") + expect_snapshot( + error = TRUE, + mcmc( + m, + chains = 1, + n_samples = 1, + warmup = 0, + sampler = adaptive_hmc(), + verbose = FALSE + ) + ) + + # proposals that are fine, but rejected anyway + z <- normal(0, 1) + m <- model(z, precision = "single") + expect_ok(mcmc( + m, + adaptive_hmc( + epsilon = 100, + ), + chains = 2, + n_samples = 5, + warmup = 0, + verbose = FALSE + )) +}) + +test_that("extra_samples works with adaptive_hmc()", { + skip_if_not(check_tf_version()) + + # set up model + a <- normal(0, 1) + m <- model(a) + + draws <- mcmc( + m, + warmup = 10, + n_samples = 10, + verbose = FALSE, + sampler = adaptive_hmc() + ) + + more_draws <- extra_samples(draws, 20, verbose = FALSE) + + expect_true(inherits(more_draws, "greta_mcmc_list")) + expect_true(coda::niter(more_draws) == 30) + expect_true(coda::nchain(more_draws) == 4) +}) diff --git a/tests/testthat/test_inference.R b/tests/testthat/test_inference.R index 4f883a045..a76e8babb 100644 --- a/tests/testthat/test_inference.R +++ b/tests/testthat/test_inference.R @@ -422,6 +422,9 @@ test_that("samplers print informatively", { expect_snapshot( hmc(Lmin = 1) ) + expect_snapshot( + adaptive_hmc(max_leapfrog_steps = 1L) + ) # # check print sees changed parameters # out <- capture_output(hmc(Lmin = 1), TRUE) diff --git a/tests/testthat/test_posteriors_binomial.R b/tests/testthat/test_posteriors_binomial.R index 92d4a0879..5772e7d79 100644 --- a/tests/testthat/test_posteriors_binomial.R +++ b/tests/testthat/test_posteriors_binomial.R @@ -29,3 +29,35 @@ test_that("posterior is correct (binomial)", { suppressWarnings(test <- ks.test(samples, comparison)) expect_gte(test$p.value, 0.01) }) + +test_that("posterior is correct (binomial) with adaptive_hmc", { + skip_if_not(check_tf_version()) + skip_on_os("windows") + # analytic solution to the posterior of the paramter of a binomial + # distribution, with uniform prior + n <- 100 + pos <- rbinom(1, n, runif(1)) + theta <- uniform(0, 1) + distribution(pos) <- binomial(n, theta) + m <- model(theta) + + draws <- get_enough_draws(m, adaptive_hmc(), 2000, verbose = FALSE) + + samples <- as.matrix(draws) + + # analytic solution to posterior is beta(1 + pos, 1 + N - pos) + shape1 <- 1 + pos + shape2 <- 1 + n - pos + + # qq plot against true quantiles + quants <- (1:99) / 100 + q_target <- qbeta(quants, shape1, shape2) + q_est <- quantile(samples, quants) + plot(q_target ~ q_est, main = "binomial posterior") + abline(0, 1) + + n_draws <- round(coda::effectiveSize(draws)) + comparison <- rbeta(n_draws, shape1, shape2) + suppressWarnings(test <- ks.test(samples, comparison)) + expect_gte(test$p.value, 0.01) +}) diff --git a/tests/testthat/test_posteriors_chi_squared.R b/tests/testthat/test_posteriors_chi_squared.R index 1332d347b..63d9e8423 100644 --- a/tests/testthat/test_posteriors_chi_squared.R +++ b/tests/testthat/test_posteriors_chi_squared.R @@ -19,3 +19,25 @@ test_that("samplers are unbiased for chi-squared", { expect_gte(stat$p.value, 0.01) }) + +test_that("samplers are unbiased for chi-squared", { + skip_if_not(check_tf_version()) + skip_on_os("windows") + df <- 5 + x <- chi_squared(df) + iid <- function(n) rchisq(n, df) + + chi_squared_checked <- check_samples( + x = x, + iid_function = iid, + sampler = adaptive_hmc() + ) + + # do the plotting + qqplot_checked_samples(chi_squared_checked) + + # do a formal hypothesis test + stat <- ks_test_mcmc_vs_iid(chi_squared_checked) + + expect_gte(stat$p.value, 0.01) +}) diff --git a/tests/testthat/test_posteriors_geweke.R b/tests/testthat/test_posteriors_geweke.R index 40a6cac5c..8df71b6d5 100644 --- a/tests/testthat/test_posteriors_geweke.R +++ b/tests/testthat/test_posteriors_geweke.R @@ -1,38 +1,36 @@ Sys.setenv("RELEASE_CANDIDATE" = "false") -test_that("samplers pass geweke tests", { +# run geweke tests on this model: +# theta ~ normal(mu1, sd1) +# x[i] ~ normal(theta, sd2) +# for i in N + +n <- 10 + +mu1 <- 0 +sd1 <- 2 +sd2 <- 1 + +# prior (n draws) +p_theta <- function(n) { + rnorm(n, mu1, sd1) +} + +# likelihood +p_x_bar_theta <- function(theta) { + rnorm(n, theta, sd2) +} + +# define the greta model (single precision for slice sampler) +x <- as_data(rep(0, n)) +greta_theta <- normal(mu1, sd1) +distribution(x) <- normal(greta_theta, sd2) +model <- model(greta_theta, precision = "single") + +test_that("hmc sampler passes geweke tests", { skip_if_not(check_tf_version()) - skip_if_not_release() - # nolint start - # run geweke tests on this model: - # theta ~ normal(mu1, sd1) - # x[i] ~ normal(theta, sd2) - # for i in N - # nolint end - - n <- 10 - mu1 <- rnorm(1, 0, 3) - sd1 <- rlnorm(1) - sd2 <- rlnorm(1) - - # prior (n draws) - p_theta <- function(n) { - rnorm(n, mu1, sd1) - } - - # likelihood - p_x_bar_theta <- function(theta) { - rnorm(n, theta, sd2) - } - - # define the greta model (single precision for slice sampler) - x <- as_data(rep(0, n)) - greta_theta <- normal(mu1, sd1) - distribution(x) <- normal(greta_theta, sd2) - model <- model(greta_theta, precision = "single") - # run tests on all available samplers geweke_hmc <- check_geweke( sampler = hmc(), @@ -48,8 +46,13 @@ test_that("samplers pass geweke tests", { geweke_stat_hmc <- geweke_ks(geweke_hmc) testthat::expect_gte(geweke_stat_hmc$p.value, 0.005) +}) - geweke_hmc_rwmh <- check_geweke( +test_that("rwmh sampler passes geweke tests", { + skip_if_not(check_tf_version()) + skip_if_not_release() + + geweke_rwmh <- check_geweke( sampler = rwmh(), model = model, data = x, @@ -59,13 +62,19 @@ test_that("samplers pass geweke tests", { thin = 5 ) - geweke_qq(geweke_hmc_rwmh, title = "RWMH Geweke test") + geweke_qq(geweke_rwmh, title = "RWMH Geweke test") - geweke_stat_rwmh <- geweke_ks(geweke_hmc_rwmh) + geweke_stat_rwmh <- geweke_ks(geweke_rwmh) testthat::expect_gte(geweke_stat_rwmh$p.value, 0.005) +}) - geweke_hmc_slice <- check_geweke( +test_that("slice sampler passes geweke tests", { + skip_if_not(check_tf_version()) + + skip_if_not_release() + + geweke_slice <- check_geweke( sampler = slice(), model = model, data = x, @@ -73,7 +82,33 @@ test_that("samplers pass geweke tests", { p_x_bar_theta = p_x_bar_theta ) - geweke_qq(geweke_hmc_slice, title = "slice sampler Geweke test") + geweke_qq(geweke_slice, title = "slice sampler Geweke test") + + geweke_stat_slice <- geweke_ks(geweke_slice) + + testthat::expect_gte(geweke_stat_slice$p.value, 0.005) +}) + +test_that("adaptive hmc sampler passes geweke tests", { + skip_if_not(check_tf_version()) + + skip_if_not_release() + + geweke_adaptive_hmc <- check_geweke( + sampler = adaptive_hmc(), + model = model, + data = x, + p_theta = p_theta, + p_x_bar_theta = p_x_bar_theta, + chains = 10, + niter = 1000, + warmup = 10000, + thin = 5 + ) + + geweke_qq(geweke_adaptive_hmc, title = "adaptive hmc sampler Geweke test") + + geweke_stat_adaptive_hmc <- geweke_ks(geweke_adaptive_hmc) - testthat::expect_gte(geweke_hmc_slice$p.value, 0.005) + testthat::expect_gte(geweke_stat_adaptive_hmc$p.value, 0.005) })