Why?

December 2, 2010

Random variable generation (Pt 2 of 3)

Filed under: AMCMC, R — Tags: , , , , — csgillespie @ 5:44 pm

Acceptance-rejection methods

This post is based on chapter 1.4 of Advanced Markov Chain Monte Carlo.

Another method of generating random variates from distributions is to use acceptance-rejection methods. Basically to generate a random number from f(x), we generate a RN from an envelope distribution g(x), where \sup f(x)/g(x) \le M \le \infty.  The acceptance-rejection algorithm is as follows:

Repeat until we generate a value from step 2:

1. Generate x from g(x) and U from Unif(0, 1)

2. If U \le \frac{f(x)}{M g(x)}, return x (as a random deviate from f(x)).

Example: the standard normal distribution

This example illustrates how we generate N(0, 1) RNs using the logistic distribution as an envelope distribution. First, note that

\displaystyle f(x) = \frac{1}{2\pi} e^{-x^2/2} \quad \mbox{and} \quad g(x) = \frac{e^{-x/s}}{s(1+ e^{-x/s})^2}

On setting s=0.648, we get M = 1.081. This method is fairly efficient and has an acceptance rate of

\displaystyle r = \frac{1}{M}\frac{\int f(x) dx}{\int g(x) dx} = \frac{1}{M} = 0.925

since both f and g are normalised densities.

R code

This example is straightforward to code:

myrnorm = function(M){
  while(1){
    u = runif(1); x = rlogis(1, scale = 0.648)
    if(u < dnorm(x)/M/dlogis(x, scale = 0.648))
      return(x)
  }
}

To check the results, we could call myrnorm a few thousand times:

hist(replicate(10000, myrnorm(1.1)), freq=FALSE)
lines(seq(-3, 3, 0.01), dnorm(seq(-3, 3, 0.01)), col=2)

Example: the standard normal distribution with a squeeze

Suppose the density f(x) is expensive to evaluate. In this scenario we can employ an easy to compute function s(x), where 0 \le s(x) \le g(x) . s(x) is called a squeeze function. In this example, we’ll use a simple rectangular function, where s(x) = 0.22 for x=-1, \ldots, 1. This is shown in the following figure:


The modified algorithm is as follows:

Repeat until we generate a value from step 2:

1. Generate x from g(x) and U from Unif(0, 1)

2. If U \le \frac{s(x)}{M g(x)} or U \le \frac{f(x)}{M g(x)}, return x (as a random deviate from f(x)).

Hence, when U \le \frac{s(x)}{M g(x)} we don’t have to compute f(x). Obviously, in this example f(x) isn’t that difficult to compute.

The Shocking Blue Green Theme. Blog at WordPress.com.

Follow

Get every new post delivered to your Inbox.