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.

November 27, 2010

Advanced Markov Chain Monte Carlo Methods (AMCMC)

Filed under: AMCMC, R — Tags: , , — csgillespie @ 4:55 pm

I’ve just received my copy of Advanced Markov Chain Monte Carlo Methods, by Liang, Liu, & Carroll. Although my PhD didn’t really involve any Bayesian methodology (and my undergrad was devoid of any Bayesian influence), I’ve found that the sort of problems I’m now tackling in systems biology demand a Bayesian/MCMC approach. There are a number of introductory books on MCMC, but not that many on advanced techniques.

This book suggests that it could be used as a possible textbook or reference guide in a one-semester statistics graduate course. I’m not entirely convinced that it would be a good textbook, but as a reference it looks very promising. A word of warning, if you’re not familiar with MCMC then you should try the Dani Gamerman MCMC book first. Some later chapters look particularly interesting, including topics on adaptive proposals, population-based MCMC methods and dynamic weightings.

Anyway, I intend to work through the book (well at least a few of the chapters) and post my results/code as I go. Well that’s the plan anyway.

 

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

Follow

Get every new post delivered to your Inbox.