View source for Acceptance-rejection method

<!--This article has been texified automatically. Since there was no Nroff source code for this article, 
the semi-automatic procedure described at https://encyclopediaofmath.org/wiki/User:Maximilian_Janisch/latexlist
was used.
If the TeX and formula formatting is correct, please remove this message and the {{TEX|semi-auto}} category.

Out of 112 formulas, 112 were replaced by TEX code.-->

{{TEX|semi-auto}}{{TEX|done}}
Introduced by J. von Neumann [[#References|[a8]]], this is the most adaptable method for sampling from complicated distributions (cf. also [[Sample|Sample]]; [[Distribution|Distribution]]). It works as follows.

Let $f ( x )$ be a given probability density (cf. also [[Density of a probability distribution|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 } } } &amp; { \text { if } m - s \leq x \leq m, } \\ { \frac { 1 } { s } - \frac { x - m } { s ^ { 2 } } } &amp; { \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 } ) } &amp; { \text { for } x \leq m, } \\ { \frac { 1 } { 2 s } \operatorname { exp } ( \frac { m - x } { s } ) } &amp; { \text { for } x \geq m. } \end{array} \right. \end{equation*}

(See also [[Laplace distribution|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|Student distribution]].)

The $t$-family contains the [[Cauchy distribution|Cauchy distribution]] for $n = 1$ and the [[Normal distribution|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 [[#References|[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 } } &amp; { \text { for } x \in ( - \infty , a - c ], } \\ { 1 } &amp; { \text { for } x \in [ a - c , a - c + b ], } \\ { \left( \frac { b - c } { x - a } \right) ^ { k + 1 } } &amp; { \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 [[#References|[a4]]] for some theory.

====References====
<table><tr><td valign="top">[a1]</td> <td valign="top">  J.H. Ahrens,   U. Dieter,   "Computer methods for sampling from gamma, beta, Poisson and binomial distributions"  ''Computing'' , '''12'''  (1974)  pp. 223–246</td></tr><tr><td valign="top">[a2]</td> <td valign="top">  L. Devroye,   "Non-uniform random variate generation" , Springer  (1986)</td></tr><tr><td valign="top">[a3]</td> <td valign="top">  U. Dieter,   "Optimal acceptance-rejection methods for sampling from various distributions"  P.R. Nelson (ed.) , ''The Frontiers of Statistical Computation, Simulation, &amp; Modeling'' , ''Amer. Ser. Math. Management Sci.''  (1987)</td></tr><tr><td valign="top">[a4]</td> <td valign="top">  U. Dieter,   "Mathematical aspects of various methods for sampling from classical distributions" , ''Proc. 1989 Winter Simulation Conference''  (1989)  pp. 477–483</td></tr><tr><td valign="top">[a5]</td> <td valign="top">  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</td></tr><tr><td valign="top">[a6]</td> <td valign="top">  D.E. Knuth,   "The art of computer programming" , '''2: Seminumerical algorithms''' , Addison-Wesley  (1998)  (Edition: Third)</td></tr><tr><td valign="top">[a7]</td> <td valign="top">  J. Leydold,   "Automatic sampling with the ratio-of-uniform method"  ''ACM Trans. Math. Software'' , '''26'''  (2000)  pp. 78–88</td></tr><tr><td valign="top">[a8]</td> <td valign="top">  J. von Neumann,   "Various techniques used in connection with random digits. Monte Carlo methods"  ''Nat. Bureau Standards'' , '''12'''  (1951)  pp. 36–38</td></tr><tr><td valign="top">[a9]</td> <td valign="top">  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)</td></tr><tr><td valign="top">[a10]</td> <td valign="top">  E. Stadlober,   "The ratio of uniforms approach for generating discrete random varates"  ''J. Comput. Appl. Math.'' , '''31''' :  1  (1990)  pp. 181–189</td></tr></table>