arms_gibbs: Gibbs sampling using ARMS

Description

This function uses ARMS (see also arms) to sample from a multivariate target distribution specified by its (potentially unnormalised) log density using Gibbs sampling. The function updates each argument to the log pdf in turn using ARMS, returning a matrix of samples.

The arguments to this function have the same meaning as for arms, except here they are recycled along the dimension of previous, rather than from sample to sample.

Usage

arms_gibbs(n_samples, previous, log_pdf, lower, upper, initial = NULL,
  n_initial = 10, convex = 0, max_points = 100, metropolis = TRUE,
  include_n_evaluations = FALSE, show_progress = FALSE)

Arguments

n_samples

Number of samples to return.

Starting value for the Gibbs sampler. Vectors of this length are passed to log_pdf as the first argument.

log_pdf

Potentially unnormalised log density of target distribution.

lower

Lower bound of the support of the target distribution.

upper

Upper bound of the support of the target distribution.

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.

show_progress

If TRUE, a progress bar is shown.

Value

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 classic 8schools example from RStan
# https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started
schools_data <- list(
  J = 8,
  y = c(28,  8, -3,  7, -1,  1, 18, 12),
  sigma = c(15, 10, 16, 11,  9, 11, 10, 18)
)

log_pdf <- function(params, p) {
  mu <- params[1]
  tau <- params[2]
  eta <- params[3 : (3 + schools_data$J - 1)]
  output <- (
    sum(dnorm(eta, 0, 1, log = TRUE)) +
    sum(dnorm(
      schools_data$y,
      mu + tau * eta,
      schools_data$sigma,
      log = TRUE
    ))
  )
  return(output)
}

x_start <- c(0, 0, rep(0, schools_data$J))
names(x_start) <- c(
  'mu',
  'tau',
  paste0('eta', 1 : schools_data$J)
)

samples <- arms_gibbs(
  250,
  x_start,
  log_pdf,
  c(-1000, 0, rep(-1000, schools_data$J)),
  1000,
  metropolis = FALSE
)
print(colMeans(samples[51 : 250, ]))
# }