library(ggplot2)
library(patchwork)
set.seed(536)
theme_set(theme_bw(base_size = 14))
# Utility function for density comparison plots
plot_comparison <- function(samples, true_density, title = "Generated vs True Distribution",
x_range = NULL) {
if (is.null(x_range)) {
x_range <- range(samples)
}
df <- data.frame(x = samples)
ggplot(df, aes(x = x)) +
geom_histogram(aes(y = after_stat(density)),
bins = 50, fill = "lightblue", color = "black", alpha = 0.7) +
stat_function(fun = true_density, color = "red", linewidth = 1.2,
xlim = x_range) +
labs(title = title, x = "Value", y = "Density") +
theme(plot.title = element_text(hjust = 0.5))
}Lab 2: Accept-Reject Sampling
STAT 536 — Statistical Computing
Learning Goals
By the end of this lab session, you will:
- Implement the Accept-Reject algorithm from scratch
- Compare different instrumental distributions and their efficiency
- Implement Envelope Accept-Reject (EAR) for improved efficiency
- Apply these methods to practical problems
- Gain intuition through visualization
Setup
Start by loading the required libraries and utility functions:
Exercise 1: Your First Accept-Reject Implementation
We generate samples from a Beta(2,5) distribution using a Uniform(0,1) instrumental distribution.
Understanding the Setup
The Accept-Reject algorithm needs:
- Target density \(f(x)\), what we want to sample from (Beta(2,5))
- Instrumental density \(g(x)\), what we can easily sample from (Uniform(0,1))
- Constant \(M\) such that \(f(x) \leq M \cdot g(x)\) for all \(x\)
Finding M
set.seed(536)
# Define the Beta(2,5) density
dbeta_target <- function(x) dbeta(x, 2, 5)
# TODO: Find where Beta(2,5) achieves its maximum
# HINT: For Beta(a,b), the mode is at (a-1)/(a+b-2) when a,b > 1
mode_location <- ___ # Fill in the formula
# TODO: Find M (the maximum value of the target density)
M <- ___ # What is the density at the mode?
print(paste("Mode of Beta(2,5) is at x =", mode_location))
print(paste("Beta(2,5) density at mode =", M))
print(paste("Therefore, M =", M))
print(paste("Theoretical acceptance rate = 1/M =", round(1/M, 4)))
# Visualize the envelope
x_seq <- seq(0, 1, length.out = 1000)
df_viz <- data.frame(
x = x_seq,
target = dbeta_target(x_seq),
envelope = rep(M, length(x_seq))
)
# Create the visualization plot
# (Plot code provided; run it after filling in M)Implementing Accept-Reject
Now implement the Accept-Reject algorithm:
set.seed(536)
beta_ar_uniform <- function(n) {
samples <- numeric(n)
num_accepted <- 0
num_trials <- 0
# TODO: Calculate M
mode_val <- ___ # Use the mode formula from the previous chunk
M <- ___ # Maximum of Beta(2,5)
# Generate n accepted samples
while (num_accepted < n) {
# Step 1: Generate candidate Y from g (uniform)
Y <- ___ # HINT: Use runif()
# Step 2: Generate U for acceptance decision
U <- ___ # HINT: Another uniform random number
# Step 3: Accept or reject
# HINT: Accept if U <= f(Y)/(M*g(Y))
# Since g(Y) = 1 for uniform, this simplifies.
if (___) { # Fill in the acceptance condition
num_accepted <- num_accepted + 1
samples[num_accepted] <- Y
}
num_trials <- num_trials + 1
}
return(list(
samples = samples,
acceptance_rate = n / num_trials,
trials = num_trials,
M = M
))
}
# Test your implementation
result <- beta_ar_uniform(1000)
print(paste("Empirical acceptance rate =", round(result$acceptance_rate, 4)))
print(paste("Total trials needed for 1000 samples:", result$trials))Visualize Your Results
set.seed(536)
# Use the provided plot_comparison function
plot_comparison(result$samples, function(x) dbeta(x, 2, 5),
"Accept-Reject Beta(2,5) Samples", c(0, 1))
# Create a Q-Q plot for validation
theoretical_quantiles <- qbeta((1:1000)/(1001), 2, 5)
# TODO: Create Q-Q plot comparing theoretical_quantiles with sort(result$samples)Check Your Understanding:
- Is your empirical acceptance rate close to 1/M?
- Does the Q-Q plot show points along the diagonal?
Exercise 2: Comparing Instrumental Distributions
We see how different instrumental distributions affect efficiency for Beta(2,5).
Implement Three Different Methods
set.seed(536)
# Method 1: Uniform(0,1) - already done above
# Method 2: Beta(2,2) instrumental
beta_ar_beta22 <- function(n) {
samples <- numeric(n)
num_accepted <- 0
num_trials <- 0
# TODO: Find M for Beta(2,2) instrumental
# HINT: M = max(f(x)/g(x)) where f=Beta(2,5), g=Beta(2,2)
# You need to find where this ratio is maximized.
x_test <- seq(0.01, 0.99, by = 0.01)
ratio <- ___ # Compute dbeta(x_test, 2, 5) / dbeta(x_test, 2, 2)
M <- ___ # Find the maximum ratio
x_max <- ___ # Find where maximum occurs
print(paste("Beta(2,2) instrumental: M =", M, "at x =", x_max))
while (num_accepted < n) {
# TODO: Generate Y from Beta(2,2)
Y <- ___
U <- runif(1)
# TODO: Acceptance condition
if (___) {
num_accepted <- num_accepted + 1
samples[num_accepted] <- Y
}
num_trials <- num_trials + 1
}
return(list(samples = samples,
acceptance_rate = n / num_trials,
M = M))
}
# Method 3: Beta(1,3) instrumental
beta_ar_beta13 <- function(n) {
# TODO: Similar to Method 2, but with Beta(1,3)
# This one should work better. Why?
# Your code here...
}
# Test all three methods
n_test <- 5000
result_uniform <- beta_ar_uniform(n_test)
result_beta22 <- beta_ar_beta22(n_test)
result_beta13 <- beta_ar_beta13(n_test)Compare Efficiency
set.seed(536)
# Create comparison table
comparison_df <- data.frame(
Method = c("Uniform(0,1)", "Beta(2,2)", "Beta(1,3)"),
M = c(result_uniform$M, result_beta22$M, result_beta13$M),
Theoretical_Efficiency = ___, # HINT: 1/M for each
Empirical_Efficiency = ___ # Use acceptance_rate from each result
)
print(comparison_df)Think About:
- Which method has the highest acceptance rate?
- Why might Beta(2,2) perform worse than expected?
- What property of Beta(1,3) makes it work well?
Report. In one or two sentences, rank the three instrumentals by the acceptance rates you actually obtained, quoting the numbers, and say whether they match the theoretical \(1/M\) values.
Exercise 3: Truncated Distributions
Generate samples from a Normal distribution truncated to \([a, b]\).
Implement Truncated Normal Sampling
set.seed(536)
truncated_normal_ar <- function(n, mu = 0, sigma = 1, a = -2, b = 2) {
samples <- numeric(n)
num_accepted <- 0
num_trials <- 0
# Truncation probability for normalization
p_trunc <- pnorm(b, mu, sigma) - pnorm(a, mu, sigma)
# TODO: Find M
# HINT: The maximum of truncated normal is at:
# - mu if a < mu < b
# - a if mu <= a
# - b if mu >= b
if (mu < a) {
x_max <- ___
} else if (mu > b) {
x_max <- ___
} else {
x_max <- ___
}
# Maximum density value
M <- dnorm(x_max, mu, sigma) / p_trunc
while (num_accepted < n) {
# TODO: Generate candidate from uniform on [a, b]
Y <- ___
U <- runif(1)
# TODO: Acceptance condition
# Target density: dnorm(Y, mu, sigma) / p_trunc
# Instrumental: 1/(b-a)
if (___) {
num_accepted <- num_accepted + 1
samples[num_accepted] <- Y
}
num_trials <- num_trials + 1
}
return(list(samples = samples,
acceptance_rate = n / num_trials))
}
# Test with different truncation ranges
result_narrow <- truncated_normal_ar(1000, mu=0, sigma=1, a=-1, b=1)
result_wide <- truncated_normal_ar(1000, mu=0, sigma=1, a=-3, b=3)
result_tail <- truncated_normal_ar(1000, mu=0, sigma=1, a=2, b=4)
print(paste("Narrow [-1,1]: Acceptance =", result_narrow$acceptance_rate))
print(paste("Wide [-3,3]: Acceptance =", result_wide$acceptance_rate))
print(paste("Tail [2,4]: Acceptance =", result_tail$acceptance_rate))Observe:
How does the acceptance rate change with truncation range?
Why is sampling from the tail less efficient?
Exercise 4: Envelope Accept-Reject (EAR)
Implement EAR with both upper and lower squeeze functions. The Envelope Accept-Reject method uses BOTH lower and upper squeeze functions to minimize evaluations of the expensive target density.
Theory
For target \(f(x)\) and instrumental \(g(x)\), we introduce squeeze functions such that: \[g_l(x) \leq f(x) \leq g_u(x) \leq M \cdot g(x)\]
where:
- \(g_l(x)\) is the lower squeeze (cheap to evaluate)
- \(g_u(x)\) is the upper squeeze (cheap to evaluate)
- \(f(x)\) is the target (expensive to evaluate)
The algorithm:
- Generate \(Y \sim g\), \(U \sim \text{Uniform}(0,1)\).
- If \(U \leq g_l(Y)/(M \cdot g(Y))\), accept \(Y\) immediately (quick accept).
- Else if \(U > g_u(Y)/(M \cdot g(Y))\), reject immediately (quick reject).
- Else evaluate \(f(Y)\) and accept if \(U \leq f(Y)/(M \cdot g(Y))\).
Basic EAR Structure
set.seed(536)
beta_ear <- function(n) {
samples <- numeric(n)
num_accepted <- 0
num_trials <- 0
quick_accepts <- 0
quick_rejects <- 0
# Setup for Beta(2,5) with uniform instrumental
mode_val <- (2-1)/(2+5-2)
M <- dbeta(mode_val, 2, 5)
# Define lower squeeze function
gl <- function(x) {
# TODO: Simple parabola touching f at mode
# HINT: Use a*(x - mode_val)^2 + b
# where b = f(mode_val) and a makes it touch at x=0 and x=1
peak <- dbeta(mode_val, 2, 5)
width <- 0.15 # Adjust for tightness
# Return parabola value
___
}
# Define upper squeeze function
gu <- function(x) {
# TODO: Piecewise linear that is always above f(x)
# HINT: Two lines meeting at the mode
peak <- dbeta(mode_val, 2, 5)
if (x <= mode_val) {
return(___) # Line from (0,0) to (mode_val, peak)
} else {
return(___) # Line from (mode_val, peak) to (1,0)
}
}
while (num_accepted < n) {
Y <- runif(1)
U <- runif(1)
# Quick accept check (using lower squeeze)
if (___) { # HINT: U <= gl(Y)/M
num_accepted <- num_accepted + 1
samples[num_accepted] <- Y
quick_accepts <- quick_accepts + 1
}
# Quick reject check (using upper squeeze)
else if (___) { # HINT: U > gu(Y)/M
quick_rejects <- quick_rejects + 1
}
# Need to evaluate f(Y)
else {
if (___) { # Standard Accept-Reject check
num_accepted <- num_accepted + 1
samples[num_accepted] <- Y
}
}
num_trials <- num_trials + 1
}
print(paste("Quick accepts:", quick_accepts, "(",
round(100*quick_accepts/n, 1), "%)"))
print(paste("Quick rejects:", quick_rejects, "(",
round(100*quick_rejects/num_trials, 1), "%)"))
return(list(samples = samples,
acceptance_rate = n / num_trials))
}
# Test EAR
result_ear <- beta_ear(1000)Key Points:
- How many function evaluations did you save?
- When is EAR most beneficial?
Exercise 5: Bayesian Application
Use Accept-Reject for sampling from a posterior distribution.
Poisson-Gamma Example
set.seed(536)
# Observed data: counts with unknown rate lambda
observed_counts <- c(3, 5, 4, 6, 3, 4, 5, 7, 4, 3)
# Prior: Gamma(alpha=2, beta=1)
prior_alpha <- 2
prior_beta <- 1
# Posterior is Gamma(alpha + sum(x), beta + n)
post_alpha <- prior_alpha + sum(observed_counts)
post_beta <- prior_beta + length(observed_counts)
# True posterior for comparison
true_posterior <- function(lambda) {
dgamma(lambda, post_alpha, post_beta)
}
# Sample using Accept-Reject with shifted Gamma instrumental
bayesian_ar <- function(n) {
samples <- numeric(n)
num_accepted <- 0
# TODO: Choose good instrumental distribution
# HINT: Use a Gamma that is slightly wider than the posterior
inst_alpha <- ___ # Try post_alpha - 0.5
inst_beta <- ___ # Try post_beta - 1
# Find M
x_test <- seq(0.1, 10, by = 0.1)
ratio <- dgamma(x_test, post_alpha, post_beta) /
dgamma(x_test, inst_alpha, inst_beta)
M <- max(ratio)
while (num_accepted < n) {
# TODO: Complete the Accept-Reject steps
Y <- ___ # Sample from instrumental
U <- runif(1)
if (___) { # Acceptance condition
num_accepted <- num_accepted + 1
samples[num_accepted] <- Y
}
}
return(samples)
}
# Generate posterior samples
posterior_samples <- bayesian_ar(5000)
# Compare with true posterior
plot_comparison(posterior_samples, true_posterior,
"Posterior Sampling via Accept-Reject")Wrap-up Questions
After completing these exercises, you should be able to answer:
- Why does the choice of instrumental distribution matter so much?
- Think about the areas between f(x) and M*g(x).
- When would you use EAR instead of standard Accept-Reject?
- Consider the cost of evaluating f(x).
- What happens to Accept-Reject efficiency in high dimensions?
- If 1D efficiency is p, what is the d-dimensional efficiency?
- How do you choose a good instrumental distribution?
- What properties should it share with the target?