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.

About these ads

Leave a Comment »

No comments yet.

RSS feed for comments on this post. TrackBack URI

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

The Shocking Blue Green Theme Blog at WordPress.com.

Follow

Get every new post delivered to your Inbox.