---
title: "Rejection sampling"
author: "Ege Rubak and Jesper Møller"
date: ""
output: slidy_presentation
---

## Rejection sampling (or von Neumann sampling or acceptance-rejection method)

Rejection sampling is a general method for simulation of a random variable (or random vector) $X$ with a pdf $f\propto f_0$, assuming there is another random variable (or random vector of the same dimension) $Y$ with pdf $g\propto g_0$ so that $f_0(x)\le Mg_0(x)$ for all possible realizations $x$ of $X$, where $M$ is a positive number. 

Rejection sampling requires only that we need to know $g_0$ and $M$, and it works as follows.

Repeat generating independently

- realizations $Y=y$ and $U=u$ until $u\le f_0(y)/[Mg_0(y)]$,
- and then return $X=y$ as a simulation from $f$.

The idea is that if it is hard to simulate directly from $f$, it should be simpler to simulate from $g$. 

The algorithm works best if $f$ and $g_0$ are approximately proportional. 

If $f=c_1f_0$ and $g=c_2g_0$, where $c_1$ and $c_2$ are normalizing constants, then the probability of getting acceptance is
$$P(U\le f(Y)/[Mg_0(Y)]) = \dots = c_2/(c_1M),$$
so the best choice of $M$ is $M=c_2/c_1$ (whose value may be unknown). 

## Target density

Assume that the target density $f$ has support on $I=[-4,4]$ where it is
proportional to
$$
f_0(x; a,b) = \exp(a (x - a)^2 - b x^4),
$$
where $a,b$ are known parameters (we will use $a=0.4$ and $b=0.08$).

```{r}
f0 <- function(x, a=.4, b=.08){exp(a * (x - a)^2 - b * x^4)}
curve(f0(x), from = -4.5, to = 4.5)
```

## Proposal distribution

Assume we make uniform proposals on $I=[-4,4]$, i.e. $g(x)=1/8$ on $I$ and zero
outside $I$. Note that $f_0(x) \le 3.1$ on $I$, so we can use $M = 3.1/(1/8) = 24.8$:
```{r}
curve(f0(x), from = -4.5, to = 4.5)
M <- 24.8
g <- function(x){dunif(x, -4, 4)}
curve(M*g(x), add = TRUE, col = "red")
```


## Running the algorithm

Now we can generate proposals from $g$ (`runif()`) and note whether they are accepted as realizations from the target distribution or not:
```{r}
N <- 10000
y <- runif(N, -4, 4)
p_accept <- f0(y)/(M*g(y))
u <- runif(N, 0, 1)
keep <- u<p_accept
head(round(cbind(y, percent = 100*p_accept, keep), 3))
```

## Results

The resulting acceptance rate is
```{r}
mean(keep)
```

A histogram of the samples:
```{r}
hist(y[keep], prob = TRUE, col = "gray")
norm_const <- integrate(f0, -4, 4)$value
curve(f0(x)/norm_const, col = "red", add = TRUE)
```