Solution 4

Exercise 4.1

Consider the Laplace family of distribution, \(\mathsf{Laplace}(\nu, \tau)\), with density \[\begin{align*} g(x; \nu, \tau) = \frac{1}{2\tau} \exp\left(- \frac{|x-\nu|}{\tau}\right), \qquad \nu \in \mathbb{R}, \tau > 0 \end{align*}\] as a candidate distribution for rejection sampling from \(\mathsf{Gauss}(0,1)\).

  1. Provide an inversion sampling algorithm to generate from \(\mathsf{Laplace}(\nu, \tau)\).
  2. Can you use the proposal to generate from a standard Gaussian? for Student-\(t\) with 1 degree of freedom? Justify your answer.
  3. Consider as proposal a location-scale version of the Student-t with \(3\) degrees of freedom. Find the optimal location and scale parameters and the upper bound \(C\) for your choice.
  4. Use the accept-reject to simulate 1000 independent observations and compute the empirical acceptance rate.

Solution.

  1. The distribution function is \[\begin{align*} F(x) = \begin{cases} \frac{1}{2} \exp\left(\frac{x-\nu}{\tau}\right) & x \leq \nu\\ 1 - \frac{1}{2} \exp\left(-\frac{x-\nu}{\tau}\right) & x >\nu\\ \end{cases} \end{align*}\] and using the quantile transform, set \(X=\nu + \tau \log(2U)\) if \(U \leq 0.5\) and \(X=\nu - \tau\log(2-2U)\) if \(U > 0.5\) for \(U \sim \mathsf{unif}(0,1)\).
  2. The Gaussian has lighter tail than the Laplace, so this won’t work. The Cauchy distribution would be a suitable candidate, albeit too heavy tailed.
  3. The optimal value for the location of the Student-\(t\) would be \(\nu\) (e.g., zero for the standard Laplace). We compute the optimal scale via in the code below: \[ \mathrm{argmin}_{\sigma \in \mathbb{R}_{+}}\mathrm{argmax}_{x \in \mathbb{R}} \{\log f(x) - \log g(x; \sigma)\} \]
#' Laplace density
dlaplace <- function(x, loc = 0, scale = 1, log = FALSE){
 stopifnot(scale > 0)
 logdens <-  -log(2*scale) - abs(x-loc)/scale
 if(log){
 return(logdens)
 } else{
 return(exp(logdens))
 }
}
dstudent <- function(x, loc = 0, scale = 1, df = 1, log = FALSE){
  logdens <- -log(scale) + dt(x = (x - loc)/scale, df = df, log = TRUE)
   if(log){
    return(logdens)
   } else{
   return(exp(logdens))
   }
}
# For each value of the scale sigma,
# find the minimum value of x (typically at zero)
opt <- optimize(f = function(sigma){
  optimize(f = function(x){
  dlaplace(x, log = TRUE) - 
    dstudent(x, scale = sigma, df = 3, log = TRUE)}, 
  maximum = TRUE, 
  interval = c(-100, 100))$objective
}, interval = c(0.1,10))
(C <- exp(opt$objective))
[1] 1.261368
(sigma <- opt$minimum)
[1] 0.9272381
# Simulate from accept-reject
ntrials <- 1.1*C*1000
# Simulate from location-scale student
candidate <- sigma*rt(n = ntrials, df = 3)
# Compute log of acceptance rate
logR <- dlaplace(candidate, log = TRUE) - 
  dstudent(candidate, scale = sigma, df = 3, log = TRUE) 
samp <- candidate[logR >= log(C) + -rexp(ntrials)]
# Monte Carlo estimator of the acceptance rate
ntrials/length(samp)
[1] 1.261368
# Plot density
library(ggplot2)
ggplot(data = data.frame(x = samp),
       mapping = aes(x = x)) +
  geom_density() +
  theme_classic()

  1. The Monte Carlo acceptance rate is 1.26, compared with the analytical bound found via numerical optimization of 1.26, to two significant digits.