diff --git a/R/adaptation.R b/R/adaptation.R index 2b81ced..ec120b8 100644 --- a/R/adaptation.R +++ b/R/adaptation.R @@ -9,8 +9,8 @@ #' @param kappa Decay rate exponent in `[0.5, 1]` for adaptation learning rate. #' #' @return List of functions with entries -#' * `initialize`, a function for initializing adapter state at beginning of -#' chain, +#' * `initialize`, a function for initializing adapter state and proposal +#' parameters at beginning of chain, #' * `update` a function for updating adapter state and proposal parameters on #' each chain iteration, #' * `finalize` a function for performing any final updates to adapter state and @@ -27,24 +27,22 @@ #' grad_log_density = function(x) -x #' ) #' proposal <- barker_proposal(target_distribution) -#' adapter <- scale_adapter( -#' proposal, -#' initial_scale = 1., target_accept_prob = 0.4 -#' ) +#' adapter <- scale_adapter(initial_scale = 1., target_accept_prob = 0.4) +#' adapter$initialize(proposal, chain_state(c(0, 0))) scale_adapter <- function( - proposal, initial_scale = NULL, target_accept_prob = NULL, kappa = 0.6) { + initial_scale = NULL, target_accept_prob = NULL, kappa = 0.6) { log_scale <- NULL - if (is.null(target_accept_prob)) { - target_accept_prob <- proposal$default_target_accept_prob() - } - initialize <- function(initial_state) { + initialize <- function(proposal, initial_state) { if (is.null(initial_scale)) { initial_scale <- proposal$default_initial_scale(initial_state$dimension()) } log_scale <<- log(initial_scale) proposal$update(scale = initial_scale) } - update <- function(sample_index, state_and_statistics) { + update <- function(proposal, sample_index, state_and_statistics) { + if (is.null(target_accept_prob)) { + target_accept_prob <- proposal$default_target_accept_prob() + } gamma <- sample_index^(-kappa) accept_prob <- state_and_statistics$statistics$accept_prob log_scale <<- log_scale + gamma * (accept_prob - target_accept_prob) @@ -53,7 +51,7 @@ scale_adapter <- function( list( initialize = initialize, update = update, - finalize = function() {}, + finalize = NULL, state = function() list(log_scale = log_scale) ) } @@ -74,15 +72,16 @@ scale_adapter <- function( #' grad_log_density = function(x) -x #' ) #' proposal <- barker_proposal(target_distribution) -#' adapter <- variance_adapter(proposal) -variance_adapter <- function(proposal, kappa = 0.6) { +#' adapter <- variance_adapter() +#' adapter$initialize(proposal, chain_state(c(0, 0))) +variance_adapter <- function(kappa = 0.6) { mean_estimate <- NULL variance_estimate <- NULL - initialize <- function(initial_state) { + initialize <- function(proposal, initial_state) { mean_estimate <<- initial_state$position() variance_estimate <<- rep(1., initial_state$dimension()) } - update <- function(sample_index, state_and_statistics) { + update <- function(proposal, sample_index, state_and_statistics) { gamma <- sample_index^(-kappa) position <- state_and_statistics$state$position() mean_estimate <<- mean_estimate + gamma * (position - mean_estimate) @@ -124,20 +123,23 @@ variance_adapter <- function(proposal, kappa = 0.6) { #' grad_log_density = function(x) -x #' ) #' proposal <- barker_proposal(target_distribution) -#' adapter <- robust_shape_adapter( -#' proposal, -#' initial_scale = 1., -#' target_accept_prob = 0.4 -#' ) +#' adapter <- robust_shape_adapter(initial_scale = 1., target_accept_prob = 0.4) +#' adapter$initialize(proposal, chain_state(c(0, 0))) robust_shape_adapter <- function( - proposal, initial_scale, target_accept_prob = 0.4, kappa = 0.6) { + initial_scale = NULL, target_accept_prob = NULL, kappa = 0.6) { rlang::check_installed("ramcmc", reason = "to use this function") shape <- NULL - initialize <- function(initial_state) { + initialize <- function(proposal, initial_state) { + if (is.null(initial_scale)) { + initial_scale <- proposal$default_initial_scale(initial_state$dimension()) + } shape <<- initial_scale * diag(initial_state$dimension()) proposal$update(shape = shape) } - update <- function(sample_index, state_and_statistics) { + update <- function(proposal, sample_index, state_and_statistics) { + if (is.null(target_accept_prob)) { + target_accept_prob <- proposal$default_target_accept_prob() + } momentum <- state_and_statistics$proposed_state$momentum() accept_prob <- state_and_statistics$statistics$accept_prob shape <<- ramcmc::adapt_S( diff --git a/R/chains.R b/R/chains.R index b164924..7abbf80 100644 --- a/R/chains.R +++ b/R/chains.R @@ -163,7 +163,7 @@ chain_loop <- function( statistic_names) { progress_bar <- get_progress_bar(use_progress_bar, n_iteration, stage_name) for (adapter in adapters) { - adapter$initialize(state) + adapter$initialize(proposal, state) } if (record_traces_and_statistics) { trace_names <- names(unlist(trace_function(state))) @@ -176,25 +176,25 @@ chain_loop <- function( traces <- NULL statistics <- NULL } - for (s in seq_len(n_iteration)) { + for (chain_iteration in seq_len(n_iteration)) { state_and_statistics <- sample_metropolis_hastings( state, target_distribution, proposal ) for (adapter in adapters) { - adapter$update(s + 1, state_and_statistics) + adapter$update(proposal, chain_iteration + 1, state_and_statistics) } state <- state_and_statistics$state if (record_traces_and_statistics) { - traces[s, ] <- unlist(trace_function(state)) + traces[chain_iteration, ] <- unlist(trace_function(state)) adapter_states <- lapply(adapters, function(a) a$state()) - statistics[s, ] <- unlist( + statistics[chain_iteration, ] <- unlist( c(state_and_statistics$statistics, adapter_states) ) } if (!is.null(progress_bar)) progress_bar$tick() } for (adapter in adapters) { - if (!is.null(adapter$finalize)) adapter$finalize() + if (!is.null(adapter$finalize)) adapter$finalize(proposal) } list(final_state = state, traces = traces, statistics = statistics) } diff --git a/README.Rmd b/README.Rmd index 424d19a..ad1f428 100644 --- a/README.Rmd +++ b/README.Rmd @@ -59,7 +59,7 @@ results <- sample_chain( initial_state = rnorm(dimension), n_warm_up_iteration = 1000, n_main_iteration = 1000, - adapters = list(scale_adapter(proposal), variance_adapter(proposal)) + adapters = list(scale_adapter(), variance_adapter()) ) mean_accept_prob <- mean(results$statistics[, "accept_prob"]) adapted_shape <- proposal$parameters()$shape diff --git a/README.md b/README.md index 6c07526..6352590 100644 --- a/README.md +++ b/README.md @@ -56,7 +56,7 @@ results <- sample_chain( initial_state = rnorm(dimension), n_warm_up_iteration = 1000, n_main_iteration = 1000, - adapters = list(scale_adapter(proposal), variance_adapter(proposal)) + adapters = list(scale_adapter(), variance_adapter()) ) mean_accept_prob <- mean(results$statistics[, "accept_prob"]) adapted_shape <- proposal$parameters()$shape diff --git a/man/robust_shape_adapter.Rd b/man/robust_shape_adapter.Rd index 9f8cac4..ea0fa9e 100644 --- a/man/robust_shape_adapter.Rd +++ b/man/robust_shape_adapter.Rd @@ -6,16 +6,12 @@ Metropolis algorithm of Vihola (2012).} \usage{ robust_shape_adapter( - proposal, - initial_scale, - target_accept_prob = 0.4, + initial_scale = NULL, + target_accept_prob = NULL, kappa = 0.6 ) } \arguments{ -\item{proposal}{Proposal object to adapt. Must define an \code{update} function -which accepts a parameter \code{scale} for setting scale parameter of proposal.} - \item{initial_scale}{Initial value to use for scale parameter. If not set explicitly a proposal and dimension dependent default will be used.} @@ -27,8 +23,8 @@ chain. If not set a proposal dependent default will be used.} \value{ List of functions with entries \itemize{ -\item \code{initialize}, a function for initializing adapter state at beginning of -chain, +\item \code{initialize}, a function for initializing adapter state and proposal +parameters at beginning of chain, \item \code{update} a function for updating adapter state and proposal parameters on each chain iteration, \item \code{finalize} a function for performing any final updates to adapter state and @@ -47,11 +43,8 @@ target_distribution <- list( grad_log_density = function(x) -x ) proposal <- barker_proposal(target_distribution) -adapter <- robust_shape_adapter( - proposal, - initial_scale = 1., - target_accept_prob = 0.4 -) +adapter <- robust_shape_adapter(initial_scale = 1., target_accept_prob = 0.4) +adapter$initialize(proposal, chain_state(c(0, 0))) } \references{ Vihola, M. (2012). Robust adaptive Metropolis algorithm with diff --git a/man/scale_adapter.Rd b/man/scale_adapter.Rd index 9943a85..76fc84f 100644 --- a/man/scale_adapter.Rd +++ b/man/scale_adapter.Rd @@ -4,17 +4,9 @@ \alias{scale_adapter} \title{Create object to adapt proposal scale to coerce average acceptance rate.} \usage{ -scale_adapter( - proposal, - initial_scale = NULL, - target_accept_prob = NULL, - kappa = 0.6 -) +scale_adapter(initial_scale = NULL, target_accept_prob = NULL, kappa = 0.6) } \arguments{ -\item{proposal}{Proposal object to adapt. Must define an \code{update} function -which accepts a parameter \code{scale} for setting scale parameter of proposal.} - \item{initial_scale}{Initial value to use for scale parameter. If not set explicitly a proposal and dimension dependent default will be used.} @@ -22,12 +14,15 @@ explicitly a proposal and dimension dependent default will be used.} chain. If not set a proposal dependent default will be used.} \item{kappa}{Decay rate exponent in \verb{[0.5, 1]} for adaptation learning rate.} + +\item{proposal}{Proposal object to adapt. Must define an \code{update} function +which accepts a parameter \code{scale} for setting scale parameter of proposal.} } \value{ List of functions with entries \itemize{ -\item \code{initialize}, a function for initializing adapter state at beginning of -chain, +\item \code{initialize}, a function for initializing adapter state and proposal +parameters at beginning of chain, \item \code{update} a function for updating adapter state and proposal parameters on each chain iteration, \item \code{finalize} a function for performing any final updates to adapter state and @@ -46,8 +41,6 @@ target_distribution <- list( grad_log_density = function(x) -x ) proposal <- barker_proposal(target_distribution) -adapter <- scale_adapter( - proposal, - initial_scale = 1., target_accept_prob = 0.4 -) +adapter <- scale_adapter(initial_scale = 1., target_accept_prob = 0.4) +adapter$initialize(proposal, chain_state(c(0, 0))) } diff --git a/man/variance_adapter.Rd b/man/variance_adapter.Rd index b6cee22..fadc395 100644 --- a/man/variance_adapter.Rd +++ b/man/variance_adapter.Rd @@ -5,19 +5,19 @@ \title{Create object to adapt proposal with per dimension scales based on estimates of target distribution variances.} \usage{ -variance_adapter(proposal, kappa = 0.6) +variance_adapter(kappa = 0.6) } \arguments{ +\item{kappa}{Decay rate exponent in \verb{[0.5, 1]} for adaptation learning rate.} + \item{proposal}{Proposal object to adapt. Must define an \code{update} function which accepts a parameter \code{shape} for setting shape parameter of proposal.} - -\item{kappa}{Decay rate exponent in \verb{[0.5, 1]} for adaptation learning rate.} } \value{ List of functions with entries \itemize{ -\item \code{initialize}, a function for initializing adapter state at beginning of -chain, +\item \code{initialize}, a function for initializing adapter state and proposal +parameters at beginning of chain, \item \code{update} a function for updating adapter state and proposal parameters on each chain iteration, \item \code{finalize} a function for performing any final updates to adapter state and @@ -37,5 +37,6 @@ target_distribution <- list( grad_log_density = function(x) -x ) proposal <- barker_proposal(target_distribution) -adapter <- variance_adapter(proposal) +adapter <- variance_adapter() +adapter$initialize(proposal, chain_state(c(0, 0))) } diff --git a/tests/testthat/test-adaptation.R b/tests/testthat/test-adaptation.R index 9eed6c8..a689db5 100644 --- a/tests/testthat/test-adaptation.R +++ b/tests/testthat/test-adaptation.R @@ -24,7 +24,8 @@ dummy_proposal_with_scale_parameter <- function(scale = NULL) { dummy_proposal_with_shape_parameter <- function(shape = NULL) { list( update = function(shape) shape <<- shape, - parameters = function() list(shape = shape) + parameters = function() list(shape = shape), + default_initial_scale = function(dimension) 1 / sqrt(dimension) ) } @@ -40,23 +41,22 @@ for (target_accept_prob in c(0.2, 0.4, 0.6)) { target_accept_prob, initial_scale, kappa ), { - proposal <- dummy_proposal_with_scale_parameter() adapter <- scale_adapter( - proposal = proposal, initial_scale = initial_scale, target_accept_prob = target_accept_prob, kappa = kappa ) check_adapter(adapter) - adapter$initialize(initial_state = chain_state(rep(0, dimension))) + proposal <- dummy_proposal_with_scale_parameter() + adapter$initialize(proposal, chain_state(rep(0, dimension))) adapter_state <- adapter$state() expect_named(adapter_state, "log_scale") expect_length(adapter_state$log_scale, 1) old_scale <- initial_scale # If accept probability higher than target scale should be increased - for (s in 1:2) { + for (sample_index in 1:2) { adapter$update( - s, list(statistics = list(accept_prob = target_accept_prob + 0.1)) + proposal, sample_index, list(statistics = list(accept_prob = target_accept_prob + 0.1)) ) expect_type(adapter$state(), "list") scale <- proposal$parameters()$scale @@ -64,9 +64,9 @@ for (target_accept_prob in c(0.2, 0.4, 0.6)) { old_scale <- scale } # If accept probability lower than target scale should be decreased - for (s in 3:4) { + for (sample_index in 3:4) { adapter$update( - s, list(statistics = list(accept_prob = target_accept_prob - 0.1)) + proposal, sample_index, list(statistics = list(accept_prob = target_accept_prob - 0.1)) ) scale <- proposal$parameters()$scale expect_lt(scale, old_scale) @@ -75,10 +75,10 @@ for (target_accept_prob in c(0.2, 0.4, 0.6)) { # For a smooth decreasing relation between accept probability and # scale should adapt over long run to give close to target accept # probability - for (s in 5:2000) { + for (sample_index in 5:2000) { accept_prob <- exp(-scale) adapter$update( - s, list(statistics = list(accept_prob = accept_prob)) + proposal, sample_index, list(statistics = list(accept_prob = accept_prob)) ) scale <- proposal$parameters()$scale } @@ -96,10 +96,10 @@ for (dimension in c(1L, 2L, 5L)) { dimension ), { - proposal <- dummy_proposal_with_scale_parameter() - adapter <- scale_adapter(proposal) + adapter <- scale_adapter() check_adapter(adapter) - adapter$initialize(initial_state = chain_state(rep(0, dimension))) + proposal <- dummy_proposal_with_scale_parameter() + adapter$initialize(proposal, chain_state(rep(0, dimension))) adapter_state <- adapter$state() expect_named(adapter_state, "log_scale") expect_length(adapter_state$log_scale, 1) @@ -121,13 +121,13 @@ for (dimension in c(1L, 2L, 5L)) { ), { proposal <- dummy_proposal_with_shape_parameter() - adapter <- variance_adapter(proposal = proposal, kappa = kappa) + adapter <- variance_adapter(kappa = kappa) check_adapter(adapter) withr::local_seed(default_seed()) target_scales <- exp(2 * rnorm(dimension)) position <- rnorm(dimension) * target_scales state <- chain_state(position = position) - adapter$initialize(state) + adapter$initialize(proposal, state) adapter_state <- adapter$state() expect_named( adapter_state, @@ -146,7 +146,7 @@ for (dimension in c(1L, 2L, 5L)) { + sqrt(1 - correlation^2) * target_scales * rnorm(dimension) ) state$update(position = position) - adapter$update(s, list(state = state)) + adapter$update(proposal, s, list(state = state)) } expect_equal( proposal$parameters()$shape, @@ -178,14 +178,13 @@ for (dimension in c(1L, 2L, 3L)) { ) proposal <- random_walk_proposal(target_distribution) adapter <- robust_shape_adapter( - proposal = proposal, initial_scale = 1., kappa = 2 / 3, target_accept_prob = target_accept_prob ) check_adapter(adapter) state <- chain_state(position = rnorm(dimension)) - adapter$initialize(state) + adapter$initialize(proposal, state) adapter_state <- adapter$state() expect_named(adapter_state, "shape") expect_nrow(adapter_state$shape, dimension) @@ -195,7 +194,7 @@ for (dimension in c(1L, 2L, 3L)) { state_and_statistics <- sample_metropolis_hastings( state, target_distribution, proposal ) - adapter$update(sample_index, state_and_statistics) + adapter$update(proposal, sample_index, state_and_statistics) state <- state_and_statistics$state mean_accept_prob <- mean_accept_prob + ( state_and_statistics$statistics$accept_prob - mean_accept_prob @@ -215,3 +214,23 @@ for (dimension in c(1L, 2L, 3L)) { ) } } + +for (dimension in c(1L, 2L, 5L)) { + test_that( + sprintf( + "Robust shape adapter with only proposal specified works in dimension %i", + dimension + ), + { + adapter <- robust_shape_adapter() + check_adapter(adapter) + proposal <- dummy_proposal_with_shape_parameter() + adapter$initialize(proposal, chain_state(rep(0, dimension))) + adapter_state <- adapter$state() + expect_named(adapter_state, "shape") + expect_nrow(adapter_state$shape, dimension) + expect_ncol(adapter_state$shape, dimension) + expect_equal(adapter_state$shape, diag(dimension) / sqrt(dimension)) + } + ) +} diff --git a/tests/testthat/test-chains.R b/tests/testthat/test-chains.R index 1e1cebc..727f8fc 100644 --- a/tests/testthat/test-chains.R +++ b/tests/testthat/test-chains.R @@ -23,9 +23,7 @@ for (n_warm_up_iteration in c(0, 1, 10)) { target_distribution <- standard_normal_target_distribution() barker_proposal(target_distribution) proposal <- barker_proposal(target_distribution) - adapters <- list( - scale_adapter(proposal, initial_scale = 1.) - ) + adapters <- list(scale_adapter(initial_scale = 1.)) withr::with_seed(default_seed(), { position <- rnorm(dimension) })