--- title: "Package `qslice`" author: "Matthew J. Heiner, David B. Dahl, and Samuel B. Johnson" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Qslice package} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library("qslice") ``` 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)} \, .$$ ```{r target} ## 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. ```{r step} ## 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 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) ``` # 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: ```{r pseudo1} 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. ```{r vis_pseudo1, fig.width=4.5, fig.height=3} utility_pseudo(pseudo = pseu, log_target = ltarget, type = "function", utility_type = "AUC", plot = TRUE) ``` 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. ```{r qslice1, fig.show='hold'} ## 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 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 ## Null hypothesis: Samples are iid Gamma(alpha, 1) ks.test(x_sample[seq(1, n_iter, by = 10)], pgamma, shape = alpha) ``` ## 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). ```{r pseudo2, fig.width=4.5, fig.height=3} 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. ```{r pseudo3} 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) names(pseu_opt$pseudo) pseu_opt$pseudo$txt ``` We conclude by running the quantile slice sampler again with an optimized pseudo-target. ```{r qslice2, fig.show='hold'} ## 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 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 ## Null hypothesis: Samples are iid Gamma(alpha, 1) ks.test(x_sample[seq(1, n_iter, by = 10)], pgamma, shape = alpha) ``` ## 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, ```{r diagnostics_arg, eval = FALSE} slice_quantile(..., diagnostics = TRUE) ``` or through global options ```{r diagnostics_options, eval = FALSE} options(qslice.diagnostics = TRUE) ``` The status of the global option can be checked by running ```{r diagnostics_getOption, eval = FALSE} 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 - Heiner, M. J., Johnson, S. B., Christensen, J. R., and Dahl, D. B. (2026+), "Quantile Slice Sampling," *arXiv preprint arXiv:2407.12608*. - Li, Y. and Walker, S. G. (2023), "A latent slice sampling algorithm," *Computational Statistics and Data Analysis*, 179, 107652. - Neal, R. M. (2003), "Slice sampling," *The Annals of Statistics*, 31, 705-767. - Nishihara, R., Murray, I., and Adams, R. P. (2014), "Parallel MCMC with Generalized Elliptical Slice Sampling," *Journal of Machine Learning Research*, 15, 2087-2112.