Merge pull request #65 from UCL/mmg/default-proposal

R/adaptation.R

@@ -31,11 +31,7 @@
 #' @export
 #'
 #' @examples
-#' target_distribution <- list(
-#'   log_density = function(x) -sum(x^2) / 2,
-#'   grad_log_density = function(x) -x
-#' )
-#' proposal <- barker_proposal(target_distribution)
+#' proposal <- barker_proposal()
 #' adapter <- scale_adapter(initial_scale = 1., target_accept_prob = 0.4)
 #' adapter$initialize(proposal, chain_state(c(0, 0)))
 scale_adapter <- function(
@@ -69,11 +65,7 @@ scale_adapter <- function(
 #' @export
 #'
 #' @examples
-#' target_distribution <- list(
-#'   log_density = function(x) -sum(x^2) / 2,
-#'   grad_log_density = function(x) -x
-#' )
-#' proposal <- barker_proposal(target_distribution)
+#' proposal <- barker_proposal()
 #' adapter <- stochastic_approximation_scale_adapter(
 #'   initial_scale = 1., target_accept_prob = 0.4
 #' )
@@ -134,11 +126,7 @@ stochastic_approximation_scale_adapter <- function(
 #' @export
 #'
 #' @examples
-#' target_distribution <- list(
-#'   log_density = function(x) -sum(x^2) / 2,
-#'   grad_log_density = function(x) -x
-#' )
-#' proposal <- barker_proposal(target_distribution)
+#' proposal <- barker_proposal()
 #' adapter <- dual_averaging_scale_adapter(
 #'   initial_scale = 1., target_accept_prob = 0.4
 #' )
@@ -212,11 +200,7 @@ dual_averaging_scale_adapter <- function(
 #'
 #' @export
 #' @examples
-#' target_distribution <- list(
-#'   log_density = function(x) -sum(x^2) / 2,
-#'   grad_log_density = function(x) -x
-#' )
-#' proposal <- barker_proposal(target_distribution)
+#' proposal <- barker_proposal()
 #' adapter <- shape_adapter()
 #' adapter$initialize(proposal, chain_state(c(0, 0)))
 shape_adapter <- function(type = "covariance", kappa = 1) {
@@ -247,11 +231,7 @@ shape_adapter <- function(type = "covariance", kappa = 1) {
 #'
 #' @export
 #' @examples
-#' target_distribution <- list(
-#'   log_density = function(x) -sum(x^2) / 2,
-#'   grad_log_density = function(x) -x
-#' )
-#' proposal <- barker_proposal(target_distribution)
+#' proposal <- barker_proposal()
 #' adapter <- variance_shape_adapter()
 #' adapter$initialize(proposal, chain_state(c(0, 0)))
 variance_shape_adapter <- function(kappa = 1) {
@@ -305,11 +285,7 @@ variance_shape_adapter <- function(kappa = 1) {
 #'
 #' @export
 #' @examples
-#' target_distribution <- list(
-#'   log_density = function(x) -sum(x^2) / 2,
-#'   grad_log_density = function(x) -x
-#' )
-#' proposal <- barker_proposal(target_distribution)
+#' proposal <- barker_proposal()
 #' adapter <- covariance_shape_adapter()
 #' adapter$initialize(proposal, chain_state(c(0, 0)))
 covariance_shape_adapter <- function(kappa = 1) {
@@ -361,11 +337,7 @@ covariance_shape_adapter <- function(kappa = 1) {
 #' @export
 #'
 #' @examples
-#' target_distribution <- list(
-#'   log_density = function(x) -sum(x^2) / 2,
-#'   grad_log_density = function(x) -x
-#' )
-#' proposal <- barker_proposal(target_distribution)
+#' proposal <- barker_proposal()
 #' adapter <- robust_shape_adapter(initial_scale = 1., target_accept_prob = 0.4)
 #' adapter$initialize(proposal, chain_state(c(0, 0)))
 robust_shape_adapter <- function(

R/barker.R

@@ -96,28 +96,27 @@ log_density_ratio_barker <- function(
 #'   log_density = function(x) -sum(x^2) / 2,
 #'   gradient_log_density = function(x) -x
 #' )
-#' proposal <- barker_proposal(target_distribution, scale = 1.)
+#' proposal <- barker_proposal(scale = 1.)
 #' state <- chain_state(c(0., 0.))
-#' withr::with_seed(876287L, proposed_state <- proposal$sample(state))
-#' log_density_ratio <- proposal$log_density_ratio(state, proposed_state)
+#' 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)
 barker_proposal <- function(
-    target_distribution,
     scale = NULL,
     shape = NULL,
     sample_auxiliary = stats::rnorm,
     sample_uniform = stats::runif) {
   scale_and_shape_proposal(
-    sample = function(state, scale_and_shape) {
+    sample = function(state, target_distribution, scale_and_shape) {
       sample_barker(
         state, target_distribution, scale_and_shape, sample_auxiliary, sample_uniform
       )
     },
-    log_density_ratio = function(state, proposed_state, scale_and_shape) {
-      log_density_ratio_barker(
-        state, proposed_state, target_distribution, scale_and_shape
-      )
-    },
+    log_density_ratio = log_density_ratio_barker,
     scale = scale,
     shape = shape,
     default_target_accept_prob = 0.574,

R/chains.R

@@ -1,8 +1,8 @@
 #' Sample a Markov chain
 #'
-#' Sample a Markov chain using Metropolis-Hastings kernel with given proposal
-#' and target distributions, optionally adapting proposal parameters in warm-up
-#' stage.
+#' Sample a Markov chain using Metropolis-Hastings kernel with a user-specified
+#' target distribution and proposal (defaulting to Barker proposal), optionally
+#' adapting proposal parameters in a warm-up stage.
 #'
 #' @inheritParams sample_metropolis_hastings
 #' @param initial_state Initial chain state. Either a vector specifying just
@@ -12,6 +12,20 @@
 #'   run.
 #' @param n_main_iteration Number of main (non-adaptive) chain iterations to
 #'   run.
+#' @param proposal Proposal distribution object. Defaults to Barker proposal,
+#'   that is the output of [barker_proposal()]. Proposal objects are lists which
+#'   must minimally define entries `sample`, a function to generate sample from
+#'   proposal distribution given current chain state and `log_density_ratio`, a
+#'   function to compute log density ratio for proposal for a given pair of
+#'   current and proposed chain states. If adapters are being used to adaptively
+#'   tune the proposal scale and shape parameters, which is the default
+#'   behaviour of `sample_chain`, then additionally the list must also define
+#'   entries: `update` a function for updating parameters of proposal,
+#'   `parameters` a function for getting current proposal parameter values,
+#'   `default_target_accept_prob` a function for getting proposal specific
+#'   default target acceptance probability for scale adaptation and
+#'   `default_initial_scale` a function for getting proposal and dimension
+#'   dependent default initial value for scale parameter.
 #' @param adapters List of adapters to tune proposal parameters during warm-up.
 #'   Defaults to using list with instances of [scale_adapter()] and
 #'   [shape_adapter()], corresponding to respectively, adapting the scale to
@@ -49,25 +63,20 @@
 #'   log_density = function(x) -sum(x^2) / 2,
 #'   gradient_log_density = function(x) -x
 #' )
-#' proposal <- barker_proposal(target_distribution, scale = 1.)
-#' n_warm_up_iteration <- 1000
-#' n_main_iteration <- 1000
 #' withr::with_seed(876287L, {
-#'   initial_state <- chain_state(stats::rnorm(2))
 #'   results <- sample_chain(
 #'     target_distribution,
-#'     proposal,
-#'     initial_state,
-#'     n_warm_up_iteration,
-#'     n_main_iteration
+#'     initial_state = stats::rnorm(2),
+#'     n_warm_up_iteration = 1000,
+#'     n_main_iteration = 1000
 #'   )
 #' })
 sample_chain <- function(
     target_distribution,
-    proposal,
     initial_state,
     n_warm_up_iteration,
     n_main_iteration,
+    proposal = barker_proposal(),
     adapters = list(scale_adapter(), shape_adapter()),
     trace_function = NULL,
     show_progress_bar = TRUE,

R/hamiltonian.R

@@ -88,23 +88,28 @@ log_density_ratio_hamiltonian <- function(
 #' )
 #'
 #' # Proposal with fixed number of leapfrog steps
-#' proposal <- hamiltonian_proposal(target_distribution, scale = 1., n_step = 5)
+#' proposal <- hamiltonian_proposal(scale = 1., n_step = 5)
 #' state <- chain_state(c(0., 0.))
-#' withr::with_seed(876287L, proposed_state <- proposal$sample(state))
-#' log_density_ratio <- proposal$log_density_ratio(state, proposed_state)
+#' 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)
 #'
 #' # Proposal with number of steps randomly sampled uniformly from 5:10
 #' sample_uniform_int <- function(lower, upper) {
 #'   lower + sample.int(upper - lower + 1, 1) - 1
 #' }
 #' proposal <- hamiltonian_proposal(
-#'   target_distribution,
 #'   scale = 1.,
 #'   n_step = c(5, 10),
 #'   sample_n_step = function(n_step) sample_uniform_int(n_step[1], n_step[2])
 #' )
-#' withr::with_seed(876287L, proposed_state <- proposal$sample(state))
+#' withr::with_seed(
+#'   876287L, proposed_state <- proposal$sample(state, target_distribution)
+#' )
 #'
 #' # Proposal with partial momentum refreshment
 #' partial_momentum_update <- function(state, phi = pi / 4) {
@@ -116,23 +121,21 @@ log_density_ratio_hamiltonian <- function(
 #'   }
 #' }
 #' proposal <- hamiltonian_proposal(
-#'   target_distribution,
 #'   scale = 1.,
 #'   n_step = 1,
 #'   sample_auxiliary = partial_momentum_update
 #' )
-#' withr::with_seed(876287L, {
-#'   proposed_state <- proposal$sample(state)
-#' })
+#' withr::with_seed(
+#'   876287L, proposed_state <- proposal$sample(state, target_distribution)
+#' )
 hamiltonian_proposal <- function(
-    target_distribution,
     n_step,
     scale = NULL,
     shape = NULL,
     sample_auxiliary = function(state) stats::rnorm(state$dimension()),
     sample_n_step = NULL) {
   scale_and_shape_proposal(
-    sample = function(state, scale_and_shape) {
+    sample = function(state, target_distribution, scale_and_shape) {
       sample_hamiltonian(
         state,
         target_distribution,
@@ -142,11 +145,7 @@ hamiltonian_proposal <- function(
         sample_n_step
       )
     },
-    log_density_ratio = function(state, proposed_state, scale_and_shape) {
-      log_density_ratio_hamiltonian(
-        state, proposed_state, target_distribution, scale_and_shape
-      )
-    },
+    log_density_ratio = log_density_ratio_hamiltonian,
     scale = scale,
     shape = shape,
     default_target_accept_prob = 0.8,

R/kernels.R

@@ -30,11 +30,11 @@ sample_metropolis_hastings <- function(
     target_distribution,
     proposal,
     sample_uniform = stats::runif) {
-  proposed_state <- proposal$sample(state)
+  proposed_state <- proposal$sample(state, target_distribution)
   log_accept_ratio <- (
     proposed_state$log_density(target_distribution) -
       state$log_density(target_distribution) +
-      proposal$log_density_ratio(state, proposed_state)
+      proposal$log_density_ratio(state, proposed_state, target_distribution)
   )
   accept_prob <- if (is.nan(log_accept_ratio)) 0 else min_1_exp(log_accept_ratio)
   if (sample_uniform(1) < accept_prob) {

R/langevin.R

@@ -54,27 +54,26 @@ log_density_ratio_langevin <- function(
 #'   log_density = function(x) -sum(x^2) / 2,
 #'   gradient_log_density = function(x) -x
 #' )
-#' proposal <- langevin_proposal(target_distribution, scale = 1.)
+#' proposal <- langevin_proposal(scale = 1.)
 #' state <- chain_state(c(0., 0.))
-#' withr::with_seed(876287L, proposed_state <- proposal$sample(state))
-#' log_density_ratio <- proposal$log_density_ratio(state, proposed_state)
+#' 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)
 langevin_proposal <- function(
-    target_distribution,
     scale = NULL,
     shape = NULL,
     sample_auxiliary = stats::rnorm) {
   scale_and_shape_proposal(
-    sample = function(state, scale_and_shape) {
+    sample = function(state, target_distribution, scale_and_shape) {
       sample_langevin(
         state, target_distribution, scale_and_shape, sample_auxiliary
       )
     },
-    log_density_ratio = function(state, proposed_state, scale_and_shape) {
-      log_density_ratio_langevin(
-        state, proposed_state, target_distribution, scale_and_shape
-      )
-    },
+    log_density_ratio = log_density_ratio_langevin,
     scale = scale,
     shape = shape,
     default_target_accept_prob = 0.574,

R/proposal.R

@@ -22,10 +22,13 @@ scale_and_shape_proposal <- function(
   scale <- scale
   shape <- shape
   list(
-    sample = function(state) sample(state, get_shape_matrix(scale, shape)),
-    log_density_ratio = function(state, proposed_state) {
+    sample = function(state, target_distribution) {
+      sample(state, target_distribution, get_shape_matrix(scale, shape))
+    },
+    log_density_ratio = function(state, proposed_state, target_distribution) {
+      shape_matrix <- get_shape_matrix(scale, shape)
       log_density_ratio(
-        state, proposed_state, get_shape_matrix(scale, shape)
+        state, proposed_state, target_distribution, shape_matrix
       )
     },
     update = function(scale = NULL, shape = NULL) {

R/random_walk.R

@@ -13,6 +13,18 @@ sample_random_walk <- function(
   chain_state(position = position, momentum = momentum)
 }
 
+#' Compute logarithm of random_walk proposal density ratio.
+#'
+#' @inherit log_density_ratio_barker return params
+#'
+#' @keywords internal
+log_density_ratio_random_walk <- function(
+    state,
+    proposed_state,
+    target_distribution,
+    scale_and_shape) {
+  0
+}
 
 #' Create a new random walk proposal object.
 #'
@@ -21,26 +33,27 @@ sample_random_walk <- function(
 #' @export
 #'
 #' @examples
-#' target_distribution <- list(
-#'   log_density = function(x) -sum(x^2) / 2
-#' )
-#' proposal <- random_walk_proposal(target_distribution, scale = 1.)
+#' target_distribution <- list(log_density = function(x) -sum(x^2) / 2)
+#' proposal <- random_walk_proposal(scale = 1.)
 #' state <- chain_state(c(0., 0.))
-#' withr::with_seed(876287L, proposed_state <- proposal$sample(state))
-#' log_density_ratio <- proposal$log_density_ratio(state, proposed_state)
+#' 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)
 random_walk_proposal <- function(
-    target_distribution,
     scale = NULL,
     shape = NULL,
     sample_auxiliary = stats::rnorm) {
   scale_and_shape_proposal(
-    sample = function(state, scale_and_shape) {
+    sample = function(state, target_distribution, scale_and_shape) {
       sample_random_walk(
         state, target_distribution, scale_and_shape, sample_auxiliary
       )
     },
-    log_density_ratio = function(state, proposed_state, scale_and_shape) 0,
+    log_density_ratio = log_density_ratio_random_walk,
     scale = scale,
     shape = shape,
     default_target_accept_prob = 0.234,