Doug Bates’s Gibbs Sampler Example

This is an example from Dr. Doug Bates’s slides Julia for R Programmers.

The task is to create a Gibbs sampler for the density \[ f(x, y) = k x^2 exp(-xy^2 - y^2 + 2y - 4x), x > 0 \] using the conditional distributions \[ \begin{eqnarray*} X | Y &\sim& \Gamma \left( 3, \frac{1}{y^2 + 4} \right) \\ Y | X &\sim& N \left( \frac{1}{1+x}, \frac{1}{2(1+x)} \right). \end{eqnarray*} \]

This is an R function for the simple Gibbs sampler

library(Matrix)
Rgibbs <- function(N, thin) {
  mat <- matrix(0, nrow=N, ncol=2)
  x <- y <- 0
  for (i in 1:N) {
    for (j in 1:thin) {
      x <- rgamma(1, 3, y * y + 4) # 3rd arg is rate
      y <- rnorm(1, 1 / (x + 1), 1 / sqrt(2 * (x + 1)))
    }
    mat[i,] <- c(x, y)
  }
  mat
}

Generate a bivariate sample of size 10,000 with a thinning of 500. How long does it take?

system.time(Rgibbs(10000, 500))

##    user  system elapsed 
##  40.243   2.385  42.638

Does compiler package help?

library(compiler)
cmpRgibbs <- cmpfun(Rgibbs)
system.time(cmpRgibbs(10000, 500))

##    user  system elapsed 
##  37.261   0.048  37.315

Show system information and clean up.

sessionInfo()

## R version 3.1.1 (2014-07-10)
## Platform: x86_64-apple-darwin13.1.0 (64-bit)
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] compiler  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
## [1] Matrix_1.1-4
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.8    evaluate_0.5.5  formatR_1.0     grid_3.1.1     
##  [5] htmltools_0.2.6 knitr_1.8       lattice_0.20-29 rmarkdown_0.4.2
##  [9] stringr_0.6.2   tools_3.1.1     yaml_2.1.13