Slice sampler

Reference

library(progress)

Example 1

Sample from the standard normal density,

\[ f(x) \propto \mathrm{exp}(-x^2/2), ~~ -\infty < x < \infty. \]

set.seed(0)

# total number of iterations for slice sampler
N = 10000

# initializing the space for the samples
us = rep(1, N)  # auxiliary variable
xs = rep(1, N)

# slice sampler
for(i in 2:N){
  us[i] = runif(1, min = 0, max = exp(-0.5*(xs[i-1]^2)))
  xs[i] = runif(1, min = -sqrt(-2*log(us[i])), max = sqrt(-2*log(us[i])))
}

true = seq(-5, 5, 0.01)
true_f = dnorm(true)

hist(xs, freq=F, breaks=200, main = "Histogram of slice samples", xlab = "x")
lines(true, true_f, type="l", col="red")

legend("topright", legend = "True density", col = "red", lty = 1)

Example 2

Sample from exponential(1),

\[ f(x) \propto \mathrm{exp}(-x), ~~ 0 < x < \infty. \]

set.seed(0)

# total number of iterations for slice sampler
N = 10000

# initializing the space for the samples
us = rep(1, N)  # auxiliary variable
xs = rep(1, N)

# slice sampler
for(i in 2:N){
  us[i] = runif(1, min = 0, max = exp(-xs[i-1]))
  xs[i] = runif(1, min = 0, max = -log(us[i]))
}

true = seq(0, 10, 0.01)
true_f = dgamma(true, shape = 1, scale = 1)

hist(xs, freq=F, breaks=200, main = "Histogram of slice samples", xlab=expression(x))
lines(true, true_f, type="l", col="red")

legend("topright", legend = "True density", col = "red", lty = 1)

Example 3 (Numerical method)

Let \(Y_i \overset{ind}{\sim} N(\theta, 1)\) and \(\theta \sim \mathrm{exp}(1)\). Then,

\[ P(\theta|y) \propto \prod_{i=1}^{n} \left [ N(y_i; \theta, 1) \right ] \mathrm{exp}(\theta; 1). \]

Let \(\theta = 0.2\). Since it is hard to derive the inverse function for the slice sampler from,

\[ \begin{aligned} u|\theta &\sim U(0, P(\theta|y)), \\ \theta|u &\sim U \left \{ \theta: u < P(\theta|y) \right \}, \end{aligned} \]

we solve it numerically using the function uniroot() for the two end points of \(\theta\).

set.seed(0)

### Observations
n = 500
true_theta = 0.2
y = rnorm(n, mean = true_theta, sd = 1)

### PDF of the posterior dstn
f = function(theta, y. = y){
    exp(sum(dnorm(y., mean = theta, sd = 1, log = TRUE)) + dexp(abs(theta), log = TRUE))
}

### Function to numerically find endpoints
A = function(u, xx, f. = f){
  left_endpoint = uniroot(function(x) f.(x) - u, c(-10^10, xx)) 
  right_endpoint = uniroot(function(x) f.(x) - u, c(10^10, xx))
  return(c(left_endpoint$root, right_endpoint$root))
}

### Slice sampler
slice_sampler = function(N, init_theta, target, A){
  ## Initialization
  u = rep(NA, N)
  theta = rep(NA, N)
  theta[1] = init_theta
  u[1] = runif(1, 0, target(theta[1]))
  
  pb = progress_bar$new(
    format = "  MCMC [:bar] :percent eta: :eta",
    total = N-1, clear = FALSE, width= 60)
  
  ## MCMC algorithm iterations
  for(i in 2:N){
    pb$tick()
    
    u[i] = runif(1, 0, target(theta[i-1]))
    endpoints = A(u[i], theta[i-1])
    theta[i] = runif(1, endpoints[1], endpoints[2])
  }
  
  return(theta)
}

set.seed(0)

N = 5000  # number of iterations for the MCMC algorithm
ptm = proc.time()  # to time the algorithm
theta = slice_sampler(N, mean(y), f, A)

proc.time() - ptm 

   user  system elapsed 
 15.960   0.537  19.359 

### Traceplot
plot(1:N, theta, type = "l", xlab = "Iterations", ylab = expression(theta))

### Posterior mean
mean(theta)

[1] 0.1980653

Example 4 (Augmentation method)

Let \(Y_i \overset{ind}{\sim} N(\theta, 1)\) and \(\theta \sim \mathrm{exp}(1)\). Then,

\[ P(\theta|y) \propto \prod_{i=1}^{n} \left [ N(y_i; \theta, 1) \right ] \mathrm{exp}(\theta; 1). \]

Unlike example 3, we will now use augmentation,

\[ P(u, \theta) \propto P(\theta)I(0 < u < P(y|\theta)). \]

The full conditional distributions are now,

\[ \begin{aligned} u|\theta, y &\sim U(0, P(y|\theta)), \\ \theta|u, y &\sim P(\theta)I(u < P(y|\theta)). \end{aligned} \]

If \(P(\theta)\) is unimodal, the conditional dstn for \(\theta|u, y\) is equivalent to

\[ P(\theta)I(\theta_L(u) < \theta < \theta_U(u)), \]

for some bounds \(\theta_L(u)\) and \(\theta_U(u)\) which depend on \(u\). To learn these bounds we can sample from \(P(\theta)\) and update the bounds. For example, if \(\theta^{(t-1)}\) is our current value in the chain, we know \(u < P(y|\theta^{(t-1)})\), or equivalently,

\[ \theta_L(u) < \theta^{(t-1)} < \theta_U(u). \]

Letting \(u^{(t)}\) be the current value for the auxiliary variable and setting \(\theta_L(u^{(t)})\) and \(\theta_U(u^{(t)})\) for the lower and upper bounds of the support for \(\theta\) respectively, the slice sampler algorithm goes as follows.

For \(t = 1, ..., T\),

1. Sample \(\theta^*\) from

\[ \theta^* \sim P(\theta)I(\theta_L(u^{(t)}) < \theta < \theta_U(u^{(t)})). \]

2. Set \(\theta^{(t)} = \theta^*\) if

\[ u^{(t)} < P(y|\theta^*), \]

otherwise set

\[ \begin{aligned} \theta_L(u^{(t)}) &= \theta^* ~~ \mathrm{if} ~~ \theta^* < \theta^{(t-1)}, \\ \theta_U(u^{(t)}) &= \theta^* ~~ \mathrm{if} ~~ \theta^* > \theta^{(t-1)}, \end{aligned} \]

and return to Step (1).

### Slice sampler
aug_slice_sampler = function(N, init_theta, like, qprior){ 
  ## Initialization
  u = rep(NA, N)
  theta = rep(NA, N)
  theta[1] = init_theta
  u[1] = runif(1, 0, like(theta[1]))
  
  pb = progress_bar$new(
    format = "  MCMC [:bar] :percent eta: :eta",
    total = N-1, clear = FALSE, width= 60)
  
  ## MCMC algorithm iterations
  for (i in 2:N){
    pb$tick()
    
    u[i] = runif(1, 0, like(theta[i - 1])) 
    success = FALSE
    endpoints = 0:1
    
    while(!success){
    # Inverse CDF
    up = runif(1, endpoints[1], endpoints[2]) 
    theta[i] = qprior(up)
    if (u[i] < like(theta[i])) { 
      success = TRUE
    } else {
      if(theta[i] < theta[i - 1]){endpoints[1] = up}
      else{endpoints[2] = up}
      }
    } 
  }
  return(theta) 
}

set.seed(0)

### Observations
n = 500
true_theta = 0.2
y = rnorm(n, mean = true_theta, sd = 1)

### Likelihood function for a given theta
like = function(tht) exp(sum(dnorm(y, mean = tht, sd = 1, log = TRUE)))

### inverse CDF of prior
qprior = function(tht) qexp(tht)

set.seed(0)

N = 5000  # number of iterations for the MCMC algorithm
ptm = proc.time()  # to time the algorithm
theta = aug_slice_sampler(N, mean(y), like, qprior)

proc.time() - ptm 

   user  system elapsed 
  5.552   0.289   8.261 

We can see that this is much faster than the slice sampler with numerical method.

### Traceplot
plot(1:N, theta, type = "l", xlab = "Iterations", ylab = expression(theta))

### Posterior mean
mean(theta)

[1] 0.1984192