Random variable generation (Pt 2 of 3)

December 2, 2010
By

(This article was first published on Why? » R, and kindly contributed to R-bloggers)

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}{2pi} 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.


To leave a comment for the author, please follow the link and comment on his blog: Why? » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Tags: , , , ,

Comments are closed.