Acceptance-rejection method

Introduced by J. von Neumann [a8], this is the most adaptable method for sampling from complicated distributions (cf. also Sample; Distribution). It works as follows.

Let $f ( x )$ be a given probability density (cf. also Density of a probability distribution), and let $h ( x )$ be a function such that $f ( x ) \leq h ( x )$ within the range of $f ( x )$. If the integral of $h ( x )$ over this range is a finite number $\alpha$, then $g ( x ) = h ( x ) / \alpha$ is a probability density function, and the following procedure is valid:

1) Take a random sample $X$ from the distribution with probability density $g ( x ) = h ( x ) / \alpha$.

2) Generate a uniform random deviate $U$ between zero and one. If $U \leq f ( X ) / h ( X )$, accept $X$ as a sample from the distribution $f ( x )$. Otherwise reject $X$ and go back to Step 1.

The ease of the method depends on the following properties of the hat function $h ( x )$:

A) One has to select a hat function $h ( x )$ from which it is easy to sample. Examples are given below.

B) The parameters of the hat function have to be determined in such a way that the area $\alpha$ below $h ( x )$ becomes minimal.

Optimal hat functions can be calculated by analytical methods. Assume that the hat function touches $f ( x )$ at two points $L$ (left) and $R$ (right), where $L < R$. Furthermore, suppose that the hat function depends on two parameters, called $m$ and $s$. Thus

\begin{equation} \tag{a1} f ( L ) = \alpha g ( L ; m , s ) , f ( R ) = \alpha g ( R ; m , s ), \end{equation}

and $f ( x ) \leq \alpha g ( x ; m , s )$ for all other $x$. Since $L$ and $R$ are local maxima of $f ( x ) / \operatorname { g } ( x ; m , s )$, one has the necessary conditions

\begin{equation} \tag{a2} \frac { f ^ { \prime } ( L ) } { f ( L ) } = \frac { g ^ { \prime } ( L ; m , s ) } { g ( L ; m , s ) } \end{equation}

and

\begin{equation} \tag{a3} \frac { f ^ { \prime } ( R ) } { f ( R ) } = \frac { g ^ { \prime } ( R ; m , s ) } { g ( R ; m , s ) }. \end{equation}

If $L$ and $R$ are uniquely determined, they should satisfy the sufficient conditions

\begin{equation*} \frac { f ^ { \prime } ( L ) } { f ( L ) } < \frac { g ^ { \prime } ( L ; m , s ) } { g ( L ; m , s ) } , \frac { f ^ { \prime } ( R ) } { f ( R ) } < \frac { g ^ { \prime } ( R ; m , s ) } { g ( R ; m , s ) }. \end{equation*}

Otherwise, the first derivative of $\operatorname { ln } ( f ( x ) / g ( x ; m , s ) )$ has to be discussed in detail.

Equations (a1), (a2) and (a3) are four equations for the determination of $L$, $R$, $m$, and $s$. Assuming that $L$, $R$ and $m$ can be expressed as functions of $s$, one has to minimize

\begin{equation} \tag{a4} \alpha ( s ) = \frac { f ( L ( s ) ) } { g ( L ( s ) ; m ( s ) , s ) } = \frac { f ( R ( s ) ) } { g ( R ( s ) ; m ( s ) , s ) }. \end{equation}

This leads to the necessary conditions

\begin{equation*} - \frac { d } { d s } \operatorname { ln } \alpha ( s ) = - \frac { d } { d L } \operatorname { ln } \frac { f ( L ) } { g ( L ; m , s ) } \frac { d L } { d s } + \end{equation*}

\begin{equation*} + \frac { d } { d m } \operatorname { ln } g ( L ; m , s ) \frac { d m } { d s } + \frac { d } { d s } \operatorname { ln } g ( L ; m , s ) = 0 , - \frac { d } { d s } \operatorname { ln } \alpha ( s ) = - \frac { d } { d R } \operatorname { ln } \frac { f ( R ) } { g ( R ; m , s ) } \frac { d R } { d s }+ \end{equation*}

\begin{equation*} + \frac { d } { d m } \operatorname { ln } g ( R ; m , s ) \frac { d m } { d s } + \frac { d } { d s } \operatorname { ln } g ( R ; m , s ) = 0. \end{equation*}

The first expression in each line is zero, by (a2) and (a3). Solving both equations for $dm / \ ds$, and comparing, yields the fundamental relation

\begin{equation} \tag{a5} \frac { d \operatorname { ln } g ( L ; m , s ) } { d m } \frac { d \operatorname { ln } g ( R ; m , s ) } { d s }= \end{equation}

\begin{equation*} = \frac { d \operatorname { ln } g ( R ; m , s ) } { d m } \frac { d \operatorname { ln } g ( L ; m , s ) } { d s }. \end{equation*}

(a1), (a2), (a3) and (a5) contain five conditions for finding candidates $L$, $R$, $m$, $s$ and $\alpha$.

Examples.

Triangular hat functions.

\begin{equation*} g ( x ; m , s ) = \left\{ \begin{array} { l l } { \frac { 1 } { s } - \frac { m - x } { s ^ { 2 } } } & { \text { if } m - s \leq x \leq m, } \\ { \frac { 1 } { s } - \frac { x - m } { s ^ { 2 } } } & { \text { if } m \leq x \leq m + s. } \end{array} \right. \end{equation*}

Samples may be obtained as $X \leftarrow m + s ( U _ { 1 } + U _ { 2 } - 1 )$. The fundamental identity (a5) leads to $s = R - L$. With its help all constants $L$, $R$, $m$, $s$, and $\alpha$ can be determined.

Double exponential hat functions.

The double exponential (or Laplace) distribution is given by

\begin{equation*} g ( x ; m , s ) = \left\{ \begin{array} { l l } { \frac { 1 } { 2 s } \operatorname { exp } ( \frac { x - m } { s } ) } & { \text { for } x \leq m, } \\ { \frac { 1 } { 2 s } \operatorname { exp } ( \frac { m - x } { s } ) } & { \text { for } x \geq m. } \end{array} \right. \end{equation*}

(See also Laplace distribution.)

Samples are obtained as $X \leftarrow m + T s E$, where $E$ is a standard exponential deviate and $T$ a random sign $\pm$. The fundamental identity (a5) leads to $2 s = R - L$.

Student-$t$ hat functions.

Its probability density function is given as

\begin{equation*} t _ { n } ( x ) = \frac { c _ { n } } { s } \left( 1 + \frac { ( x - m ) ^ { 2 } } { s ^ { 2 } n } \right) ^ { - ( n + 1 ) / 2 } \end{equation*}

for $- \infty < x < \infty$, $n \geq 1$, where

\begin{equation*} c _ { n } = \frac { 1 } { \sqrt { n } B \left( \frac { n } { 2 } , \frac { 1 } { 2 } \right) } = \frac { \Gamma \left( \frac { n + 1 } { 2 } \right) } { \sqrt { n \pi } \Gamma \left( \frac { n } { 2 } \right) }. \end{equation*}

(See also Student distribution.)

The $t$-family contains the Cauchy distribution for $n = 1$ and the normal distribution for $n \rightarrow \infty$ as extreme cases. For $n = 2$ and $n = 3$ special sampling methods are available: If $U$ denotes a $( 0,1 )$-uniform random variable, then $X \leftarrow ( U - 1 / 2 ) / ( \sqrt { ( U - U ^ { 2 } ) } / 2 )$ samples from $t_2$ and values from the $t_3$-distribution can be generated efficiently by the ratio-of-uniforms method of A.J. Kinderman and J.F. Monahan [a5], which will be discussed subsequently. The fundamental identity (a5) yields $s ^ { 2 } = ( R - m ) ( m - L )$. As before, $L$, $R$, $s$, $m$, and $\alpha$ can be determined explicitly.

The ratio-of-uniforms method was introduced by Kinderman and Monahan for sampling from a density $f ( x )$. First the table mountain-function is constructed: Let $a$, $b$ and $c$ be real numbers and let $k$ be equal to $1$, but $k = 2$ might be another possible choice for some densities. Then

\begin{equation*} y = \left\{ \begin{array} { l l } { \left( \frac { c } { a - x } \right) ^ { k + 1 } } & { \text { for } x \in ( - \infty , a - c ], } \\ { 1 } & { \text { for } x \in [ a - c , a - c + b ], } \\ { \left( \frac { b - c } { x - a } \right) ^ { k + 1 } } & { \text { for } x \in [ a - c + b , \infty ]. } \end{array} \right. \end{equation*}

Samples from the area below the table mountain-function are obtained by the transformation

\begin{equation*} X = \alpha + \frac { b V - c } { U ^ { 1 / k } } , Y = U ^ { 1 / k }. \end{equation*}

The table mountain-function is taken as a hat function for $f ( x ) / f$, where $f = \operatorname { max } f ( x )$. In this case the fundamental identity leads to $f ( L ) = f ( R )$ for calculating optimal constants.

Squeeze functions.

Step 2 of the acceptance-rejection method can be improved if some lower bound

\begin{equation*} b ( x ) \leq q ( x ) = \frac { f ( x ) } { h ( x ) } , \text { for all } - \infty < x < \infty, \end{equation*}

is known and easy to calculate. Step 2 then changes to:

2') Generate a uniform random deviate $U$ between zero and one. If $U \leq b ( X )$, accept $X$ as a sample from the target distribution $f ( x )$. If $U \geq f ( X ) / h ( X )$, reject $X$ and go back to $1$. Otherwise accept $X$.

Squeeze functions have been constructed in many procedures. See [a4] for some theory.

References

[a1]	J.H. Ahrens, U. Dieter, "Computer methods for sampling from gamma, beta, Poisson and binomial distributions" Computing , 12 (1974) pp. 223–246
[a2]	L. Devroye, "Non-uniform random variate generation" , Springer (1986)
[a3]	U. Dieter, "Optimal acceptance-rejection methods for sampling from various distributions" P.R. Nelson (ed.) , The Frontiers of Statistical Computation, Simulation, & Modeling , Amer. Ser. Math. Management Sci. (1987)
[a4]	U. Dieter, "Mathematical aspects of various methods for sampling from classical distributions" , Proc. 1989 Winter Simulation Conference (1989) pp. 477–483
[a5]	A.J. Kinderman, J.F. Monahan, "Computer generation of random variables using the ratio of uniform deviates" ACM Trans. Math. Software , 3 (1977) pp. 257–260
[a6]	D.E. Knuth, "The art of computer programming" , 2: Seminumerical algorithms , Addison-Wesley (1998) (Edition: Third)
[a7]	J. Leydold, "Automatic sampling with the ratio-of-uniform method" ACM Trans. Math. Software , 26 (2000) pp. 78–88
[a8]	J. von Neumann, "Various techniques used in connection with random digits. Monte Carlo methods" Nat. Bureau Standards , 12 (1951) pp. 36–38
[a9]	E. Stadlober, "Sampling from Poisson, binomial and hypergeometric distributions: Ratio of uniforms as a simple and fast alternative" Math.–Statist. Sekt. 303 Forschungsgesellschaft Joanneum, Graz, Austria (1989)
[a10]	E. Stadlober, "The ratio of uniforms approach for generating discrete random varates" J. Comput. Appl. Math. , 31 : 1 (1990) pp. 181–189