**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 $latex 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)"))

Not bad!

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