Package qslice

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

2026-10-21

The qslice R package implements the quantile slice sampler of Heiner et al. (2026+) along with other popular slice samplers. 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, as well as independence Metropolis-Hastings and random-walk Metropolis samplers.

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 \(V \mid X = x \sim \text{Uniform}(0, g(x))\). This results in a joint density for \(X\) and \(V\) that is proportional to a constant and yields a marginal density of \(X\) that is proportional to \(g\).

Uniform slice samplers first draw from the conditional distribution of \(V\) given above, and then from the conditional of \(X \mid V = v\), which is uniform on the set \(A = \{ x : v < g(x) \}\), known as the slice region. If \(A\) is available analytically, a custom ideal slice sampler can be implemented. The samplers available in the qslice package include general-purpose hybrid steps that repeatedly sample until a draw on the slice region is found. Each uses some variant of the bracket shrinkage procedure of Neal (2003).

Example

Sampling methods in qslice require 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 \(\alpha\) and rate 1):

\[ g(x) = x^{\alpha - 1} \exp(-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 accept x as the only argument. In this 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 \(\alpha\) 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 often 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 <- 3e3
x_sample <- numeric(n_iter)
state <- list(x = 0.5)  # initialize
n_eval <- 0  # count target evaluations

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

## Check samples
n_eval / n_iter  # target evaluations per MCMC iteration
#> [1] 6.617333
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) (thin for better approximation)
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.035794, p-value = 0.8368
#> alternative hypothesis: two-sided

Quantile slice sampler

The quantile slice sampler uses an approximation to the target density, called a pseudo-target. Let \(\hat{G}(\cdot)\) represent the distribution function (CDF) of the pseudo-target with inverse (i.e., quantile) function \(\hat{G}^{-1}(\cdot)\) and density \(\hat{g}(\cdot)\). We define a new random variable \(\psi = \hat{G}(X) \in (0,1)\). The original target can be written as \[g(x) = \frac{g(x)}{\hat{g}(x)}\hat{g}(x) \, ,\] which by the probability integral transform becomes \[h(\psi) = \frac{g(\hat{G}^{-1}(\psi))}{\hat{g}(\hat{G}^{-1}(\psi))} \, 1_{(0 \le \psi \le 1)} \, .\] Quantile slice sampling is equivalent to uniform slice sampling \(\psi\) with target \(h\), and transforming samples back to \(X\), which have stationary density proportional to \(g\).

If the pseudo-target approximates the target well, \(h(\psi)\) will be close to a constant function, yielding an efficient quantile slice sampler.

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, df = 1), 
                    lb = 0)

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

  1. The support of the original target should be equal to or a subset of the support of the pseudo-target. Note the use of lb = 0 in pseudo_list to restrict the lower bound of the pseudo-target’s support.
  2. The tails of the pseudo-target should be at least as heavy as those of the original target. For this reason, we recommend the Student-\(t\) distribution.

We visualize \(g(x)\), \(\hat{g}(x)\) and \(h(\psi)\) below.

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

#> [1] 0.5073589

The function utility_pseudo measures how close \(h(\cdot)\) is to a constant function. Here we use AUC \(\in (0, 1]\), defined as the area under the curve \(h(\psi) / \max_z(h(z))\). Alternatively, use utility_type = "MSW" to calculate the expected total width of the slice region in the \(\psi\) space, which is the interval \((0,1)\).

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 <- 3e3
x_sample <- psi_sample <- numeric(n_iter)
state <- list(x = 0.5)  # initialize
n_eval <- 0  # count target evaluations

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

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

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, xlim = c(0,1))

utility_shrinkslice(u = psi_sample, utility_type = "AUC") # estimate AUC from transformed samples
#> [1] 0.5070984

## 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.057492, p-value = 0.2746
#> alternative hypothesis: two-sided

Optimal pseudo-targets and tuning

The quantile slice sampler works even with crude approximations to the target. We can also find a pseudo-target within the Student-\(t\) family that maximizes AUC (or MSW).

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

Alternatively, we can tune the sampler using initial samples from a pilot run to specify a pseudo-target.

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


names(pseu_opt)
#> [1] "pseudo"          "utility"         "utility_type"    "opt"            
#> [5] "nbins"           "utility_options" "tol_opt"
names(pseu_opt$pseudo)
#> [1] "d"      "ld"     "q"      "p"      "txt"    "params" "lb"     "ub"    
#> [9] "family"
pseu_opt$pseudo$txt
#> [1] "t(loc = 0.87, sc = 2.18, df = 5) I(0 < x < Inf)"

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

## Set up MCMC
n_iter <- 3e3
x_sample <- psi_sample <- numeric(n_iter)
state <- list(x = 0.5)  # initialize
n_eval <- 0  # count target evaluations

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

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

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, xlim = c(0,1))

utility_shrinkslice(u = psi_sample, utility_type = "AUC") # estimate AUC from transformed samples
#> [1] 0.7234131

## 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.03881, p-value = 0.7569
#> alternative hypothesis: two-sided

Diagnostic mode

Samplers in qslice and pseudo-target utility functions have an optional diagnostic mode that returns detailed error messages tracking proposals, target evaluations, and steps. Diagnostic mode can be set directly in the call to each sampler, for example,

slice_quantile(..., diagnostics = TRUE)

or through global options

options(qslice.diagnostics = TRUE)

The status of the global option can be checked by running

getOption("qslice.diagnostics")

Because diagnostic mode validates all calls to target densities, pseudo-target densities, and pseudo-target quantile functions, as well as inputs and outputs at each step, diagnostic mode can significantly slow sampling. We recommend its use only for testing, validation, and debugging.

References