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:
Suppose we would like to estimate \(S= E[x^2]\) where \(X\) has density proportional to \(q(x) = e^ \frac{-|x|^3}{3}\ \)
Let the target function be \(f(x) = e^ \frac{-|x|^3}{3}\ \)
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)
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
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