Biostat M280: Doug Bates’s Gibbs Sampler Example Using R

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.146   1.452  41.884

Does compiler package help?

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

##    user  system elapsed 
##  37.497   1.463  39.055

Show system information and clean up.

sessionInfo()

## R version 3.2.3 (2015-12-10)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X 10.11.2 (El Capitan)
## 
## 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.2-3
## 
## loaded via a namespace (and not attached):
##  [1] tools_3.2.3     htmltools_0.2.6 yaml_2.1.13     rmarkdown_0.9.2
##  [5] grid_3.2.3      knitr_1.11      stringr_0.6.2   digest_0.6.8   
##  [9] lattice_0.20-33 evaluate_0.8