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.
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:
lb = 0
to set the
pseudo-target’s lower bound in the call to
pseudo_list
.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.107652
Neal, R. M. (2003), “Slice sampling,” The Annals of Statistics, 31, 705-767. https://doi.org/10.1214/aos/1056562461
Nishihara, 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
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.