/
simple_example.R
107 lines (97 loc) · 3.83 KB
/
simple_example.R
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
95
96
97
98
99
100
101
102
103
104
105
106
library(ggplot2)
library(dplyr)
library(tidyr)
theme_set(bayesplot::theme_default(base_family = "sans"))
source("code/monitornew.R")
set.seed(666)
nreps = 1000
nchains=4;
niter = 1000;
npars=4
sd = matrix(data = 1.0,nrow=nchains, ncol= npars-2);
sd[1,1] = 1/3;
means = c(2,0,0,0)
ar_param = 0.3
Rhat = data.frame(rhat = rep(NA,2*npars*nreps),
version = rep(NA,2*npars*nreps),
par = rep(NA,2*npars*nreps),
neff = rep(NA,2*npars*nreps),
neff_bulk = rep(NA,2*npars*nreps),
neff_tail = rep(NA,2*npars*nreps))
sims = array(data=0.0,dim = c(niter,nchains,npars))
for (rep in 1:nreps) {
for(chain in 1:nchains) {
for (par in 1:2) {
sims[,chain,par] = arima.sim(n=niter, list(ar=ar_param),sd=sqrt(1-ar_param^2))*sd[chain,par]
}
sims[,chain,3] = (arima.sim(n=niter, list(ar=ar_param),sd=sqrt(1-ar_param^2)))/arima.sim(n=niter, list(ar=ar_param),sd=sqrt(1-ar_param^2)) + means[chain]
sims[,chain,4] = arima.sim(n=niter, list(ar=ar_param),sd=sqrt(1-ar_param^2))/arima.sim(n=niter, list(ar=ar_param),sd=sqrt(1-ar_param^2))
}
index_start = 2*npars*(rep-1) + 1;
index_end = 2*npars*rep;
nothing <- capture.output(monitor_old <- rstan::monitor(sims,warmup=0))
nothing <- capture.output(monitor_new <- monitor(sims,warmup=0))
par_names <- c("Finite mean, different variance", "Finite mean, same variance", "Infinite mean, different location","Infinite mean, same location")
Rhat$rhat[index_start:index_end] = c(monitor_old[,10],monitor_new$Rhat)
Rhat$version[index_start:index_end] = c(rep("old",npars),rep("new",npars))
Rhat$par[index_start:index_end] = c(par_names,par_names)
Rhat$neff[index_start:index_end] = c(monitor_old[,9],monitor_new$Bulk_ESS)
Rhat$neff_bulk[index_start:index_end] = c(rep(NA,npars),monitor_new$Bulk_ESS)
Rhat$neff_tail[index_start:index_end] = c(rep(NA,npars), monitor_new$Tail_ESS)
}
Rhat %>% ggplot(aes(rhat,colour=NULL, fill=version)) +
geom_histogram(position="identity", alpha=0.9) +
facet_wrap(~par,nrow = 2,ncol = 2, scales = "free") +
scale_fill_manual(
name=paste("Version of Rhat"),
breaks=c("new", "old"),
labels=c("This paper","Gelman et al. (2013)"),
values=c("new" = "#d1e1ec", "old" ="#a25050")
) +
guides(colour=FALSE) +
labs(y="", x="Rhat") +
xlim(c(0.99, 1.15)) +
theme(
text = element_text(size=16),
strip.text = element_text(size=16),
axis.title = element_text(size=16),
axis.text = element_text(size=16),
axis.line = element_line(),
axis.line.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
)
ggsave(file="code/simple_rhat_compare.png",width=12, height = 7.5, units="in")
ggsave(file="paper/graphics/simple_rhat_compare.png",width=12, height = 7.5, units="in")
Rhat %>% ggplot(aes(neff,colour=version, fill=version)) +
geom_histogram() +
facet_wrap(~par,nrow = 2,ncol = 2) +
labs(y = "") +
theme_minimal() +
theme(
strip.text = element_text(size=16),
axis.title = element_text(size=16),
axis.text = element_text(size=16),
axis.line.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
)
ggsave(file="code/simple_neff_compare.png",width=12, height = 7.5, units="in")
Rhat %>%
select(neff_bulk,neff_tail,par) %>%
gather(`neff_bulk`,`neff_tail`,key="type",value="neff") %>%
drop_na() %>%
ggplot(aes(x=neff,colour=type,fill=type)) +
geom_histogram() +
facet_wrap(~par,nrow = 2,ncol = 2) +
labs(y = "") +
theme_minimal() +
theme(
strip.text = element_text(size=16),
axis.title = element_text(size=16),
axis.text = element_text(size=16),
axis.line.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
)
ggsave(file="code/simple_bulk_vs_tail.png",width=12, height = 7.5, units="in")