slice_quantile_mv {qslice}

Multivariate Quantile Slice Sampler

Description

Quantile slice sampler for a random vector (Heiner et al., 2024+). The pseudo-target is specified through independent univariate distributions.

Usage

slice_quantile_mv(x, log_target, pseudo)

Arguments

x

The current state (as a numeric vector).

log_target

A function taking numeric vector that evaluates the log-target density, returning a numeric scalar.

pseudo

List of length equal to the number of dimensions in x. Each element is itself a list that specifies the pseudo-target for the corresponding dimension with functions ld that evaluates the log density, p that evaluates the CDF, and q that evaluates the quantile (inverse-CDF) function.

Value

A list containing three elements: "x" is the new state, "u" is the value of the CDF of the psuedo-target associated with the returned value, inverse CDF method, and "nEvaluations is the number of evaluations of the target function used to obtain the new state.

References

Heiner, M. J., Johnson, S. B., Christensen, J. R., and Dahl, D. B. (2024+), "Quantile Slice Sampling," arXiv preprint arXiv:###

Examples

lf <- function(x) dbeta(x[1], 3, 4, log = TRUE) + dbeta(x[2], 5, 3, log = TRUE)
ps_shsc <- list(c(2, 2), c(2, 1))
ps <- list(
  list(ld = function(x) dbeta(x, ps_shsc[[1]][1], ps_shsc[[1]][2], log = TRUE),
       p = function(x) pbeta(x, ps_shsc[[1]][1], ps_shsc[[1]][2]),
       q = function(x) qbeta(x, ps_shsc[[1]][1], ps_shsc[[1]][2]) ),
  list(ld = function(x) dbeta(x, ps_shsc[[2]][1], ps_shsc[[2]][2], log = TRUE),
       p = function(x) pbeta(x, ps_shsc[[2]][1], ps_shsc[[2]][2]),
       q = function(x) qbeta(x, ps_shsc[[2]][1], ps_shsc[[2]][2]) )
  )
n_iter <- 10 # set to 1e4 for more complete illustration
draws <- matrix(0.2, nrow = n_iter, ncol = 2)
draws_u <- draws
draws_u[1,] <- sapply(1:length(ps), function(k) ps[[k]]$p(draws[1,k]))
nEvaluations <- 0L
for (i in seq.int(2, n_iter)) {
  out <- slice_quantile_mv(draws[i - 1, ], log_target = lf, pseudo = ps)
  draws[i,] <- out$x
  draws_u[i,] <- out$u
  nEvaluations <- nEvaluations + out$nEvaluations
  cat(i, '\r')
}
nEvaluations / (nrow(draws) - 1)
plot(draws[,1], draws[,2], xlim = c(0, 1))
hist(draws[,1], freq = FALSE); curve(dbeta(x, 3, 4), col = "blue", add = TRUE)
hist(draws[,2], freq = FALSE); curve(dbeta(x, 5, 3), col = "blue", add = TRUE)
plot(draws_u[,1], draws_u[,2], xlim = c(0, 1))
hist(draws_u[,1], freq = FALSE)
hist(draws_u[,2], freq = FALSE)
auc(u = draws_u[,1])
auc(u = draws_u[,2])