Overview

Rejection sampling is a method of drawing points from a target distribution by sampling from an easier distribution and correcting the sampling probability through the rejection of some candidates.

Suppose we are interested in drawing points from \(f(x)\). Let \(g\) denote another density from which we know how to sample and for which we can easily calculate \(g(x)\). Let \(e\) denote an envelope such that \(e(x) = \frac{g(x)}{\alpha}\geq f(x) \forall x\) for which \(f(x) > 0\) for a given constant \(\alpha \leq 1\).

Then rejection sampling can be completed in the following steps:

  1. Sample \(Y\) ~ \(g\)
  2. Sample \(U\) ~ \(Unif(0,1)\)
  3. Reject \(Y\) if \(U > \frac{f(Y)}{e(Y)}\), otherwise keep \(Y\)
  4. Repeat for desired sample size

Demonstration

Suppose we would like to estimate \(S= E[x^2]\) where \(X\) has density proportional to \(q(x) = e^ \frac{-|x|^3}{3}\ \)

Target Function

Let the target function be \(f(x) = e^ \frac{-|x|^3}{3}\ \)

Envelope Function

Consider using the standard normal distribution \(N(0,1)\)
\(g = \frac{1}{\sqrt{2\pi}}e^{\frac{-x^2}{2}}\)

To determine \(e\), we must find \(\alpha\) such that \(\frac{g(x)}{\alpha} > f(x)\). Choose \(\alpha = 0.3\).

f <- function(x) {         
  exp(-abs(x^3)/3)
}
e <- function(x) {     
  dnorm(x)/0.3
}
curve(e(x), from = -2, to = 2, col = "blue")
curve(f(x), add = T)

Rejection Sampling

Follow the steps listed above for a sample size of 100,000. Below are a few points drawn using this method.

n = 100000
y = rnorm(n, 0, 1)
u = runif(n)
x = y[u < f(y)/e(y)]
head(x)
## [1]  1.14850878 -0.16326236 -0.17638245  0.04891202  0.05516292  0.65279403

The acceptance ratio is

length(x)/n
## [1] 0.7733


Calculate \(E[x^2]\)

Now that we have our points from the sample, square each accepted x, and take the mean to get \(E[x^2]\).

mean(x^2)
## [1] 0.7801127