Add vignette demoing adjusting Barker proposal noise distribution (#67)

* Add vignette demoing changing Barker proposal noise distribution

* Decouple descriptions of alt proposals from barker_proposal

* Add bimodal Barker proposal convenience function

* Fix linter errors

Author: matt-graham

NAMESPACE

@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 export(barker_proposal)
+export(bimodal_barker_proposal)
 export(chain_state)
 export(covariance_shape_adapter)
 export(dual_averaging_scale_adapter)

R/barker.R

@@ -64,11 +64,13 @@ log_density_ratio_barker <- function(
 
 #' Create a new Barker proposal object.
 #'
-#' Returns a list with function to sample from the proposal, evaluate the log
-#' density ratio for a state pair for the proposal and update the proposal
-#' parameters. The proposal has two parameters `scale` and `shape`. At least one
-#' of `scale` and `shape` must be set before sampling from the proposal or
-#' evaluating the log density ratio.
+#' The Barker proposal is a gradient-based proposal inspired by the Barker
+#' accept-reject rule and proposed in Livingstone and Zanella (2022). It offers
+#' improved robustness compared to alternative gradient-based proposals such as
+#' Langevin proposals.
+#'
+#' For more details see the vignette:
+#' \code{vignette("barker-proposal", package = "rmcmc")}
 #'
 #' @inheritParams sample_barker
 #' @param scale Scale parameter of proposal distribution. A non-negative scalar
@@ -89,6 +91,11 @@ log_density_ratio_barker <- function(
 #' * `default_initial_scale`: a function which given a dimension gives a default
 #'   value to use for the initial proposal scale parameter.
 #'
+#' @references Livingstone, S., & Zanella, G. (2022). The Barker proposal:
+#'   combining robustness and efficiency in gradient-based MCMC. _Journal of the
+#'   Royal Statistical Society Series B: Statistical Methodology_, 84(2),
+#'   496-523.
+#'
 #' @export
 #'
 #' @examples
@@ -123,3 +130,63 @@ barker_proposal <- function(
     default_initial_scale = function(dimension) dimension^(-1 / 6)
   )
 }
+
+
+#' Create a new Barker proposal object with bimodal noise distribution.
+#'
+#' Convenience function for creating a Barker proposal with bimodal auxiliary
+#' noise variable distribution, corresponding to equally-weighted normal
+#' components with shared variance `sigma` and means `±sqrt(1 - sigma^2)`.
+#' This choice of noise distribution was suggested in Vogrinc et al. (2023) and
+#' found to give improved performance over the default choice of a standard
+#' normal auxiliary noise distribution in a range of targets.
+#'
+#' For more details see the vignette:
+#' \code{vignette("adjusting-noise-distribution", package = "rmcmc")}
+#'
+#' @inherit barker_proposal params return
+#'
+#' @param sigma Standard deviation of equally-weighted normal components in
+#'   bimodal auxiliary noise distribution, with corresponding means of
+#'   `±sqrt(1 - sigma^2)`.
+#'
+#' @references Vogrinc, J., Livingstone, S., & Zanella, G. (2023). Optimal
+#'   design of the Barker proposal and other locally balanced
+#'   Metropolis–Hastings algorithms. _Biometrika_, 110(3), 579-595.
+#'
+#' @seealso [barker_proposal()]
+#'
+#' @export
+#'
+#' @examples
+#' target_distribution <- list(
+#'   log_density = function(x) -sum(x^2) / 2,
+#'   gradient_log_density = function(x) -x
+#' )
+#' proposal <- bimodal_barker_proposal(scale = 1.)
+#' state <- chain_state(c(0., 0.))
+#' withr::with_seed(
+#'   876287L, proposed_state <- proposal$sample(state, target_distribution)
+#' )
+#' log_density_ratio <- proposal$log_density_ratio(
+#'   state, proposed_state, target_distribution
+#' )
+#' proposal$update(scale = 0.5)
+bimodal_barker_proposal <- function(
+    sigma = 0.1,
+    scale = NULL,
+    shape = NULL,
+    sample_uniform = stats::runif) {
+  sample_bimodal <- function(dimension) {
+    return(
+      sample(c(-1, 1), dimension, TRUE) * sqrt(1 - sigma^2) +
+        stats::rnorm(dimension) * sigma
+    )
+  }
+  barker_proposal(
+    scale = scale,
+    shape = shape,
+    sample_auxiliary = sample_bimodal,
+    sample_uniform = sample_uniform
+  )
+}

R/hamiltonian.R

@@ -59,7 +59,14 @@ log_density_ratio_hamiltonian <- function(
 
 #' Create a new Hamiltonian proposal object.
 #'
-#' @inherit barker_proposal return params description
+#' The Hamiltonian proposal augments the target distribution with normally
+#' distributed auxiliary momenta variables and simulates the dynamics for a
+#' Hamiltonian function corresponding to the negative logarithm of the density
+#' of the resulting joint target distribution using a leapfrog integrator, with
+#' the proposed new state being the forward integrate state with momenta negated
+#' to ensure reversibility.
+#'
+#' @inherit barker_proposal return params
 #'
 #' @param n_step Number of leapfrog steps to simulate Hamiltonian dynamics for
 #'   in each proposed move, or parameter passed to function specified by
@@ -79,6 +86,11 @@ log_density_ratio_hamiltonian <- function(
 #'   be used to specify parameter(s) of distribution to sample number of steps
 #'   from.
 #'
+#' @references Duane, S., Kennedy, A. D., Pendleton, B. J., & Roweth, D. (1987).
+#'   Hybrid Monte Carlo. _Physics Letters B_, 195(2), 216-222.
+#' @references Neal, R. M. (2011). MCMC Using Hamiltonian Dynamics. In _Handbook
+#'   of Markov Chain Monte Carlo_ (pp. 113-162). Chapman and Hall/CRC.
+#'
 #' @export
 #'
 #' @examples

R/langevin.R

@@ -45,7 +45,17 @@ log_density_ratio_langevin <- function(
 
 #' Create a new Langevin proposal object.
 #'
-#' @inherit barker_proposal return params description
+#' The Langevin proposal is a gradient-based proposal corresponding to a
+#' Euler-Maruyama time discretisation of a Langevin diffusion.
+#'
+#' @inherit barker_proposal return params
+#'
+#' @references Besag, J. (1994). "Comments on "Representations of knowledge in
+#'   complex systems" by U. Grenander and MI Miller". _Journal of the Royal
+#'   Statistical Society, Series B_. 56: 591–592.
+#' @references Roberts, G. O., & Tweedie, R. L. (1996). Exponential convergence
+#'   of Langevin distributions and their discrete approximations. _Bernoulli_ 2
+#'   (4), 341 - 363.
 #'
 #' @export
 #'

R/random_walk.R

@@ -26,9 +26,12 @@ log_density_ratio_random_walk <- function(
   0
 }
 
-#' Create a new random walk proposal object.
+#' Create a new (Gaussian) random walk proposal object.
 #'
-#' @inherit barker_proposal return params description
+#' The Gaussian random walk proposal samples a new proposed state by perturbing
+#' the current state with zero-mean normally distributed noise.
+#'
+#' @inherit barker_proposal return params
 #'
 #' @export
 #'