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