arms_gibbs: Gibbs sampling using ARMS

In armspp: Adaptive Rejection Metropolis Sampling (ARMS) via 'Rcpp'

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

See Also

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

Examples

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

armspp documentation built on May 24, 2019, 5:05 p.m.