Accept-Reject Algorithm

Author

Jiachen Ai; UID: 206182615

Problem Statement

You are given a probability density function (PDF) \(f(x)\) that you wish to sample from:

\[ f(x) = \frac{3}{8} x^2, \quad \text{for } x \in [0, 2]. \]

  1. Proposal Distribution: Use a uniform distribution \(g(x) = \frac{1}{2}\) on the interval \([0, 2]\) as the proposal distribution.

Answer:

The proposal distribution \(g(x)\) is a uniform distribution on the interval \([0, 2]\), which can be expressed as: \[ g(x) = \frac{1}{2}, \quad \text{for } x \in [0, 2]. \]

  1. Determine \(c\): Find an appropriate value for \(c\) such that \(c g(x) \geq f(x)\) for all \(x \in [0, 2]\).

Answer:

To determine the value of \(c\), we need to find the maximum value of \(\frac{f(x)}{g(x)}\) for \(x \in [0, 2]\).

Given: \(f(x) = \frac{3}{8} x^2\) and \(g(x) = \frac{1}{2}\)

We have: \(\frac{f(x)}{g(x)} = \frac{\frac{3}{8} x^2}{\frac{1}{2}} = \frac{3}{4} x^2\)

To find the maximum value of \(\frac{f(x)}{g(x)}\), we take the derivative with respect to \(x\) and set it equal to zero: \(\frac{d}{dx} \left( \frac{3}{4} x^2 \right) = \frac{3}{2} x = 0\)

Solving for \(x\), we get \(x = 0\). Since \(\frac{f(x)}{g(x)}\) is monotonically increasing for \(x \in [0, 2]\), the maximum value occurs at \(x = 2\). Therefore, the maximum value of \(\frac{f(x)}{g(x)}\) is: \(\frac{f(2)}{g(2)} = \frac{\frac{3}{8} \cdot 2^2}{\frac{1}{2}} = 3\)

Thus, we choose \(c = 3\) such that \(c g(x) \geq f(x)\) for all \(x \in [0, 2]\).

  1. Implement the Accept-Reject Algorithm: Write out the steps of the accept-reject algorithm for generating random samples from \(f(x)\) using the proposal distribution \(g(x)\) and the value of \(c\) you determined.

Answer:

The Accept-Reject Algorithm for generating random samples from \(f(x)\) using the proposal distribution \(g(x)\) and \(c = 3\) is as follows:

  1. Generate a random sample \(x\) from the uniform distribution \(g(x) = \frac{1}{2}\) on the interval \([0, 2]\).
  2. Generate a random sample \(u\) from the uniform distribution on the interval \([0, 1]\).
  3. Calculate the acceptance probability \(A(x) = \frac{f(x)}{c g(x)} = \frac{\frac{3}{8} x^2}{3 \cdot \frac{1}{2}} = \frac{1}{4} x^2\).
  4. If \(u \leq A(x)\), accept the sample \(x\); otherwise, reject the sample.
  5. Repeat steps 1-4 until the desired number of samples is generated.

The acceptance rate of the accept-reject algorithm can be calculated as the ratio of accepted samples to the total number of samples generated.

  1. Simulate Samples (Programming): Using R, simulate 1000 random samples from the distribution \(f(x)\). Plot a histogram of the generated samples along with the true density \(f(x)\).

Answer:

First, we define the accept-reject function in R that implements the accept-reject algorithm using the proposal distribution \(g(x)\) which is uniform on \([0, 2]\) and the value of \(c = 3\).

accept_reject <- function(fx, n) {
  
  # Simulates a sample of size n from the pdf fx 
  # using the accept-reject algorithm
  x <- numeric(n)
  count <- 0
  iterations <- 0  # track the number of iterations
  # create a list to store the result
  list <- list()
  
  while (count < n) {
    
    # Generate a random sample from the proposal distribution 
    # g(x) which is uniform on [0, 2]
    temp <- runif(1, 0, 2) 
    
    # Generate uniform random value to decide 
    # whether to accept or reject the candidate sample
    y <- runif(1, 0, 1)
    
    iterations <- iterations + 1  # Increment the iteration counter
    
    if (y < fx(temp)) {
      count <- count + 1
      x[count] <- temp
    }
  }
  
  return(list(samples = x, iterations = iterations))
}

Next, we simulate 100000 random samples from the distribution \(f(x)\).

# Define the acceptance probability function A(x) = f(x) / c g(x)
Ax <- function(x) {
  return(1/4 * x^2)
}

# Define the true density function f(x)
fx <- function(x) {
  return(3/8 * x^2)
}


# Simulate 100000 random samples from the distribution f(x)
set.seed(2024)
samples <- accept_reject(Ax, n = 100000)

Finally, we plot a histogram of the generated samples along with the true density \(f(x)\).

# Plot histogram of the generated samples
hist(samples$samples, freq = FALSE, breaks = 50, col = "lightblue", 
     main = "Histogram of Generated Samples vs. True Density",
     xlab = "x", ylab = "Density")

# Add the true density f(x) to the plot
curve(fx(x), add = TRUE, col = "red", lwd = 2)

Evaluate Efficiency:

The acceptance rate of the accept-reject algorithm is calculated above. The choice of \(c\) directly affects the efficiency of the algorithm because it determines how closely the proposal distribution approximates the target distribution.

Answer:

The acceptance rate of the accept-reject algorithm can be calculated as the ratio of accepted samples to the total number of samples generated. In this case, the acceptance rate can be calculated as follows:

# Calculate the acceptance rate
acceptance_rate <- length(samples$samples) / samples$iterations
acceptance_rate
[1] 0.3343755

Bonus - Alternative Proposal Distribution: Suggest an alternative proposal distribution \(g(x)\) that could improve the acceptance rate of the algorithm. Justify why this choice might be better.

Answer:

In order to improve the acceptance rate of the accept-reject algorithm, we need to choose a proposal distribution \(g(x)\) that closely approximates the target distribution \(f(x)\).

Here I suggest using the Beta distribution with parameters \(\alpha = 2\) and \(\beta = 1\) as the proposal distribution. The Beta distribution is a flexible distribution that can be adjusted to closely match the shape of the target distribution. In this case, the Beta(2, 1) distribution can be rescaled to the interval \([0, 2]\) to serve as a proposal distribution for generating samples from the target distribution \(f(x)\).

# Accept-Reject function using Beta distribution
accept_reject <- function(fx, n) {
  
  # Simulates a sample of size n from the pdf fx 
  # using the accept-reject algorithm
  x <- numeric(n)
  count <- 0
  iterations <- 0  # track the number of iterations
  
  while (count < n) {
    
    # Generate a random sample from the Beta(3, 1) distribution, 
    # and rescale it to the interval [0, 2]
    temp <- rbeta(1, 2, 1) * 2
    
    # Generate uniform random value to decide 
    # whether to accept or reject the candidate sample
    y <- runif(1, 0, 1)
    
    iterations <- iterations + 1  # Increment the iteration counter
    
    # Check if we accept the sample based on the acceptance probability
    if (y < fx(temp)) {
      count <- count + 1
      x[count] <- temp
    }
  }
  
  return(list(samples = x, iterations = iterations))
}

Next, we simulate 100000 random samples from the distribution \(f(x)\).

# Define the acceptance probability function A(x) = f(x) / c g(x)
Ax <- function(x) {
  return(1/4 * x^2)
}

# Simulate 100000 random samples from the distribution f(x)
set.seed(2024)
result <- accept_reject(Ax, n = 100000)

Finally, we plot a histogram of the generated samples along with the true density \(f(x)\).

# Plot histogram of the generated samples
hist(result$samples, freq = FALSE, breaks = 50, col = "lightblue", 
     main = "Histogram of Generated Samples vs. True Density",
     xlab = "x", ylab = "Density")

# Add the true density f(x) to the plot
curve(fx(x), add = TRUE, col = "red", lwd = 2)

Evaluate Efficiency:

# Calculate the acceptance rate
acceptance_rate <- length(result$samples) / result$iterations
acceptance_rate
[1] 0.498671