In search of a random gamma variate…

March 16, 2010
By

(This article was first published on Stats raving mad » R, and kindly contributed to R-bloggers)

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 sumlimits_{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  kin 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(leftlfloor k rightrfloor ,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 )-xlambda }}{{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 $latex 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!

To leave a comment for the author, please follow the link and comment on his blog: Stats raving mad » 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.