Provide per proposal defaults for target acceptance rate and initial …

…scale (#42)

* Provide proposal dependent defaults for target acceptance rate and initial scale

* Update test to cover case using default initial_scale

R/adaptation.R

@@ -2,9 +2,10 @@
 #'
 #' @param proposal Proposal object to adapt. Must define an `update` function
 #'   which accepts a parameter `scale` for setting scale parameter of proposal.
-#' @param initial_scale Initial value to use for scale parameter.
+#' @param initial_scale Initial value to use for scale parameter. If not set
+#'   explicitly a proposal and dimension dependent default will be used.
 #' @param target_accept_prob Target value for average accept probability for
-#'   chain.
+#'   chain. If not set a proposal dependent default will be used.
 #' @param kappa Decay rate exponent in `[0.5, 1]` for adaptation learning rate.
 #'
 #' @return List of functions with entries
@@ -31,9 +32,15 @@
 #'   initial_scale = 1., target_accept_prob = 0.4
 #' )
 scale_adapter <- function(
-    proposal, initial_scale, target_accept_prob = 0.4, kappa = 0.6) {
+    proposal, 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) {
+    if (is.null(initial_scale)) {
+      initial_scale <- proposal$default_initial_scale(initial_state$dimension())
+    }
     log_scale <<- log(initial_scale)
     proposal$update(scale = initial_scale)
   }

R/barker.R

@@ -87,6 +87,10 @@ log_density_ratio_barker <- function(
 #'   for a given pair of current and proposed chain states,
 #' * `update`: a function to update parameters of proposal,
 #' * `parameters`: a function to return list of current parameter values.
+#' * `default_target_accept_prob`: a function returning the default target
+#'   acceptance rate to use for any scale adaptation.
+#' * `default_initial_scale`: a function which given a dimension gives a default
+#'   value to use for the initial proposal scale parameter.
 #'
 #' @export
 #'
@@ -118,6 +122,8 @@ barker_proposal <- function(
       )
     },
     scale = scale,
-    shape = shape
+    shape = shape,
+    default_target_accept_prob = 0.4,
+    default_initial_scale = function(dimension) 2.38 / (dimension)^(1 / 3)
   )
 }

R/langevin.R

@@ -80,6 +80,8 @@ langevin_proposal <- function(
       )
     },
     scale = scale,
-    shape = shape
+    shape = shape,
+    default_target_accept_prob = 0.574,
+    default_initial_scale = function(dimension) 2.38 / (dimension)^(1 / 3)
   )
 }

R/proposal.R

@@ -20,7 +20,12 @@ get_shape_matrix <- function(scale, shape) {
   }
 }
 
-scale_and_shape_proposal <- function(sample, log_density_ratio, scale, shape) {
+scale_and_shape_proposal <- function(
+    sample,
+    log_density_ratio,
+    scale, shape,
+    default_target_accept_prob,
+    default_initial_scale) {
   scale <- scale
   shape <- shape
   list(
@@ -38,6 +43,8 @@ scale_and_shape_proposal <- function(sample, log_density_ratio, scale, shape) {
         shape <<- shape
       }
     },
-    parameters = function() list(scale = scale, shape = shape)
+    parameters = function() list(scale = scale, shape = shape),
+    default_target_accept_prob = function() default_target_accept_prob,
+    default_initial_scale = default_initial_scale
   )
 }

R/random_walk.R

@@ -42,6 +42,8 @@ random_walk_proposal <- function(
     },
     log_density_ratio = function(state, proposed_state, scale_and_shape) 0,
     scale = scale,
-    shape = shape
+    shape = shape,
+    default_target_accept_prob = 0.234,
+    default_initial_scale = function(dimension) 2.38 / (dimension)^(1 / 2)
   )
 }

README.Rmd

@@ -47,26 +47,19 @@ library(rmcmc)
 
 set.seed(876287L)
 dimension <- 3
-scales <- exp(2 * rnorm(dimension))
+scales <- exp(rnorm(dimension))
 target_distribution <- list(
   log_density = function(x) -sum((x / scales)^2) / 2,
   gradient_log_density = function(x) -x / scales^2
 )
 proposal <- barker_proposal(target_distribution)
-adapters <- list(
-  scale_adapter(proposal, initial_scale = 1., target_accept_prob = 0.4),
-  variance_adapter(proposal)
-)
-n_warm_up_iteration <- 1000
-n_main_iteration <- 1000
-initial_state <- chain_state(rnorm(dimension))
 results <- sample_chain(
   target_distribution = target_distribution,
   proposal = proposal,
-  initial_state = initial_state,
-  n_warm_up_iteration = n_warm_up_iteration,
-  n_main_iteration = n_main_iteration,
-  adapters = adapters
+  initial_state = rnorm(dimension),
+  n_warm_up_iteration = 1000,
+  n_main_iteration = 1000,
+  adapters = list(scale_adapter(proposal), variance_adapter(proposal))
 )
 mean_accept_prob <- mean(results$statistics[, "accept_prob"])
 adapted_shape <- proposal$parameters()$shape

README.md

@@ -44,26 +44,19 @@ library(rmcmc)
 
 set.seed(876287L)
 dimension <- 3
-scales <- exp(2 * rnorm(dimension))
+scales <- exp(rnorm(dimension))
 target_distribution <- list(
   log_density = function(x) -sum((x / scales)^2) / 2,
   gradient_log_density = function(x) -x / scales^2
 )
 proposal <- barker_proposal(target_distribution)
-adapters <- list(
-  scale_adapter(proposal, initial_scale = 1., target_accept_prob = 0.4),
-  variance_adapter(proposal)
-)
-n_warm_up_iteration <- 1000
-n_main_iteration <- 1000
-initial_state <- chain_state(rnorm(dimension))
 results <- sample_chain(
   target_distribution = target_distribution,
   proposal = proposal,
-  initial_state = initial_state,
-  n_warm_up_iteration = n_warm_up_iteration,
-  n_main_iteration = n_main_iteration,
-  adapters = adapters
+  initial_state = rnorm(dimension),
+  n_warm_up_iteration = 1000,
+  n_main_iteration = 1000,
+  adapters = list(scale_adapter(proposal), variance_adapter(proposal))
 )
 mean_accept_prob <- mean(results$statistics[, "accept_prob"])
 adapted_shape <- proposal$parameters()$shape
@@ -73,7 +66,7 @@ cat(
   sprintf("Adapter scale est.: %s", toString(adapted_shape)),
   sep = "\n"
 )
-#> Average acceptance probability is 0.46
-#> True target scales: 2.26617033226883, 1.89818769776724, 0.0767505506297473
-#> Adapter scale est.: 1.77277384748788, 1.71554065105575, 0.0804144979270686
+#> Average acceptance probability is 0.40
+#> True target scales: 1.50538046096953, 1.37774732725824, 0.277038897322645
+#> Adapter scale est.: 1.35010920408606, 1.5140138215658, 0.248974800274054
 ```

tests/testthat/test-adaptation.R

@@ -15,7 +15,9 @@ check_adapter <- function(adapter) {
 dummy_proposal_with_scale_parameter <- function(scale = NULL) {
   list(
     update = function(scale) scale <<- scale,
-    parameters = function() list(scale = scale)
+    parameters = function() list(scale = scale),
+    default_target_accept_prob = function() 0.234,
+    default_initial_scale = function(dimension) 1 / sqrt(dimension)
   )
 }
 
@@ -46,7 +48,7 @@ for (target_accept_prob in c(0.2, 0.4, 0.6)) {
             kappa = kappa
           )
           check_adapter(adapter)
-          adapter$initialize(initial_state = NULL)
+          adapter$initialize(initial_state = chain_state(rep(0, dimension)))
           adapter_state <- adapter$state()
           expect_named(adapter_state, "log_scale")
           expect_length(adapter_state$log_scale, 1)
@@ -87,6 +89,25 @@ for (target_accept_prob in c(0.2, 0.4, 0.6)) {
   }
 }
 
+for (dimension in c(1L, 2L, 5L)) {
+  test_that(
+    sprintf(
+      "Scale adapter with only proposal specified works in dimension %i",
+      dimension
+    ),
+    {
+      proposal <- dummy_proposal_with_scale_parameter()
+      adapter <- scale_adapter(proposal)
+      check_adapter(adapter)
+      adapter$initialize(initial_state = chain_state(rep(0, dimension)))
+      adapter_state <- adapter$state()
+      expect_named(adapter_state, "log_scale")
+      expect_length(adapter_state$log_scale, 1)
+      expect_equal(adapter_state$log_scale, -log(dimension) / 2)
+    }
+  )
+}
+
 for (dimension in c(1L, 2L, 5L)) {
   for (kappa in c(0.5, 0.6, 0.8)) {
     for (correlation in c(0.0, 0.5, 0.8)) {

tests/testthat/test-proposal.R

@@ -28,7 +28,9 @@ for (dimension in c(1, 2)) {
       chain_state(position = state$position() + offset)
     }
     log_density_ratio <- function(state, proposed_state, scale_and_shape) 0
-    proposal <- scale_and_shape_proposal(sample, log_density_ratio, NULL, NULL)
+    proposal <- scale_and_shape_proposal(
+      sample, log_density_ratio, NULL, NULL, NULL, function(d) NULL
+    )
     state <- chain_state(position = rnorm(dimension))
     expect_error(proposal$sample(state), "must be set")
     expect_identical(proposal$parameters(), list(scale = NULL, shape = NULL))
@@ -40,5 +42,7 @@ for (dimension in c(1, 2)) {
     check_chain_state(proposed_state)
     expect_all_different(state$position(), proposed_state$position())
     expect_identical(proposal$log_density(state, proposed_state), 0)
+    expect_identical(proposal$default_target_accept_prob(), NULL)
+    expect_identical(proposal$default_initial_scale(dimension), NULL)
   })
 }