acceptance-rejection method

The acceptance-rejection methodis an algorithm for generating random samples froman arbitrary probability distribution, givenas ingredients random samples from a related distribution and theuniform distribution.

The acceptance-rejection method’s chief advantage over the inverse CDFmethod of generating random numbersis that it requires neither the cumulative distribution functionnor its inverse to be computed. So in many cases it can run faster.

Set-up

• Let $X$ be a random variable with some other probability distributionthat we know how to draw samples from — that is, generate on a computer.

• Let $U$ be a random variable uniformly distributed on the interval $[0,1]$.

• Let $Y$ be the random variable that we want to be able to generate.Assume $Y$ has a probability distribution thatis absolutely continuous to the probability distribution for $X$,with density $\rho$.

• Further assume that the density $\rho$ is bounded above bya constant $c$.So $\rho (x)\le c$ for all $x$ in the range of $X$;and necessarily $c\ge 1$.

In most applications, both $X$and $Y$ will be continuous random variableswith densities $g$ and $f$ respectively.In that case we have $\rho (x)=f(x)/g(x)$,and we require $f(x)\le cg(x)$.

(If $g(x)=0$, then set $\rho (x)=0$.Note that we cannot have $f(x)>0$ and $g(x)=0$ simultaneouslyon a set of positive measure, since $Y$ has a distributionabsolutely continuous with respect to that of $X$.)

The random variables $X$ and $Y$ can be multi-variate.

Algorithm

The procedure to generate a value for $Y$ is:

1. Generate a value for $X$.

2. Generate a value for $U$.

3. If $U\le \rho (X)/c$, then set $Y=X$ (“accept”).

4. Otherwise, go back to step 1 (“reject”),repeating until we obtain a value of $Y$ in step 3.

Intuitive explanation

When we generate $X$ and $U$ as prescribed in the algorithm,we are in fact picking the point $(X,cU)$ in the rectangular boxbelow. And the test $U\le \rho (X)/c$ determinesthat point lies below the graph of $\rho$.It seems plausible that if we keep only the points thatfall under the graph of the density $\rho$, and ignore the points above,then the distribution of the abscissashould have density $\rho$.

Figure 1: Acceptance and rejection regions for a density

The acceptance-rejection method works more efficientlyas the distribution of $X$ and $Y$ become similar enough— that is, $\rho (x)$ and its upper bound $c$ are close to one.This makes the rejection region smaller, and so the algorithmis likely to go through fewer repetitions discarding the rejects.

Justification

We now prove that the acceptance-rejection method works.

Let $X_n$, for $n\in \mathbb{N}$, be independentrandom variables representing the samples, all with law $G$.Let $U_n$, for $n\in \mathbb{N}$,be independent random variables, with the uniform distribution over $[0,1]$,and independent from $X_n$.

Let

$$N=inf\{n\in \mathbb{N}:U_n\le \frac{\rho (X_n)}{c}\}$$

be the number of draws (for $X_n$ and $U_n$) taken by the algorithmbefore acceptance.Then we must show that $Y=X_N$ has the correct distribution:it should be distributed with density $\rho (x)$ with respect to $dG(x)$.

We have, by independence,

$$\mathbb{P}(N\ge n)=\mathbb{P}\left(\bigcap _{k=1}^{n-1}\{U_k>\frac{\rho (X_k)}{c}\}\right)=\mathbb{P}{(U_1>\frac{\rho (X_1)}{c})}^{n-1}.$$

We can calculate the last probability explicitly.Letting $H$ be the law of $Z_1=\rho (X_1)/c$,and using the independence of $U_1$ from $Z_1$, we find:

$\mathbb{P}(U_1>Z_1)$$={\displaystyle {\int }_0^1}{\displaystyle {\int }_z^1}𝑑u𝑑H(z)={\displaystyle {\int }_0^1}(1-z)𝑑H(z)$

$=1-𝔼Z_1$

$=1-{\displaystyle \int \frac{\rho (x)}{c}𝑑G(x)}=1-{\displaystyle \frac{1}{c}}.$

From the equation $N={\sum }_{n=1}^{\mathrm{\infty }}\mathrm{𝟏}\{N\ge n\}$,we take expectations of both sides, evaluatingthe resulting geometric series:

$$𝔼N=\sum _{n=1}^{\mathrm{\infty }}\mathbb{P}(N\ge n)=\sum _{n=1}^{\mathrm{\infty }}{(1-\frac{1}{c})}^{n-1}=c.$$

Thus almost surely,and the algorithm terminates, on average, after drawing $c$ samples.

Finally, for all measurable sets $B$,we have

$\mathbb{P}(Y\in B)$$={\displaystyle \sum _{n=1}^{\mathrm{\infty }}}\mathbb{P}(\{Y\in B\}\cap \{N=n\})$

$={\displaystyle \sum _{n=1}^{\mathrm{\infty }}}\mathbb{P}(X_n\in B,U_n\le {\displaystyle \frac{\rho (X_n)}{c}},N\ge n)$

$={\displaystyle \sum _{n=1}^{\mathrm{\infty }}}\mathbb{P}(N\ge n)\mathbb{P}(X_n\in B,U_n\le {\displaystyle \frac{\rho (X_n)}{c}})$

$={\displaystyle \sum _{n=1}^{\mathrm{\infty }}}\mathbb{P}(N\ge n){\displaystyle {\int }_B}{\displaystyle {\int }_0^{\rho (x)/c}}dudG(x)$

$=\left({\displaystyle {\int }_B}{\displaystyle \frac{\rho (x)}{c}}dG(x)\right){\displaystyle \sum _{n=1}^{\mathrm{\infty }}}\mathbb{P}(N\ge n)$

$={\displaystyle {\int }_B}\rho (x)𝑑G(x),$

as is to be shown.