-
Notifications
You must be signed in to change notification settings - Fork 0
/
usage.Rmd
44 lines (33 loc) · 1.39 KB
/
usage.Rmd
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
---
title: "Gibbs Sampler Example Usages"
output: html_notebook
---
```{r}
library(latex2exp)
library(tidyverse)
source("samplers.R")
```
## Using the Multivariate Normal Gibbs Sampler
```{r}
y_blue <- as.matrix(read.table("data/bluecrab.dat"))
y_orange <- as.matrix(read.table("data/orangecrab.dat"))
post_blue <- gibbs_mvn(y_blue)
post_orange <- gibbs_mvn(y_orange)
blue_theta_plot <- ggplot(data.frame(post_blue$theta), aes(x = X1, y = X2)) +
geom_point(color = "blue", shape = ".") +
labs(x = TeX('$\\theta_1$'), y = TeX('$\\theta_2$')) +
ggtitle("Posterior Means for Blue Crab Data")
orange_theta_plot <- ggplot(data.frame(post_orange$theta), aes(x = X1, y = X2)) +
geom_point(color = "orange", shape = ".") +
labs(x = TeX('$\\theta_1$'), y = TeX('$\\theta_2$')) +
ggtitle("Posterior Means for Orange Crab Data")
grid.arrange(blue_theta_plot, orange_theta_plot, ncol = 2, nrow = 1)
```
## Using the Hierarchical Normal Gibbs Sampler
```{r}
filenames <- sapply(1:8, function(x) paste0("data/school", x, ".dat"))
school_data <- sapply(filenames, scan, quiet = TRUE)
post <- gibbs_hierarchical(school_data, mu_0 = mu_0, gamma_0_sq = 5,
tau_0_sq = tau_0_sq, eta_0 = eta_0,
sigma_0_sq = sig_0_sq, nu_0 = nu_0, num_iter = 1000)
```