bayesian-learning-with-R/rjags_intro.Rmd at master · NCSU-STATS/bayesian-learning-with-R · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
---
title: "Introduction to rjags "
subtitle: "Just Another Gibbs Sampler"
author: "Bruce Campbell"
fontsize: 11pt
output: pdf_document
bibliography: BruceCampbell_ST540_HW_1.bib
---
---

```{r setup, include=FALSE,echo=FALSE}
rm(list = ls())
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(dev = 'pdf')
knitr::opts_chunk$set(cache=TRUE)
knitr::opts_chunk$set(tidy=TRUE)
knitr::opts_chunk$set(prompt=FALSE)
knitr::opts_chunk$set(fig.height=5)
knitr::opts_chunk$set(fig.width=7)
knitr::opts_chunk$set(warning=FALSE)
knitr::opts_chunk$set(message=FALSE)
knitr::opts_knit$set(root.dir = ".")
knitr::opts_chunk$set(tidy.opts=list(width.cutoff=38),tidy=TRUE)
library(latex2exp)
library(pander)
library(ggplot2)
library(GGally)

```

##

```{r}
setwd("C:/E/brucebcampbell-git/bayesian-learning-with-R")
library(rjags)
N <- 2000
x <- rnorm(N, 0, 5)

write.table(x,
            file = 'example1.data',
            row.names = FALSE,
            col.names = FALSE)

jags <- jags.model('jags_normal_model.bug',
                   data = list('x' = x,
                               'N' = N),
                   n.chains = 10,
                   n.adapt = 100)

update(jags, 1000)

jags.samples(jags,
             c('mu', 'tau'),
             4000)

```

##Coda can be used to provide debugging information for JAGS
```{r}
library(coda)

samples <- coda.samples(jags,
                        c('mu', 'tau'),
                        1000)
plot(samples)
```


### Beta Binomial Model in JAGS

```{r}
n      <- 20
Y      <- 4
a      <- 3
b      <- 1
model_string <- "model{

  # Likelihood
  Y ~ dbinom(theta,n)

  # Prior
  theta ~ dbeta(a, b)
}"
model <- jags.model(textConnection(model_string),
                    data = list(Y=Y,n=n,a=a,b=b))
update(model, 10000, progress.bar="none"); # Burnin for 10000 samples

samp <- coda.samples(model,
        variable.names=c("theta"),
        n.iter=20000, progress.bar="none")

summary(samp)
plot(samp)
```