MCMC software

class: center, middle, inverse, title-slide

# MCMC software
### Math 315: Bayesian Statistics

---

# Summary of MCMC

.large[
- With the combination of Gibbs and Metropolis-Hastings
sampling we can fit virtually any model

- In some cases Bayesian computing is actually preferable
to maximum likelihood analysis

- In most cases Bayesian computing is slower

- However, it is often worth the wait for
improved uncertainty quantification and interpretability

- In all cases it is important to carefully monitor convergence 
]

.footnote[Slide adapted from Brian Reich]

---

# Why Just Another Gibbs Sampler (JAGS)?

.large[
- You can fit virtually any model

- You can call JAGS from R, allowing for plotting and
data manipulation in R

- It runs on all platforms: LINUX, Mac, Windows

- There is a lot of help online

- R has many built in packages for convergence diagnostics
]

.footnote[Slide adapted from Brian Reich]

---

# How does JAGS work?

.large[
- You specify the model by declaring the likelihood and priors

- JAGS then sets up the MCMC sampler
    + e.g., works out the full conditional distributions for all parameters

- It returns MCMC samples in a matrix or array

- It also automatically produces posterior summaries like
means, credible sets, and convergence diagnostics

- User’s manual: http://people.stat.sc.edu/hansont/stat740/jags_user_manual.pdf
]

.footnote[Slide adapted from Brian Reich]

---

# Steps for any analysis using JAGS

.large[
1. Specify the model as a string

2. Compile the model using `jags.model()`

3. Draw burn-in samples using `update()`

4. Draw posterior samples using `coda.samples()`

5. Inspect/summarize the results using the `plot()` and `summary()`
]

---

# Example: Rocket launches

.large[
`$Y_i =$` 0 or 1 (failure/success)

.bold[Likelihood:] `$Y_i \overset{\rm iid}{\sim}{\rm Bernoulli} (\theta)$`

.bold[Prior:] Elicitation leads to uniform on (0.1, 0.9)

.bold[Posterior:]

`$$p(\theta | y) \propto \begin{cases}
\theta^3 (1-\theta)^8 & \text{if } 0.1 < \theta < 0.9\\
0 & \text{otherwise.}\end{cases}$$`
]

---

# 0. Load rjags and data

```r
library(rjags)

# Load data
launches <- read.table("https://www4.stat.ncsu.edu/~wilson/BR/table21.txt", header = TRUE)
y <- launches$Outcome
n <- nrow(launches)
```

---

# 1. Specify the model as a string

```r
model_string <- textConnection("model{

# Likelihood
  for(i in 1:n) {
    y[i] ~ dbern(theta)
  }
  
  # Prior
  theta ~ dunif(a, b)
  
}")
```

.content-box-yellow[
.large[
Be sure to use `textConnection()` if you are not saving this model string as a separate text file.
]
]

---
# 2. Compile the MCMC code

```r
inits <- list(theta = 0.5)
data  <- list(y = y, n = n, a = 0.1, b = 0.9)
model <- jags.model(model_string, data = data, 
                    inits = inits, n.chains = 1)
```

```
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 11
##    Unobserved stochastic nodes: 1
##    Total graph size: 15
## 
## Initializing model
```

.content-box-yellow[
.large[
`inits` and `data` should be **lists**
]
]

---

# 3. Draw burn-in samples

```r
update(model, 1000)
```

.content-box-yellow[
.large[
Add `progress.bar = "none"` to silence the progress bar in your .Rmd files]
]

---

# 4. Draw posterior samples

```r
samples <- coda.samples(
  model, 
  variable.names = "theta", 
  n.iter = 10000, 
  progress.bar = "none"
)
```

.content-box-yellow[
.large[
`variable.names` expects a character vector with the names of the variables to be monitored]
]

---

# 5. Inspect/summarize the results

```r
summary(samples)
```

```
## 
## Iterations = 2001:12000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##       0.314778       0.118194       0.001182       0.001574 
## 
## 2. Quantiles for each variable:
## 
##   2.5%    25%    50%    75%  97.5% 
## 0.1240 0.2245 0.3040 0.3927 0.5713
```

---

# 5. Inspect/summarize the results

```r
plot(samples)
```

![](18-jags_files/figure-html/unnamed-chunk-7-1.svg)

---

class: inverse

# Your turn

.Large[.bolder[Work through the examples on the JAGS handout with a neighbor.]]

.Large[.bolder[Check GitHub for a .Rmd template]]

.Large[.bolder[https://github.com/aloy/math315-fall2019]]