--- title: "JAGS Lab:" author: 'Name - here' output: html_document --- # Lab Overview JAGS (Just Another Gibbs Sampler) is a system that builds MCMC samplers. For a comprehensive overview of JAGS, see the user manual [http://people.stat.sc.edu/hansont/stat740/jags_user_manual.pdf](http://people.stat.sc.edu/hansont/stat740/jags_user_manual.pdf). The goal today is to install JAGS and learn a little about the basic functionality by executing code from R. ### 1. Installing JAGS 1. Download and install JAGS: [https://sourceforge.net/projects/mcmc-jags/files/](https://sourceforge.net/projects/mcmc-jags/files/) 2. In R, install the `rjags` package. ### 2. JAGS Overview JAGS is designed for inference on Bayesian models using Markov Chain Monte Carlo (MCMC) simulation. Running a model refers to generating samples from the posterior distribution of the model parameters. This takes place in five steps: 1. Definition of the model 2. Compilation 3. Initialization 4. Adaptation and burn-in 5. Monitoring The next stages of analysis are done outside of JAGS: convergence diagnostics, model criticism, and summarizing the samples must be done using other packages more suited to this task. We will use the `coda` package to assess MCMC samples. ### 3. Tutorial #### 3a Beta-Binomial Model ```{r} ## JAGS lab - STAT 532 require(rjags) # Must have previously installed package rjags. # Simulate the data from binomial distribution: set.seed(09282017) n <- 50 theta.true <- .3 y.sim <- rbinom(n, 1, theta.true) dataList = list( # Put the information into a list. y = y.sim , Ntotal = n ) # Define the model as a string object: modelString = " model { for ( i in 1:Ntotal ) { y[i] ~ dbern( theta ) } theta ~ dbeta( 1 , 1 ) } " writeLines( modelString , con="TEMPmodel.txt" ) #creates a text file object # Initialize the chains based on MLE of data. # Option: Use function that generates random values for each chain: initsList = function() { resampledY = sample( y.sim , replace=TRUE ) thetaInit = sum(resampledY)/length(resampledY) thetaInit = 0.001+0.998*thetaInit # keep away from 0,1 return( list( theta=thetaInit ) ) } initsList = function() { return( list( theta=runif(1) ) ) } initsList() ``` ##### Q1. Summarize what the `initsList` function is doing. ```{r} # Run the chains: jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList , inits=initsList , n.chains=3 , n.adapt=500 ) ``` ##### Q2. Describe the five arguments in the `jags.model()` function. ```{r} codaSamples = coda.samples( jagsModel , variable.names=c("theta") , n.iter=5000) summary(codaSamples) ``` ##### Q3. What is the coda.samples function doing in this case? ```{r} plot(codaSamples) traceplot(codaSamples) ``` ##### Q4. What do the three colors on the traceplot represent? ```{r} gelman.plot(codaSamples) HPDinterval(codaSamples) ``` #### 3b Normal Model ```{r} ### Normal Model N <- 1000 mu <- 5 tau.sq <- 1/5 # precision = 1 / sigma^2 x <- rnorm(N, mu, sqrt(1/tau.sq)) # Define the model: modelString = " model { for (i in 1:N) { x[i] ~ dnorm(mu, tau.sq) } mu ~ dnorm(0, .0001) tau.sq ~ dgamma(.1, .1) } " writeLines( modelString , con="TEMPmodel.txt" ) ``` ##### Q5. Write out the model and priors specified in this case. Note the dollar sign can be used to embed latex commands in line, $p(\mu) \sim$. Two dollar signs are used to signify equations. $$p(\tau^2) \sim $$ ```{r} jags <- jags.model('TEMPmodel.txt', data = list('x' = x, 'N' = N), n.chains = 4, n.adapt = 100) codaSamples = coda.samples( jags , variable.names=c("mu","tau.sq") , n.iter=1000) summary(codaSamples) plot(codaSamples) HPDinterval(codaSamples) traceplot(codaSamples) ``` ##### Q6. Summarize the results of the MCMC sampler.