Acceptance-Rejection Method

• In this section we discuss another method for the generation of pseudo-random numbers $y_1,y_2,\ldots$
  • that can be regarded as realizations of independent and identically distributed random variables $Y_1,Y_2\ldots$. Their distribution function is assumed to be given; it is denoted by $G$.
  • This method also requires a sequence of independent and identically distributed pseudo-random numbers $x_1,x_2,\ldots$, but we abandon the condition that they need to be uniformly distributed on $(0,1]$.
  • The only condition we impose on their distribution function $F$ is that $G$ needs to be absolutely continuous with respect to $F$ with bounded density $g(x)=dG(x)/dF(x)$,
  • i.e., for some constant $c>0$, we have

    $$g(x)\le c \qquad G(y)=\int_{-\infty}^y g(x)\,dF(x)\,,\qquad \forall\, x,y\in\mathbb{R}\,.$$ (19)

    
    

• First of all we consider the discrete case.
  • Let $a_j=j$ for all $j=0,1,\ldots$, and let ${\mathbf{p}}=(p_0,p_1,\ldots)^\top$ and ${\mathbf{q}}=(q_0,q_1,\ldots)^\top$ be two arbitrary probability functions such that for all $j=0,1,\ldots$ $p_j=0$ implies $q_j=0$.
  • Let $X:\Omega\to\{0,1,\ldots\}$ be a random variable $P(X=j)=p_j$ for all $j=0,1,\ldots$,
  • and let $c>0$ be a positive number

    $$\Bigl(\;g(j)=\;\Bigr)\qquad\frac{q_j}{p_j}\le c \mbox{for all $j\ge 0$\ such that $p_j>0$.}$$ (20)

    
    

Theorem 3.5    

• Let $(U_1,X_1),(U_2,X_2),\ldots$ be a sequence of independent and identically distributed random vectors whose components are independent. Furthermore, let $U_i$ be a $(0,1]$-uniformly distributed random variable and $X_i$ be distributed according to ${\mathbf{p}}$.
• Then
  • the random variable

    $$I=\min\Bigl\{k\ge 1: U_k<\frac{q_{X_k}}{c\,p_{X_k}}\Bigr\}$$ (21)

    
    
    is geometrically distributed with expectation $c$, i.e., $I\sim\,{\rm Geo}(c^{-1})$,
  • and the random variable $Y=X_I$ is distributed according to ${\mathbf{q}}$.

Proof

 

• By the definition of $I$ given in (21), we obtain for all $j\ge 1$
  

  $$P(I=j) = P\Bigl(U_1\ge\frac{q_{X_1}}{c\,p_{X_1}}\,,\ldots,\, U_{j-1}\ge\frac{q_{X_{j-1}}}{c\,p_{X_{j-1}}}\,,\,U_j<\frac{q_{X_j}}{c\,p_{X_j}}\Bigr)$$

  $$= P\Bigl(U_1\ge\frac{q_{X_1}}{c\,p_{X_1}}\Bigr)\,\ldots\, P\Bigl(U_... ..._{X_{j-1}}}{c\,p_{X_{j-1}}}\Bigr)\, P\Bigl(U_j<\frac{q_{X_j}}{c\,p_{X_j}}\Bigr)$$

  $$= p\,q^{j-1}\,,$$

  
  
  • where $q=1-p$ and
    

    $$p = P\Bigl(U_1<\frac{q_{X_1}}{c\,p_{X_1}}\Bigr)$$

    $$= \sum\limits_{k:\, p_k>0} P\Bigl(U_1<\frac{q_{X_1}}{c\,p_{X_1}}\mid X_1=k\Bigr)\,P(X_1=k)$$

    $$= \sum\limits_{k:\, p_k>0} P\Bigl(U_1<\frac{q_k}{c\,p_k}\Bigr)\, p_k$$

    $$= \sum\limits_{k:\, p_k>0}\frac{q_k}{c\,p_k}\; p_k = \;\frac{1}{c}\;.$$

    
    

  • This shows $I\sim\,{\rm Geo}(c^{-1})$.
• Furthermore, for any $j\ge 1$ such that $p_j>0$, we get that
  

  $$P(Y=j) = P(X_I=j) \;=\; \sum\limits_{k=1}^\infty P(X_I=j,\, I=k)$$

  $$= \sum\limits_{k=1}^\infty P(X_k=j,\,I=k) \;=\; \sum\limits_{k=1}^\infty P(X_k=j)\, q^{k-1}\,P\Bigl(U_k\le\;\frac{q_j}{cp_j}\Bigr)$$

  $$= \sum\limits_{k=1}^\infty p_j\,q^{k-1}\;\frac{q_j}{cp_j} \;=\;\frac{q_j}{c}\; \sum\limits_{k=1}^\infty q^{k-1}$$

  $$= \frac{q_j}{c}\;\frac{1}{1-q}\; \;=\; q_j\,,$$

  
  
  and for all $j\ge 1$ such that $p_j=0$
  $$P(Y=j) \;=\; \sum\limits_{k=1}^\infty P(X_k=j,\, I=k)\;\le\; \sum\limits_{k=1}^\infty P(X_k=j)\;=\; 0\,.$$
  
   
    $\Box$
  
  
  

Remarks

 

• Theorem 3.5 implies that the mean number of $F$-distributed pseudo-random numbers necessary to obtain a $G$-distributed random number is $c$.
• In case there are several alternatives for the choice of the the distribution function $F$,
  • possessing equally nice properties with respect to the generation of $F$-distributed pseudo-random numbers,
  • then one should choose the distribution function with the smallest $c$.
• Furthermore, as a consequence of Theorem 3.5,
  • the values $g(x)$ and $g(j)$ of the density in (19) and (20), respectively need only be known up to a constant factor.

In the general (i.e. not necessarily discrete) case one can proceed in a similar way. The following result will serve as foundation for constructing acceptance-rejection algorithms.

Theorem 3.6    

• Let $F,G:\mathbb{R}\to[0,1]$ be two arbitrary distribution functions such that $(\ref{def.dic.geh})$ holds.
• Let $(U_1,X_1),(U_2,X_2),\ldots$ be a sequence of independent and identically distributed random vectors whose components are independent. Furthermore, let $U_i$ be a $(0,1]$-uniformly distributed random variable and $X_i$ be distributed according to $F$.
• Then the random variable

  $$I=\min\Bigl\{k\ge 1: U_k<\frac{g(X_k)}{c}\Bigr\}$$ (22)

  
  
  is geometrically distributed with expectation $c$, i.e., $I\sim\,{\rm Geo}(c^{-1})$ and the random variable $Y=X_I$ is distributed according to $G$.

Proof

 

• Similarly to the proof of Theorem 3.5 we obtain $P(I=j)=p\,q^{j-1}$ for any $j\ge 1$ where
  

  $$p = P\Bigl(U_1<\frac{g(X_1)}{c}\Bigr)= \int_\mathbb{R}P\Bigl(U_1<\frac{g(X_1)}{c}\mid X_1=x\Bigr)\,dF(x)$$

  $$= \int_\mathbb{R}P\Bigl(U_1<\frac{g(x)}{c}\Bigr)\,dF(x)= \int_\mathbb{R}\;\frac{g(x)}{c}\;dF(x)=\frac{1}{c}\;.$$

  
  

• Furthermore, for all $y\in\mathbb{R}$ we have
  

  $$P(Y\le y) = P(X_I\le y)\;=\; \sum\limits_{k=1}^\infty P(X_I\le y,\, I=k) = \sum\limits_{k=1}^\infty P(X_k\le y,\, I=k)$$

  $$= \sum\limits_{k=1}^\infty \int_{-\infty}^y P(I=k\mid X_k=v)\, dF(v... ..._{k=1}^\infty q^{k-1}\int_{-\infty}^y P\Bigl(U_k<\;\frac{g(v)}{c}\Bigr)\, dF(v)$$

  $$= \frac{1}{1-q}\;\int_{-\infty}^y \;\frac{g(v)}{c}\; dF(v) = \int_{-\infty}^y g(v) dF(v)\;=\; G(y)\,,$$

  
  
  where $1-q=p=c^{-1}$.
  
  $\Box$

In the same way we obtain the following vectorial version of Theorem 3.6.

Theorem 3.7    

• Let $m\ge 1$ be an arbitrary but fixed natural number and let $F,G:\mathbb{R}^m\to[0,1]$ be two arbitrary distribution functions (of $m$-dimensional random vectors) and let $c>0$ be a constant such that

  $$g({\mathbf{x}})\le c \qquad G({\mathbf{y}})=\int_{(-\infty,{\mathbf{y}}]} g({\mathbf{x... ...dF({\mathbf{x}})\,,\qquad \forall\, {\mathbf{x}},{\mathbf{y}}\in\mathbb{R}^m\,.$$ (23)

  
  

• Let $(U_1,{\mathbf{X}}_1),(U_2,{\mathbf{X}}_2),\ldots$ be a sequence of independent and identically distributed random vectors whose components are also independent. Furthermore, let $U_i$ be a $(0,1]$-uniformly distributed random variable and ${\mathbf{X}}_i$ be distributed according to $F$.
• Then the random variable

  $$I=\min\Bigl\{k\ge 1: U_k<\frac{g({\mathbf{X}}_k)}{c}\Bigr\}$$ (24)

  
  
  is geometrically distributed with expectation $c$, i.e., $I\sim\,{\rm Geo}(c^{-1})$ and the random vector ${\mathbf{Y}}={\mathbf{X}}_I$ is distributed according to $G$.

Examples

 

1. Uniform distribution on bounded Borel sets 
   • Let the random vector ${\mathbf{X}}:\Omega\to\mathbb{R}^m$ (with distribution function $F$) be uniformly distributed on the square $(-1,1]^m$ and let $B\in\mathcal{B}((-1,1]^m$ be an arbitrary Borel subset of $(-1,1]^m$ of positive Lebesgue measure $\vert B\vert$.
   • Then the distribution function $G:\mathbb{R}^m\to[0,1]$ given by
     $$G({\mathbf{y}})=\int_{(-\infty,{\mathbf{y}}]}\;\frac{{1\hspace{-1... ...}{\vert B\vert}\;dF({\mathbf{x}})\,,\qquad\forall\,{\mathbf{y}}\in\mathbb{R}^m$$
     is absolutely continuous with respect to $F$ and we obtain for the (Radon-Nikodym) density
     $g:\mathbb{R}^m\to[0,\infty)$ that
     $$g({\mathbf{x}})=\frac{{1\hspace{-1mm}{\rm I}}({\mathbf{x}}\in B)}{\vert B\vert}\le c=\vert B\vert^{-1}$$   and$$\qquad \frac{g({\mathbf{x}})}{c}\;={1\hspace{-1mm}{\rm I}}({\mathbf{x}}\in B)\,, \qquad\forall\,{\mathbf{x}}\in\mathbb{R}^m\,.$$

   • By Theorem 3.7 we can now in the following way generate pseudo-random vectors ${\mathbf{y}}_1,{\mathbf{y}}_2,\ldots$ that are uniformly distributed on $B$.

     1.

     Generate $m$ pseudo-random numbers $u_1,\ldots,u_m$ that are uniformly distributed on the interval $(0,1]$.

     2.

     If $(2u_1-1,\ldots,2u_m-1)^\top\not\in B$, then return to step 1.

     3.

     Otherwise put ${\mathbf{y}}=(2u_1-1,\ldots,2u_m-1)^\top$.

   
   

2. Normal distribution 
   • As an alternative to the Box-Muller algorithm discussed in Section 3.2.1 we will now introduce another method to generate normally distributed pseudo-random numbers,
     • which is often called the polar method.
     • Notice that the polar method avoids calculating the trigonometric functions in (12).
   • Let the random vector $(V_1,V_2)$ be uniformly distributed on the unit circle $B$, where
     $B=\{(x_1,x_2)\in\mathbb{R}^2:\,x_1^2+x_2^2\le 1\}$.
   • Then, the random vector $(Y_1,Y_2)$ where
     $$Y_1=\sqrt{-2 \log(V_1^2+V_2^2)}\; \frac{V_1}{\sqrt{V_1^2+V_2^2}}\,,\qquad Y_2=\sqrt{-2 \log(V_1^2+V_2^2)}\; \frac{V_2}{\sqrt{V_1^2+V_2^2}}$$
     is N $({\mathbf{o}},{\mathbf{I}})$-distributed, i.e., $Y_1$, $Y_2$ are independent and N$(0,1)$-distributed random variables. This can be seen as follows.
     • By the substitution
       $$v_1=r\cos \theta\,,\qquad v_2=r\sin \theta\,,$$

     • i.e. by a transformation into polar coordinates we obtain for arbitrary $y_1,y_2\in\mathbb{R}$
       

       $${P(Y_1\le y_1,\, Y_2\le y_2)}$$

       $$= \frac{1}{\pi}\;\int_B {1\hspace{-1mm}{\rm I}}\Bigl(\frac{v_1\sqrt... ...v_2 \sqrt{-2 \log(v_1^2+v_2^2)}}{\sqrt{v_1^2+v_2^2}}\,\le y_2\Bigr)\,d(v_1,v_2)$$

       $$= \frac{1}{\pi}\;\int_0^{2\pi}\int_0^1 r\,{1\hspace{-1mm}{\rm I}}\B... ...\theta\,\le y_1\,,\, \sqrt{-2 \log(r^2)}\; \sin\theta\le y_2\Bigr)\,dr\,d\theta$$

       $$= \frac{1}{2}\;\frac{1}{2\pi}\;\int_0^{2\pi}\int_0^\infty {1\hspace... ...}\cos\theta\le y_1,\, \sqrt{x}\sin\theta\le y_2\Bigr)\;e^{-x/2}\,dx\,d\theta\,,$$

       
       

     • where the last equality results from the following substitution:
       $$x=-2 \log(r^2)$$   bzw.$$\qquad -\;\frac{1}{2}\; e^{-x/2}\,dx=2r\,dr\,.$$

     • By the same argument that was used to verify formula (12) in Section 3.2.1 one can check that the last term can be written as the product $F(y_1)F(y_2)$ of two N$(0,1)$-distribution functions.
   • The pseudo-random numbers $y_1,\ldots,y_{2n}$ with
     $$y_{2k-1}=\sqrt{-2 \log(v_{2k-1}^2+v_{2k}^2)}\; \frac{v_{2k-1}}{\s... ...sqrt{-2 \log(v_{2k-1}^2+v_{2k}^2)}\; \frac{v_{2k}}{\sqrt{v_{2k-1}^2+v_{2k}^2}}$$
     • can thus be regarded as realizations of independent and N$(0,1)$-distributed random variables,
     • if $(v_1,v_2),\ldots,(v_{2n-1},v_{2n})$ are realizations of the random variables $(V_1,V_2),\ldots,(V_{2n-1},V_{2n})$ that are independent and uniformly distributed on the unit circle
       $$B=\{(x_1,x_2)\in\mathbb{R}^2:\,x_1^2+x_2^2\le 1\}\,.$$

     • Those can be generated via acceptance-rejection sampling as explained in the last example.