Question 1

Suppose that you can easily generate random variables from the density \(g(x)\), \(x > 0\). How would you generate random variables from density \[f(x) = \frac{g(x)}{m (1 + x)} \,,\] where \[ m = \int_0^\infty \frac{g(y)}{1 + y} dy \,. \]

Question 1

Rejection sampling

• \(\sup_x f(x) / g(x) = 1 / m\)
• Acceptance region: \(U \leq 1 / (1 + X)\)

Alogrithm:

1. Sample \(X\) according to \(g(x)\)
2. Sample \(U \sim \mathcal{U}(0,1)\)
3. If \(U \leq 1 / (1 + X)\), accept, else reject

Question 2

Generate a sample from the density \[ f(x) = \frac{1}{\log(2)(x + 1)(x + 2)}, \quad x > 0 \] using an acceptance-rejection (AR) method and the density \[ g(x) = \frac{1}{(1 + x)^2} \,. \] A direct inverse transformation method was used in Exercise 3, 2f. Discuss the pros and cons of each method.

Question 2

Question 2

Acceptance region

\[ \begin{aligned} c &= \sup_{x > 0} \frac{f(x)}{g(x)} \\ &= \sup_{x > 0} \frac{1}{\log(2)} \frac{(x + 1)^2}{(x + 1) (x + 2)} \\ &= \frac{1}{\log(2)} \sup_{x > 0} \frac{x + 1}{x + 2} \\ &= \frac{1}{\log(2)} \end{aligned} \]

Acceptance region: \(U \leq \frac{X + 1}{X + 2}\)

Question 2

Question 2

f <- function(x){
    1 / (log(2) * (x + 1) * (x + 2))
}

g <- function(x){
    1 / (1 + x)^2
}

c <- 1 / log(2) 

G_inv <- function(p){
    p / (1 - p)
}

F_inv <- function(p){
    (1 - 2^p) / (2^(p - 1) - 1)
}

N <- 10^6
set.seed(123) 

U_1 <- runif(N) 
U_2 <- runif(N) 
X_1 <- F_inv(U_1) 
X_2 <- G_inv(U_1) 
X_2 <- X_2[U_2 <= f(X_2) / (c * g(X_2))] 

ar <- length(X_2) / N 
th_ar <- 1 / c
print(rbind(th_ar, ar)) 

           [,1]
th_ar 0.6931472
ar    0.6925490

curve(f, from = 0, to = F_inv(.975),  col = "black", lwd = 2, lty = 2)  
lines(density(X_1[X_1 < F_inv(.975)]), col = "red", lwd = 2, xlab = "", main = "") 
lines(density(X_2[X_2 < F_inv(.975)]), col = "blue", lwd = 2) 

legend("topright", 
    legend = c("pdf", "Inv. transform", "Rejection sampling"), 
    col = c("black", "red", "blue"), 
    lty = c(2,  1, 1), 
    lwd = rep(2, 3)) 

Question 3

Let \(X\) be a standard normal random variable. Let \(Y = X − 1\), conditionally on \(X\) being between \(1\) and \(2\). We would like to generate random variables with density the same as \(Y\)

\[ f(y) = \frac{\exp \left\{- \frac{(1 + y)^2}{2} \right\}}{\sqrt{2 \pi} (\Psi(2) - \Psi(1))} \]

The natural thing to do is to sample standard normal random variables and accept only those that lie between \(1\) and \(2\) and return \(X − 1\) for the accepted values. As this is wasteful, an obvious way to double the acceptance rate is to accept the variables that lie between \(−2\) and \(−1\) as well as those between \(1\) and \(2\) and return \(|X| − 1\) for those values. An alternative is using an acceptance-rejection method based on the density

\[ g(y) = \frac{a \exp\{-ay \}}{1 - \exp\{-a\}} \,, \] for \(0 < y < 1.\) Find the best value for \(a\). Is this method (using the optimal \(a\)) better than the previous one?

Question 3

Rejection sampling

1. Find a goood upper bound \(c(a)\) for \(\sup_{y \in (0, 1)} \, f(y) / g(y)\)
2. Minimize \(c(a)\) with respect to \(a\)
3. Set up rejection sampling as per usual

Question 3

Upper bound \(c(a)\)

1. \[ \begin{aligned} \underset{y \in (0,1)}{\sup} \, \frac{f(y)}{g(y)} &= \underbrace{\frac{1 - \exp\{-a\}}{a \sqrt{2 \pi} (\Psi(2) - \Psi(1))}}_{C} \underset{y \in (0,1)}{\sup} \, \exp \left\{ - \frac{(1 + y)^2}{2} + ay \right\} \\ &\leq C \underset{y \in \Re}{\sup} \, \exp \left\{ - \frac{(1 + y)^2}{2} + ay \right\} \\ \end{aligned} \]

2. \[ \begin{aligned} \arg\underset{y \in \Re}{ \min} \, \exp \left\{ - \frac{(1 + y)^2}{2} + ay \right\} &= \arg\underset{y \in \Re}{ \min} \, - \frac{(1 + y)^2}{2} + ay \\ &= \arg\underset{y \in \Re}{\textrm{solve}} \, \{a - 1 - y = 0 \} \end{aligned} \]

3. \[ \underset{y \in (0,1)}{\sup} \, \frac{f(y)}{g(y)} \leq C \exp\left\{\frac{a^2}{2} - a \right\} \equiv c(a) \]

Question 3

Minimising \(c(a)\)

Use R

c_a <- function(a){
    C <- (1 - exp(-a)) / (a * sqrt(2 * pi) * (pnorm(2) - pnorm(1))) 
    
    C * exp(a^2 / 2 - a) 
}

opt_a <- optimize(c_a, lower = 0, upper = 10, maximum = FALSE)$minimum

cat("Optimal a: ", opt_a) 

Optimal a:  1.387895

Question 3

Inverse transform sampling from envelope

\[ \begin{aligned} G(z) &= \int_0^z \frac{a \exp\{-xa\}}{1 - \exp\{-a\}} dx \\ &= -\left[\frac{ \exp\{-xa\}}{1 - \exp\{-a\}} \right]_0^z \\ &= \frac{1 - \exp\{-xa\} }{1 - \exp\{-a\}} \end{aligned} \]

\[ \begin{aligned} G(x)^{-1} &= - \frac{\log \left(1 - p(1 - \exp\{-a\}) \right)}{a} \end{aligned} \]

Question 3

Business as usual

f <- function(y){
    ds <- rep(0, length = length(y))
    inds <- (y > 0) & (y < 1)
    ds[inds] <- dnorm(1 + y[inds]) / (pnorm(2) - pnorm(1)) 
    ds
}

g <- function(y, a){
    ds <- rep(0, length = length(y))
    inds <- (y > 0) & (y < 1)
    ds[inds] <-  a * exp(-a * y[inds]) / (1 - exp(-a)) 
    ds
}

G_inv <- function(p, a){
    - log(1 - p * (1 - exp(-a))) / a 
}

N <- 10^5 
set.seed(123)

# Rejection sampling 
U_1 <- runif(N) 
U_2 <- runif(N)
X <- G_inv(U_1, opt_a) 
Y <- X[U_2 < exp(-(1 + X)^2 / 2 + opt_a * X + opt_a - opt_a^2 / 2)]

# Natural method 1 & 2
X_2 <- rnorm(N) 
Y_2 <- X_2[(X_2 < 2) & (X_2 > 1)] - 1 
X_3 <- abs(X_2)
Y_3 <- X_3[(X_3 < 2) & (X_3 > 1)] - 1 

cat("Rejection sampling AR: ", length(Y) / N, ", Natural method AR: ", length(Y_2) / N, ", Natural method 2 AR: ", length(Y_3) / N)

Rejection sampling AR:  0.96499 , Natural method AR:  0.13852 , Natural method 2 AR:  0.27333

curve(f, from = 0, to = 1,  col = "black", lwd = 2, lty = 2)  
lines(density(Y), col = "red", lwd = 2, xlab = "", main = "") 
lines(density(Y_2), col = "blue", lwd = 2) 
lines(density(Y_3), col = "green", lwd = 2) 

legend("topright", 
    legend = c("pdf", "Rejection sampling", "NM 1", "NM 2"), 
    col = c("black", "red", "blue", "green"), 
    lty = c(2, 1, 1, 1), 
    lwd = rep(2, 4))