From 20425e3a2cc9edce51f25251234bc12b3f56fa4b Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Mon, 6 Jan 2025 16:29:32 +0000 Subject: [PATCH 1/4] Add vignette demoing changing Barker proposal noise distribution --- vignettes/adjusting-noise-distribution.Rmd | 224 +++++++++++++++++++++ 1 file changed, 224 insertions(+) create mode 100644 vignettes/adjusting-noise-distribution.Rmd diff --git a/vignettes/adjusting-noise-distribution.Rmd b/vignettes/adjusting-noise-distribution.Rmd new file mode 100644 index 0000000..38f3a7e --- /dev/null +++ b/vignettes/adjusting-noise-distribution.Rmd @@ -0,0 +1,224 @@ +--- +title: "Adjusting the noise distribution in the Barker proposal" +output: rmarkdown::html_vignette +bibliography: references.bib +link-citations: true +vignette: > + %\VignetteIndexEntry{Adjusting the noise distribution in the Barker proposal} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +The Barker proposal implementation in the `rmcmc` package in `barker_proposal()` by default uses a standard normal distribution for the auxiliary noise variables generated in the proposal, as suggested by @livingstone2022barker. +The analysis in @vogrinc2023optimal however implies that alternative choices of noise distribution can give improved performance in some situations. +This vignette demonstrates how to adjust the Barker proposal noise distribution. + +```{r setup} +library(rmcmc) +``` + +## Example target distribution + +As a simple example of a target distribution, we consider a $D$-dimensional distribution with product form (independent dimensions) and 'hyperbolic' marginals + +$$\pi(x) \propto \exp\left(-\sum_{d=1}^D(\delta^2 + x_d^2)^{\frac{1}{2}}\right).$$ +The gradient of the corresponding log density function can be derived as +$$\nabla(\log \pi)(x)_d = -x_d / (\delta^2 + x_d^2)^{\frac{1}{2}}.$$ +This can be implemented in R as + +```{r} +delta_sq <- 0.1 +target_distribution <- list( + log_density = function(x) -sum(sqrt(delta_sq + x^2)), + gradient_log_density = function(x) -x / sqrt(delta_sq + x^2) +) +``` + +Here we use $\delta^2 = `r delta_sq`$, corresponding to Example 4 in @vogrinc2023optimal. + +## Creating Barker proposal with a custom noise distribution + +The `barker_proposal()` function accepts an optional argument `sample_auxiliary`, +which can be used to specify the distribution the auxiliary noise variables are drawn from. +The argument should be passed a function accepting a single integer argument, +corresponding to the target distribution dimension, +and returning a vector of auxiliary noise variable samples, +one per target distribution dimension. +By default this function is set to `stats::rnorm()`, +thus using independent standard normal variates for the auxiliary noise variables. + +Below we create a function `sample_bimodal` which instead samples from a bimodal normal mixture, +with two equally-weighted normal components with means $\pm (1 - \sigma^2)$ and common variance $\sigma^2 = 0.01$. + +```{r} +sample_bimodal <- function(dimension) { + sigma <- 0.1 + return( + sample(c(-1, 1), dimension, TRUE) * sqrt(1 - sigma^2) + + stats::rnorm(dimension) * sigma + ) +} +``` + +This choice of bimodal distribution for the auxiliary noise variables is motivated by the analysis in @vogrinc2023optimal, +which shows that for target distributions of product form and using a Barker proposal, +the asymptotic _expect squared jump distance_ (ESJD) a measure of chain sampling performance, +is maximised when the auxiliary noise distribution is chosen to be the Rademacher distribution +(the distribution with half of its mass in atoms at -1 and +1). +The Rademacher distribution itself does not give a practical sampling algorithm, +as the resulting Markov kernel will not be irreducible. +As a more practical alternative @vogrinc2023optimal therefore suggests using the above normal mixture bimodal distribution, +which can be considered a relaxation of the Rademacher distribution. + +We now create two instances of the the Barker proposal, +the first using the default of a standard normal distribution for the auxiliary noise variables, +and the second using the bimodal distribution. + +```{r} +normal_noise_proposal <- barker_proposal() +bimodal_noise_proposal <- barker_proposal(sample_auxiliary = sample_bimodal) +# Above line is equivalent to +# bimodal_noise_proposal <- bimodal_barker_proposal() +``` + +A convenience function `bimodal_barker_proposal()` is also provided in the package to simplify creating a proposal with bimodal noise distribution as above, with optional `sigma` argument specifying the standard deviation of the normal components. + +## Sampling chains with proposals + +Now that we have the target distribution and two proposals specified, +we can sample chains to compare their performance. + +```{r} +set.seed(7861932L) +dimension <- 100 +initial_state <- rnorm(dimension) +adapters <- list(scale_adapter()) +n_warm_up_iteration <- 10000 +n_main_iteration <- 10000 +``` + +Above we specify some common settings for the `sample_chain()` calls, +namely the dimension of the target distribution, +the initial chain state (shared for both chains), +the adapters to use during the warm-up iterations (here we adapt only the proposal scale to achieve a target acceptance rate of 0.57), +and the number of warm-up and main chain iterations. + +We can now sample a chain using the bimodal noise proposal variant + +```{r} +bimodal_noise_results <- sample_chain( + target_distribution = target_distribution, + proposal = bimodal_noise_proposal, + initial_state = initial_state, + n_warm_up_iteration = n_warm_up_iteration, + n_main_iteration = n_main_iteration, + adapters = adapters +) +``` + +and similarly using the default normal noise proposal + +```{r} +normal_noise_results <- sample_chain( + target_distribution = target_distribution, + proposal = normal_noise_proposal, + initial_state = initial_state, + n_warm_up_iteration = n_warm_up_iteration, + n_main_iteration = n_main_iteration, + adapters = adapters +) +``` + +## Comparing performance of proposals + +For a realisation of a Markov chain $X_{1:N}$, +the ESJD metric used to define the notion of sampling efficiency optimized over in @vogrinc2023optimal, +can be estimated as + +$$ + \widehat{\textsf{ESJD}}(X_{1:N}) = \frac{1}{N - 1} \sum_{n=1}^{N-1} \Vert X_{n+1} - X_n \Vert^2 +$$ +with corresponding R implementation as below + +```{r} +expected_square_jumping_distance <- function(traces) { + n_iteration <- nrow(traces) + mean(rowSums(traces[2:n_iteration, ] - traces[1:n_iteration - 1, ])^2) +} +``` + +Applying this function to the chain traces generated using the Barker proposal with bimodal noise, + + +```{r} +cat( + sprintf( + "Expected square jumping distance using normal noise is %.2f", + expected_square_jumping_distance(bimodal_noise_results$traces[, 1:dimension]) + ) +) +``` + +and Barker proposal with normal noise, + +```{r} +cat( + sprintf( + "Expected square jumping distance using bimodal noise is %.2f", + expected_square_jumping_distance(normal_noise_results$traces[, 1:dimension]) + ) +) +``` + +in both cases selecting columns `1:dimension` to exclude the traced log density values, +we find that the proposal with bimodal noise gives roughly a factor two improvement in ESJD compared to using normal noise, +correlating with the results in @vogrinc2023optimal. + + +```{r} +library(posterior) +``` + +Using the `posterior` R package, we can all compare how the chains perform in terms of the estimated _effective sample size_ (ESS) of the estimate of the mean of the target log density. + +First considering the chain using the Barker proposal with bimodal noise + + +```{r} +cat( + sprintf( + "Estimated ESS of mean(target_log_density) using bimodal noise is %.0f", + ess_mean( + extract_variable( + as_draws_matrix(bimodal_noise_results$traces), "target_log_density" + ) + ) + ) +) +``` + +and similarly for the Barker proposal with normal noise + +```{r} +cat( + sprintf( + "Estimated ESS of mean(target_log_density) using normal noise is %.0f", + ess_mean( + extract_variable( + as_draws_matrix(normal_noise_results$traces), "target_log_density" + ) + ) + ) +) +``` + +we again find that the bimodal noise variant appears to show significantly improved sampling efficiency. + +## References From dce5f4b7a67877a121d6028116d273cb7f4ca13c Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Mon, 6 Jan 2025 16:30:04 +0000 Subject: [PATCH 2/4] Decouple descriptions of alt proposals from barker_proposal --- R/hamiltonian.R | 14 +++++++++++++- R/langevin.R | 12 +++++++++++- R/random_walk.R | 7 +++++-- man/hamiltonian_proposal.Rd | 18 +++++++++++++----- man/langevin_proposal.Rd | 16 +++++++++++----- man/random_walk_proposal.Rd | 9 +++------ 6 files changed, 56 insertions(+), 20 deletions(-) diff --git a/R/hamiltonian.R b/R/hamiltonian.R index fe958b9..b7d52fc 100644 --- a/R/hamiltonian.R +++ b/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 diff --git a/R/langevin.R b/R/langevin.R index 6c53280..28ff465 100644 --- a/R/langevin.R +++ b/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 #' diff --git a/R/random_walk.R b/R/random_walk.R index 934c7e6..63f1236 100644 --- a/R/random_walk.R +++ b/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 #' diff --git a/man/hamiltonian_proposal.Rd b/man/hamiltonian_proposal.Rd index a4b1f4a..f09ed6d 100644 --- a/man/hamiltonian_proposal.Rd +++ b/man/hamiltonian_proposal.Rd @@ -56,11 +56,12 @@ value to use for the initial proposal scale parameter. } } \description{ -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 \code{scale} and \code{shape}. At least one -of \code{scale} and \code{shape} must be set before sampling from the proposal or -evaluating the log density ratio. +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. } \examples{ target_distribution <- list( @@ -110,3 +111,10 @@ withr::with_seed( 876287L, proposed_state <- proposal$sample(state, target_distribution) ) } +\references{ +Duane, S., Kennedy, A. D., Pendleton, B. J., & Roweth, D. (1987). +Hybrid Monte Carlo. \emph{Physics Letters B}, 195(2), 216-222. + +Neal, R. M. (2011). MCMC Using Hamiltonian Dynamics. In \emph{Handbook +of Markov Chain Monte Carlo} (pp. 113-162). Chapman and Hall/CRC. +} diff --git a/man/langevin_proposal.Rd b/man/langevin_proposal.Rd index d6573ca..7860be6 100644 --- a/man/langevin_proposal.Rd +++ b/man/langevin_proposal.Rd @@ -33,11 +33,8 @@ value to use for the initial proposal scale parameter. } } \description{ -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 \code{scale} and \code{shape}. At least one -of \code{scale} and \code{shape} must be set before sampling from the proposal or -evaluating the log density ratio. +The Langevin proposal is a gradient-based proposal corresponding to a +Euler-Maruyama time discretisation of a Langevin diffusion. } \examples{ target_distribution <- list( @@ -54,3 +51,12 @@ log_density_ratio <- proposal$log_density_ratio( ) proposal$update(scale = 0.5) } +\references{ +Besag, J. (1994). "Comments on "Representations of knowledge in +complex systems" by U. Grenander and MI Miller". \emph{Journal of the Royal +Statistical Society, Series B}. 56: 591–592. + +Roberts, G. O., & Tweedie, R. L. (1996). Exponential convergence +of Langevin distributions and their discrete approximations. \emph{Bernoulli} 2 +(4), 341 - 363. +} diff --git a/man/random_walk_proposal.Rd b/man/random_walk_proposal.Rd index cdc5e3a..3bc00b2 100644 --- a/man/random_walk_proposal.Rd +++ b/man/random_walk_proposal.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/random_walk.R \name{random_walk_proposal} \alias{random_walk_proposal} -\title{Create a new random walk proposal object.} +\title{Create a new (Gaussian) random walk proposal object.} \usage{ random_walk_proposal( scale = NULL, @@ -37,11 +37,8 @@ value to use for the initial proposal scale parameter. } } \description{ -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 \code{scale} and \code{shape}. At least one -of \code{scale} and \code{shape} must be set before sampling from the proposal or -evaluating the log density ratio. +The Gaussian random walk proposal samples a new proposed state by perturbing +the current state with zero-mean normally distributed noise. } \examples{ target_distribution <- list(log_density = function(x) -sum(x^2) / 2) From 75d61f83c75585a89e335fc688cea126897a7a3d Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Mon, 6 Jan 2025 16:30:48 +0000 Subject: [PATCH 3/4] Add bimodal Barker proposal convenience function --- NAMESPACE | 1 + R/barker.R | 77 ++++++++++++++++++++++++++++++--- man/barker_proposal.Rd | 19 ++++++--- man/bimodal_barker_proposal.Rd | 78 ++++++++++++++++++++++++++++++++++ tests/testthat/test-barker.R | 8 ++++ 5 files changed, 173 insertions(+), 10 deletions(-) create mode 100644 man/bimodal_barker_proposal.Rd diff --git a/NAMESPACE b/NAMESPACE index 2bb2300..87e344e 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/barker.R b/R/barker.R index 342f38c..50c077d 100644 --- a/R/barker.R +++ b/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 + ) +} diff --git a/man/barker_proposal.Rd b/man/barker_proposal.Rd index 73f1d09..03edd8f 100644 --- a/man/barker_proposal.Rd +++ b/man/barker_proposal.Rd @@ -41,11 +41,14 @@ value to use for the initial proposal scale parameter. } } \description{ -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 \code{scale} and \code{shape}. At least one -of \code{scale} and \code{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. +} +\details{ +For more details see the vignette: +\code{vignette("barker-proposal", package = "rmcmc")} } \examples{ target_distribution <- list( @@ -62,3 +65,9 @@ log_density_ratio <- proposal$log_density_ratio( ) proposal$update(scale = 0.5) } +\references{ +Livingstone, S., & Zanella, G. (2022). The Barker proposal: +combining robustness and efficiency in gradient-based MCMC. \emph{Journal of the +Royal Statistical Society Series B: Statistical Methodology}, 84(2), +496-523. +} diff --git a/man/bimodal_barker_proposal.Rd b/man/bimodal_barker_proposal.Rd new file mode 100644 index 0000000..ae111c6 --- /dev/null +++ b/man/bimodal_barker_proposal.Rd @@ -0,0 +1,78 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/barker.R +\name{bimodal_barker_proposal} +\alias{bimodal_barker_proposal} +\title{Create a new Barker proposal object with bimodal noise distribution.} +\usage{ +bimodal_barker_proposal( + sigma = 0.1, + scale = NULL, + shape = NULL, + sample_uniform = stats::runif +) +} +\arguments{ +\item{sigma}{Standard deviation of equally-weighted normal components in +bimodal auxiliary noise distribution, with corresponding means of +\verb{±sqrt(1 - sigma^2)}.} + +\item{scale}{Scale parameter of proposal distribution. A non-negative scalar +value determining scale of steps proposed.} + +\item{shape}{Shape parameter of proposal distribution. Either a vector +corresponding to a diagonal shape matrix with per-dimension scaling +factors, or a matrix allowing arbitrary linear transformations.} + +\item{sample_uniform}{Function which generates a random vector from standard +uniform distribution given an integer size.} +} +\value{ +Proposal object. A list with entries +\itemize{ +\item \code{sample}: a function to generate sample from proposal distribution given +current chain state, +\item \code{log_density_ratio}: a function to compute log density ratio for proposal +for a given pair of current and proposed chain states, +\item \code{update}: a function to update parameters of proposal, +\item \code{parameters}: a function to return list of current parameter values. +\item \code{default_target_accept_prob}: a function returning the default target +acceptance rate to use for any scale adaptation. +\item \code{default_initial_scale}: a function which given a dimension gives a default +value to use for the initial proposal scale parameter. +} +} +\description{ +Convenience function for creating a Barker proposal with bimodal auxiliary +noise variable distribution, corresponding to equally-weighted normal +components with shared variance \code{sigma} and means \verb{±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. +} +\details{ +For more details see the vignette: +\code{vignette("adjusting-noise-distribution", package = "rmcmc")} +} +\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) +} +\references{ +Vogrinc, J., Livingstone, S., & Zanella, G. (2023). Optimal +design of the Barker proposal and other locally balanced +Metropolis–Hastings algorithms. \emph{Biometrika}, 110(3), 579-595. +} +\seealso{ +\code{\link[=barker_proposal]{barker_proposal()}} +} diff --git a/tests/testthat/test-barker.R b/tests/testthat/test-barker.R index 5dc1e2d..4b5a513 100644 --- a/tests/testthat/test-barker.R +++ b/tests/testthat/test-barker.R @@ -5,3 +5,11 @@ test_scale_and_shape_proposal( dimensions = c(1L, 2L), scales = c(0.5, 1., 2.) ) + +test_scale_and_shape_proposal( + bimodal_barker_proposal, + proposal_name = "Bimodal Barker proposal", + target_distribution = standard_normal_target_distribution(), + dimensions = c(1L, 2L), + scales = c(0.5, 1., 2.) +) From da2ec804f9b73ebb1f5163dc04f00a2f22464853 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Mon, 6 Jan 2025 16:48:42 +0000 Subject: [PATCH 4/4] Fix linter errors --- R/barker.R | 4 ++-- vignettes/adjusting-noise-distribution.Rmd | 8 +++----- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/R/barker.R b/R/barker.R index 50c077d..60b9a41 100644 --- a/R/barker.R +++ b/R/barker.R @@ -179,8 +179,8 @@ bimodal_barker_proposal <- function( sample_uniform = stats::runif) { sample_bimodal <- function(dimension) { return( - sample(c(-1, 1), dimension, TRUE) * sqrt(1 - sigma^2) - + stats::rnorm(dimension) * sigma + sample(c(-1, 1), dimension, TRUE) * sqrt(1 - sigma^2) + + stats::rnorm(dimension) * sigma ) } barker_proposal( diff --git a/vignettes/adjusting-noise-distribution.Rmd b/vignettes/adjusting-noise-distribution.Rmd index 38f3a7e..117557f 100644 --- a/vignettes/adjusting-noise-distribution.Rmd +++ b/vignettes/adjusting-noise-distribution.Rmd @@ -61,8 +61,8 @@ with two equally-weighted normal components with means $\pm (1 - \sigma^2)$ and sample_bimodal <- function(dimension) { sigma <- 0.1 return( - sample(c(-1, 1), dimension, TRUE) * sqrt(1 - sigma^2) - + stats::rnorm(dimension) * sigma + sample(c(-1, 1), dimension, TRUE) * sqrt(1 - sigma^2) + + stats::rnorm(dimension) * sigma ) } ``` @@ -84,11 +84,9 @@ and the second using the bimodal distribution. ```{r} normal_noise_proposal <- barker_proposal() bimodal_noise_proposal <- barker_proposal(sample_auxiliary = sample_bimodal) -# Above line is equivalent to -# bimodal_noise_proposal <- bimodal_barker_proposal() ``` -A convenience function `bimodal_barker_proposal()` is also provided in the package to simplify creating a proposal with bimodal noise distribution as above, with optional `sigma` argument specifying the standard deviation of the normal components. +A convenience function `bimodal_barker_proposal()` is also provided in the package to simplify creating a proposal with bimodal noise distribution as above, with optional `sigma` argument specifying the standard deviation of the normal components. That is the second line above can equivalently be written `bimodal_noise_proposal <- bimodal_barker_proposal()` without the need to define the `sample_bimodal` function. ## Sampling chains with proposals