arms function

Adaptive Rejection Metropolis Sampling

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).

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.
• previous: 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

Returns

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

Examples

# 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)

References

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

See Also

http://www1.maths.leeds.ac.uk/~wally.gilks/adaptive.rejection/web_page/Welcome.html