## In search of a random gamma variate…

One of the most common exersices given to Statistical Computing,Simulation or relevant classes is the generation of random numbers from a gamma distribution. At first this might seem straightforward in terms of the lifesaving relation that exponential and gamma random variables share. So, it’s easy to get a gamma random variate using the fact that

.

The code to do this is the following

rexp1 <- function(lambda, n) { u <- runif(n) x <- -log(u)/lambda } rgamma1 <- function(k, lambda) { sum(rexp1(lambda, k)) }

This works unfortunately only for the case

In the general case we got to result to more “complex” (?) simulation, hence programming. The first technique we gonna use is rejection sampling. As the proposal (or proxy or instrumental) density we set the . The key to implementation is to maximise the ratio of the two densities, ie

.

We find the maximum of the ratio along the next lines.

m<-function(x) exp(-x)*x^(-k+a)*(1/(-1+r))^k*(r^a)*gamma(k)/gamma(a) grid<-seq(0,10,by=.1) a=3.45 k=3 r=2.06 plot(m(grid)~grid,type="l",col="red") grid() ind<-which.max(m(grid)) # 6 m.max<-grid[ind] # 0.5

Analytically we can work out that the maximum is achieved at , then the actual value is .

Now, we draw variates from the integer gamma until one is accepted.

*UPDATED*

n=200000 a=3.45 k=3 r=2.06 lambda=2.71 sample<-rep(NA,n) start<-Sys.time() for (i in 1:n) { # The following is a function tha draws ONE random variate. # It's useful to have it in this form one <- function(a, lambda) { k <- floor(a) m <-m(a-k) while (TRUE) { x <- rgamma1(k,lambda-1) p.accept <- dgamma(x,a,lambda)/(m*dgamma(x,k,lambda-1)) if (runif(1)<p.accept) return(x) } } sample[i]<-one(a, lambda) } Sys.time()-start # Time difference of 25.738 secs grid2 <- seq(0, 10, length.out=100) plot(grid2, m(a-k)*dgamma(grid2,k,lambda-1), type="l", lty=2, col="red", xlab="x", ylab="Density",lwd=2) lines(grid2, dgamma(grid2,a,lambda), col="green",lwd=2) lines(density(sample),col="steelblue",lwd=2) legend("topright", col=c("red","green","steelblue"),lty=c(2,1), legend=c("m(a-k)*g(x)","sample dansity","f(x)"))

Not bad!

The code looks really ineteresting but I have not been able to get it to work. The function m is defined but the function ‘one’ uses a variable m in the calculation of p.accept. The variable m is also used later in the plot command, but this variable is not defined. Should this be the upper case M?

I think also that the time calculation is wrong because start<-Sys.time() is placed after the start of the for loop but the calculation Sys.time()-start is after the close of the for loop.

You are right Mike. I posted the draft code i had somewhere in my computer insted of the one I put together for this post. i posted the correct one now.

I think there are still a few more bugs though! The variable m is still used in the plot function after the for loop. If this is replaced by m1 then m1 needs to be re-defined outside of the ‘one’ function. Also in the 2nd ‘lines’ function the closing bracket after “steelblue” should be after the sample in the density function.

Great derivation! This is also Exercise 2.17 in our book Introducing Monte Carlo Methods with R (Springer, 2009), see the solution manual on arxiv:1001.2906. Cheers, Christian

Thanx for the comment Christian!

You (and Casella) managed to put together a great book on Monte Carlo using R. I enjoyed it & it’s good to give the solution manual freely (opposed to traditional instructor limitations).

This shpuld be obvious to me at least! I should have noticed the absense of green ;)

Check the code now. I passed to Eclipse and found it OK.

Thanks Manos, it now works fine.

Thanks again for taking the time to write the code.