arms: Adaptive Rejection Metropolis Sampling

Description

This function performs Adaptive Rejection Metropolis Sampling to sample from a target distribution specified by its (potentially unnormalised) log density. The function constructs a rejection distribution based on piecewise linear functions that envelop the log density of the target.

If the target is log-concave, the metropolis parameter can be set to FALSE, and an accept-reject sampling scheme is used which yields independent samples.

Otherwise, if metropolis is TRUE, a Metropolis-Hastings step is used to construct a Markov chain with a stationary distribution matching the target. It is possible in this case for the rejection distribution to be a poor proposal, so users should be careful to check the output matches the desired distribution.

All arguments other than n_samples, include_n_evaluations and arguments can be either vectors or lists as appropriate. If they are vectors, they will be recycled in the same manner as, e.g., rnorm. The entries of arguments may be vectors/lists and will also be recycled (see examples).

Usage

arms(n_samples, log_pdf, lower, upper, previous = (upper + lower)/2,
  initial = lower + (1:n_initial) * (upper - lower)/(n_initial + 1),
  n_initial = 10, convex = 0, max_points = max(2 * n_initial + 1,
  100), metropolis = TRUE, include_n_evaluations = FALSE,
  arguments = list())

Arguments

n_samples

Number of samples to return.

log_pdf

Potentially unnormalised log density of target distribution. Can also be a list of functions.

lower

Lower bound of the support of the target distribution.

upper

Upper bound of the support of the target distribution.

The previous value of the Markov chain to be used if metropolis = TRUE.

initial

Initial points with which to build the rejection distribution.

n_initial

Number of points used to form initial; ignored if initial provided.

convex

Convexity adjustment.

max_points

Maximum number of points to allow in the rejection distribution.

metropolis

Whether to use a Metropolis-Hastings step after rejection sampling. Not necessary if the target distribution is log concave.

include_n_evaluations

Whether to return an object specifying the number of function evaluations used.

arguments

List of additional arguments to be passed to log_pdf

Value

Vector or matrix of samples if include_n_evaluations is FALSE, otherwise a list.

References

Gilks, W. R., Best, N. G. and Tan, K. K. C. (1995) Adaptive rejection Metropolis sampling. Applied Statistics, 44, 455-472.

Examples

Run this code

# NOT RUN {
# The normal distribution, which is log concave, so metropolis can be FALSE
result <- arms(
  1000, dnorm, -1000, 1000, metropolis = FALSE,
  arguments = list(log = TRUE), include_n_evaluations = TRUE
)
print(result$n_evaluations)
hist(result$samples, freq = FALSE, br = 20)
curve(dnorm(x), min(result$samples), max(result$samples), col = 'red', add = TRUE)

# Mixture of normals: 0.4 N(-1, 1) + 0.6 N(4, 1). Not log concave.
dnormmixture <- function(x) {
  parts <- log(c(0.4, 0.6)) + dnorm(x, mean = c(-1, 4), log = TRUE)
  log(sum(exp(parts - max(parts)))) + max(parts)
}
samples <- arms(1000, dnormmixture, -1000, 1000)
hist(samples, freq = FALSE)

# List of log pdfs, demonstrating recycling of log_pdf argument
samples <- arms(
  10,
  list(
    function(x) -x ^ 2 / 2,
    function(x) -(x - 10) ^ 2 / 2
  ),
  -1000,
  1000
)
print(samples)

# Another way to achieve the above, this time with recycling in arguments
samples <- arms(
  10, dnorm, -1000, 1000,
  arguments = list(
    mean = c(0, 10), sd = 1, log = TRUE
  )
)
print(samples)
# }