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.

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))\), 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 \(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 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 \(V_t \mid X_{t-1}\) followed by \(X_t \mid V_t\) produces a sequence of draws \(\{X_t\}\) whose distribution converges to the marginal density proportional to \(g(x)\).

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 \(\alpha\) and rate 1):

\[ g(x) = x^{\alpha - 1} \exp(-x) \, 1_{(x>0)} \, .\]

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

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 \(\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 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
<- 1e3
n_iter <- numeric(n_iter + 1)
x_sample 1] <- 0.5 # initialize
x_sample[<- 0 # count target evaluations
n_eval
## Run step-and-shrink slice sampler
for (i in 2:(n_iter+1)) {
<- slice_stepping_out(x = x_sample[i-1],
state log_target = ltarget,
w = 2.0)
<- n_eval + state$nEvaluations
n_eval <- state$x
x_sample[i]
}
## Check samples
/ n_iter # target evaluations per MCMC iteration
n_eval #> [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
```

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 use the probability integral transform to 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 after transformation becomes \[h(\psi) = \frac{g(\hat{G}^{-1}(\psi))}{\hat{g}(\hat{G}^{-1}(\psi))} \, 1_{(0 \le \psi \le 1)} \, .\] If the pseudo-target approximates the target well, \(h(\psi)\) will be close to a constant function, yielding an efficient slice sampler. The quantile slice sampler uses target \(h(\psi)\) and adaptively shrinks (around the previously sampled value of \(\psi\)) to draw a sample belonging the slice region \(A_\psi \subseteq [0,1]\).

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:

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

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

- 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`

. - 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)\), \(\hat{g}(x)\) and \(h(\psi)\) 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(\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))\).

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
<- 1e3
n_iter <- psi_sample <- numeric(n_iter + 1)
x_sample 1] <- psi_sample[1] <- 0.5 # initialize
x_sample[<- 0 # count target evaluations
n_eval
## Run step-and-shrink slice sampler
for (i in 2:(n_iter+1)) {
<- slice_quantile(x = x_sample[i-1],
state log_target = ltarget,
pseudo = pseu)
<- n_eval + state$nEvaluations
n_eval <- state$x
x_sample[i] <- state$u
psi_sample[i]
}
## Check samples
/ n_iter # target evaluations per iteration of MCMC
n_eval #> [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
```

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.

```
<- pseudo_opt(log_target = ltarget,
pseu_opt 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.

```
<- pseudo_opt(samples = x_sample,
pseu_opt 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"
```

```
$pseudo$txt
pseu_opt#> [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
<- 1e3
n_iter <- psi_sample <- numeric(n_iter + 1)
x_sample 1] <- psi_sample[1] <- 0.5 # initialize
x_sample[<- 0 # count target evaluations
n_eval
## Run step-and-shrink slice sampler
for (i in 2:(n_iter+1)) {
<- slice_quantile(x = x_sample[i-1],
state log_target = ltarget,
pseudo = pseu_opt$pseudo)
<- n_eval + state$nEvaluations
n_eval <- state$x
x_sample[i] <- state$u
psi_sample[i]
}
## Check samples
/ n_iter # target evaluations per iteration of MCMC
n_eval #> [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
```

Heiner, M. J., Johnson, S. B., Christensen, J. R., and Dahl, D. B. (2024+), “Quantile Slice Sampling,”

*arXiv preprint arXiv:###*.Li, Y. and Walker, S. G. (2023), “A latent slice sampling algorithm,”

*Computational Statistics and Data Analysis*, 179, 107652. https://doi.org/10.1016/j.csda.2022.107652Neal, R. M. (2003), “Slice sampling,”

*The Annals of Statistics*, 31, 705-767. https://doi.org/10.1214/aos/1056562461Nishihara, R., Murray, I., and Adams, R. P. (2014), “Parallel MCMC with Generalized Elliptical Slice Sampling,”

*Journal of Machine Learning Research*, 15, 2087-2112. https://jmlr.org/papers/v15/nishihara14a.html