# 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

.

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 $alpha -k$, 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)"))
```

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, 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: , , , , , , , ,

## Recent popular posts

Contact us if you wish to help support R-bloggers, and place your banner here.

# Never miss an update! Subscribe to R-bloggers to receive e-mails with the latest R posts.(You will not see this message again.)

Click here to close (This popup will not appear again)