Adaptive squeezed rejection sampling is a method of drawing points from a target distribution, and goes a step further than rejection sampling by utilizing an automatic envelope generation strategy for squeezed rejection sampling.
Suppose we are interested in drawing points from \(f(x)\), a concave function. 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\). It is simpler to generate this envelope function in the log space. Take \(n\) points on \(log(f(x))\) and connect their tangent lines to determine \(log(e(x))\). This ensures that when exponentiated, the envelope function encompasses \(f(x)\).
We will also define a squeeze function \(s(x)\) such that \(s(x) \leq f(x) \forall x\) for which \(f(x) > 0\). Using the selected points from generating the envelope function, connect the points to determine \(log(s(x))\). This ensures that when exponentiated, the squeeze function is below \(f(x)\).
Then adaptive 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}\ \)
Then \(log(f(x)) = \frac{-|x|^3}{3}\ \)
Select points \((-1\), \(-\frac{1}{3} )\), \((0\), \(0)\), and \((1\), \(-\frac{1}{3} )\) from \(log(f(x))\)
By computing the tangent lines at each point, finding the points of intersection, and merging the functions, we get the log of the envelope function, \[
log(e(x)) =
\begin{cases}
0 & \text{if $-\frac{2}{3} < x < \frac{2}{3}$} \\
\frac{2}{3} - |x| & \text{$otherwise$} \\
\end{cases}
\] Exponentiate to get the envelope function, \[
e(x) =
\begin{cases}
1 & \text{if $-\frac{2}{3} < x < \frac{2}{3}$} \\
e^{\frac{2}{3}} - |x| & \text{$otherwise$} \\
\end{cases}
\]
Let the log of the squeezing function be \(log(s(x)) = -\frac{|x|}{3}\) if \(-1 < x < 1\)
Exponentiate to get the squeezing function \(s(x) = e^{-\frac{|x|}{3}}\) if \(-1 < x < 1\)
logf <- function(x) {
-abs(x^3)/3
}
loge <- function(x) {
ifelse( (x>-2/3)&(x<2/3), 0, 2/3-abs(x) );
}
logs <- function(x) {
ifelse( (x>-1)&(x<1), -abs(x)/3, NA );
}
f <- function(x) {
exp(logf(x))
}
e <- function(x) {
exp(loge(x));
}
s <- function(x) {
exp(logs(x));
}
par(mfrow=c(1,2))
curve(logf(x), from = -2, to = 2, col = "blue")
curve(loge(x), add = T)
curve(logs(x), add = T, col = "red")
abline(v = -1, lty = 3,col = "red")
abline(v = 1, lty = 3,col = "red")
curve(f(x), from = -2, to = 2, col = "blue")
curve(e(x), add = T)
curve(s(x), add = T, col = "red")
abline(v = -1, lty = 3,col = "red")
abline(v = 1, lty = 3,col = "red")
\(\int_{-\infty}^{\infty} e(x) \; dx = \int_{-\infty}^{-\frac{2}{3}} e^{\frac{2}{3} + x } \; dx + \int_{-\frac{2}{3}}^{\frac{2}{3}} e^{0} \; dx + \int_{\frac{2}{3}}^{\infty} e^{\frac{2}{3} - x } \; dx = \frac{10}{3}\)
So the normalizing constant is \(\frac{3}{10}\)
Then the CDF of \(g\) is \[
G(x) =
\begin{cases}
\frac{3}{10}e^{\frac{2}{3}+x} & \text{if $x< -\frac{2}{3}$} \\
\frac{3}{10}x + \frac{1}{2} & \text{if $-\frac{2}{3} \leq x \leq \frac{2}{3}$} \\
1-\frac{3}{10}e^{\frac{2}{3}-x} & \text{if $x > \frac{2}{3}$} \\
\end{cases}
\]
Take the inverse of \(G(x)\), \[ G^{-1}(u) = \begin{cases} log(\frac{10}{3}u)-\frac{2}{3} & \text{if $0 < u < \frac{3}{10}$} \\ \frac{10}{3}(u-\frac{1}{2}) & \text{$-\frac{3}{10} \leq u \leq \frac{7}{10}$} \\ \frac{2}{3}-log(\frac{10}{3}(1-u)) & \text{$\frac{7}{10} < u < 1$} \\ \end{cases} \]
Ginv <- function(u) {
ifelse(u<3/10, log(u*10/3)-2/3, ifelse(u>7/10, 2/3-log((1-u)*10/3),(u-1/2)*10/3));
}
Below is a function for performing rejection sampling with a sample size of \(n\) points.
# adaptive rejection sampling function
ars <- function(n) {
x <- rep(NA, n);
# number of points accepted
ct <- 0;
# number of points sampled
total <- 0;
# number of points caught by squeeze
squeeze <- 0;
while(ct < n) {
y <- Ginv(runif(1));
u <- runif(1);
# check squeeze range
if(y > -1 && y < 1) {
# under squeeze
if(u < s(y)/e(y)) {
ct <- ct + 1;
x[ct] <- y;
squeeze <- squeeze + 1;
}
# above squeeze
else {
# under f
if(u < f(y)/e(y)) {
ct <- ct + 1;
x[ct]<-y;
}
}
}
# outside squeeze but under f
else if(u < f(y)/e(y)) {
ct <- ct + 1;
x[ct] <- y;
}
total <- total + 1;
}
list(x = x, acratio_sx = squeeze/total, acratio = ct/total);
}
Choose a sample size of 100,000. Below are a few points drawn using this method.
samp_size = 100000
set.seed(920)
ars_points <- ars(samp_size)
head(ars_points$x)
## [1] 0.2950836 1.3434853 1.1435056 -1.3004348 -0.3099414 0.1288534
The theoretical evelope ratio is \(\frac{\int_{-\infty}^{\infty} f(x)\;dx}{\int_{-\infty}^{\infty} e(x)\;dx}\), the proportion of points in \(f\) that are in \(e\).
integrate(f, lower = -Inf, upper = Inf)$value / integrate(e, lower = -Inf, upper = Inf)$value
## [1] 0.7727395
For this simulation, the envelope ratio is
ars_points$acratio
## [1] 0.7754643
The theoretical squeeze ratio is \(\frac{\int_{-1}^{1} s(x)\;dx}{\int_{-\infty}^{\infty} e(x)\;dx}\), the proportion of points in \(s\) that are in \(e\).
integrate(s, lower = -1, upper = 1)$value / integrate(e, lower = -Inf, upper = Inf)$value
## [1] 0.5102436
For this simulation, the squeeze ratio is
ars_points$acratio_sx
## [1] 0.512706
Now that we have our points from the sample, square each accepted x, and take the mean to get \(E[x^2]\).
mean(ars_points$x^2)
## [1] 0.7762001