Processing math: 91%

Qslice package

Matthew J. Heiner, David B. Dahl, and Samuel B. Johnson

2024-05-29

The quantile slice sampler of Heiner et al. (2024+), as well as other popular slice samplers, are implemented in the qslice package. Several utility functions for specifying pseudo-target distributions, for diagnostics and for tuning, are also included.

Other implemented methods include the generalized elliptical (Nishihara et al., 2014), latent (Li and Walker, 2023), and stepping-out (Neal, 2003) slice samplers, and independence Metropolis-Hastings sampler.

All sampling functions in the qslice package perform a single scan that can be iteratively called as part of an MCMC routine.

Slice sampling

Let g(x) represent an unnormalized target density. Augment the target with a latent random variable VX=xUniform(0,g(x)), yielding a joint density for X and V that is proportional to a constant.

Simple slice samplers first draw from the conditional distribution of V given above, and then from the conditional of XV=v, which is uniform on the set A={x:v<g(x)} known as the slice region. If A is available analytically, a custom slice sampler can be implemented. The slice samplers available in this package repeatedly sample until a draw on the slice region is found. All use the shrinking method in Neal (2003), which we demonstrate below.

A Gibbs sampler drawing VtXt1 followed by XtVt produces a sequence of draws {Xt} whose distribution converges to the marginal density proportional to g(x).

Example

Each sampling method requires a function to evaluate the natural logarithm of a (typically unnormalized) target density. For the sake of illustration, we use a target that is proportional to a gamma density (with shape α and rate 1):

g(x)=xα1exp(x)1(x>0).

## Target is a Gamma(shape = alpha, rate = 1) distribution
alpha <- 2.5
ltarget <- function(x) ifelse(x > 0.0, (alpha - 1.0) * log(x) - x, -Inf)

Note that log-target functions in qslice rely on lexical scoping. For example, the parameter alpha is not passed as an argument to function ltarget. If the target represents a full conditional distribution that changes at each iteration of an MCMC algorithm, the user must either: 1. change the value of α within the environment in which ltarget is defined or 2. redefine the ltarget function.

We first demonstrate the stepping-out-and-shrinkage sampler of Neal (2003), a general tool that is robust to the choice of its tuning parameter w. The sampling function slice_stepping_out performs a single draw, so we embed it within an MCMC algorithm.

## Set up MCMC
n_iter <- 1e3
x_sample <- numeric(n_iter + 1)
x_sample[1] <- 0.5  # initialize
n_eval <- 0  # count target evaluations

## Run step-and-shrink slice sampler
for (i in 2:(n_iter+1)) {
  state <- slice_stepping_out(x = x_sample[i-1], 
                              log_target = ltarget,
                              w = 2.0)
  n_eval <- n_eval + state$nEvaluations
  x_sample[i] <- state$x
}

## Check samples
n_eval / n_iter  # target evaluations per MCMC iteration
#> [1] 6.679
hist(x_sample, freq = FALSE, n = 30)
curve(dgamma(x, shape = alpha), col = "blue", lwd = 2, add = TRUE)


## Null hypothesis: Samples are iid Gamma(alpha, 1)
ks.test(x_sample[seq(1, n_iter, by = 5)], pgamma, shape = alpha)
#> 
#>  Asymptotic one-sample Kolmogorov-Smirnov test
#> 
#> data:  x_sample[seq(1, n_iter, by = 5)]
#> D = 0.076299, p-value = 0.1947
#> alternative hypothesis: two-sided

Quantile slice sampler

The quantile slice sampler uses an approximation to the target density, called a pseudo-target. Let ˆG() represent the distribution function (CDF) of the pseudo-target with inverse (i.e., quantile) function ˆG1() and density ˆg(). We use the probability integral transform to define a new random variable ψ=ˆG(X)[0,1]. The original target can be written as g(x)=g(x)ˆg(x)ˆg(x), which after transformation becomes h(ψ)=g(ˆG1(ψ))ˆg(ˆG1(ψ))1(0ψ1). If the pseudo-target approximates the target well, h(ψ) will be close to a constant function, yielding an efficient slice sampler. The quantile slice sampler uses target h(ψ) and adaptively shrinks (around the previously sampled value of ψ) to draw a sample belonging the slice region Aψ[0,1].

Example

The quantile slice sampler requires a pseudo-target. This is specified in the qslice package with a list containing two functions: the log-density ld and inverse-CDF (quantile) q functions. qslice offers pseudo-target constructor functions with options for several distribution families. See help(pseudo_list). Here we use a truncated t distribution:

pseu <- pseudo_list(family = "t", 
                    params = list(loc = 0, sc = 3, degf = 1), 
                    lb = 0)

When specifying a pseudo-target, please keep in mind that:

  1. The support of the pseudo-target should match (or exceed) that of the original target. Note the use of lb = 0 to set the pseudo-target’s lower bound in the call to pseudo_list.
  2. The pseudo-target should have tails that are at least as heavy as the original target. For this reason, we recommend the Student-t distribution.

We visualize g(x), ˆg(x) and h(ψ) below.

utility_pseudo(pseudo = pseu, log_target = ltarget, 
               type = "function", 
               utility_type = "AUC", plot = TRUE)

#> [1] 0.4907685

The function utility_pseudo measures how close h() is to a constant function. Here we use AUC (0,1], defined as the area under the curve h(ψ)/max.

The function slice_quantile performs a single draw using the quantile slice sampler. It also outputs draws of \psi (called u in the output), which can be used as a diagnostic for pseudo-target fitness.

## Set up MCMC
n_iter <- 1e3
x_sample <- psi_sample <- numeric(n_iter + 1)
x_sample[1] <- psi_sample[1] <- 0.5  # initialize
n_eval <- 0  # count target evaluations

## Run step-and-shrink slice sampler
for (i in 2:(n_iter+1)) {
  state <- slice_quantile(x = x_sample[i-1], 
                          log_target = ltarget,
                          pseudo = pseu)
  n_eval <- n_eval + state$nEvaluations
  x_sample[i] <- state$x
  psi_sample[i] <- state$u
}

## Check samples
n_eval / n_iter  # target evaluations per iteration of MCMC
#> [1] 2.669

hist(x_sample, freq = FALSE, n = 30)
curve(dgamma(x, shape = alpha), col = "blue", lwd = 2, add = TRUE)

hist(psi_sample, freq = FALSE, n = 30)

auc(u = psi_sample) # calculate AUC from transformed samples
#> [1] 0.4634259

## Null hypothesis: Samples are iid Gamma(alpha, 1)
ks.test(x_sample[seq(1, n_iter, by = 10)], pgamma, shape = alpha)
#> 
#>  Asymptotic one-sample Kolmogorov-Smirnov test
#> 
#> data:  x_sample[seq(1, n_iter, by = 10)]
#> D = 0.085335, p-value = 0.4602
#> alternative hypothesis: two-sided

Optimal pseudo-targets and tuning

The quantile slice sampler is effective even with crude approximations to the target. We can also find a pseudo-target within the Student-t family that maximizes AUC.

pseu_opt <- pseudo_opt(log_target = ltarget, 
                       type = "function",
                       family = "t", degf = c(1, 5, 20), 
                       lb = 0,
                       utility_type = "AUC", plot = TRUE)

Alternatively, we can use initial samples from the target to specify a pseudo-target. This gives a simple procedure for “tuning” the sampler.

pseu_opt <- pseudo_opt(samples = x_sample, 
                       type = "samples",
                       family = "t", degf = c(1, 5, 20), 
                       lb = 0,
                       utility_type = "AUC", plot = TRUE, nbins = 20)


names(pseu_opt)
#> [1] "pseudo"       "utility"      "utility_type" "opt"          "nbins"       
#> [6] "tol_int"      "tol_opt"
names(pseu_opt$pseudo)
#> [1] "d"      "ld"     "q"      "p"      "txt"    "params" "lb"     "ub"    
#> [9] "family"
pseu_opt$pseudo$txt
#> [1] "t(loc = 1.51, sc = 1.94, degf = 20) I(0 < x < Inf)"

We conclude by running the quantile slice sampler again with an optimized pseudo-target.

## Set up MCMC
n_iter <- 1e3
x_sample <- psi_sample <- numeric(n_iter + 1)
x_sample[1] <- psi_sample[1] <- 0.5  # initialize
n_eval <- 0  # count target evaluations

## Run step-and-shrink slice sampler
for (i in 2:(n_iter+1)) {
  state <- slice_quantile(x = x_sample[i-1], 
                          log_target = ltarget,
                          pseudo = pseu_opt$pseudo)
  n_eval <- n_eval + state$nEvaluations
  x_sample[i] <- state$x
  psi_sample[i] <- state$u
}

## Check samples
n_eval / n_iter  # target evaluations per iteration of MCMC
#> [1] 2.254

hist(x_sample, freq = FALSE, n = 30)
curve(dgamma(x, shape = alpha), col = "blue", lwd = 2, add = TRUE)

hist(psi_sample, freq = FALSE, n = 30)

auc(u = psi_sample) # calculate AUC from transformed samples
#> [1] 0.5853801

## Null hypothesis: Samples are iid Gamma(alpha, 1)
ks.test(x_sample[seq(1, n_iter, by = 10)], pgamma, shape = alpha)
#> 
#>  Asymptotic one-sample Kolmogorov-Smirnov test
#> 
#> data:  x_sample[seq(1, n_iter, by = 10)]
#> D = 0.073568, p-value = 0.6513
#> alternative hypothesis: two-sided

References