The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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 \(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)\).

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 \(\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 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
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 \(\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]\).

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)\), \(\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
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

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.