Home > statistics > In search of a random gamma variate…

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

{{X}_{i}}\tilde{\ }Exp(\lambda )\Rightarrow \sum\limits_{i}{{{X}_{i}}}\tilde{\ }Ga(k,\lambda ).

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 k\in \mathbb{N}

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 Ga(\left\lfloor k \right\rfloor ,\lambda -1). The key to implementation is to maximise the ratio of the two densities, ie

\frac{f(x)}{g(x)}=\frac{{{e}^{x(-1+\lambda )-x\lambda }}{{x}^{-k+\alpha }}{{(\lambda -1)}^{-k}}{{\lambda }^{\alpha }}\Gamma (k)}{\Gamma (\alpha )}=\frac{{{e}^{-x}}{{x}^{-k+\alpha }}{{(\lambda -1)}^{-k}}{{\lambda }^{\alpha }}\Gamma (k)}{\Gamma (\alpha )}.

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 \alpha -k, then the actual value is M=\frac{f(\alpha -k)}{g(\alpha -k)}.

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!

About these ads
  1. Mike
    March 17, 2010 at 11:54 am

    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.

    • March 17, 2010 at 1:20 pm

      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.

      • Mike
        March 17, 2010 at 4:17 pm

        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.

  2. xi'an
    March 17, 2010 at 12:42 pm

    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

    • March 17, 2010 at 6:16 pm

      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).

  3. March 17, 2010 at 4:39 pm

    Mike :

    Also in the 2nd ‘lines’ function the closing bracket after “steelblue” should be after the sample in the density function.

    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.

  4. Mike
    March 17, 2010 at 6:07 pm

    Thanks Manos, it now works fine.

    Thanks again for taking the time to write the code.

  1. March 25, 2010 at 6:07 pm
  2. April 5, 2010 at 11:47 pm

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

Follow

Get every new post delivered to your Inbox.

%d bloggers like this: