-
Notifications
You must be signed in to change notification settings - Fork 0
/
grid_dispatch_B_sdvS.R
119 lines (90 loc) · 3.34 KB
/
grid_dispatch_B_sdvS.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
107
108
109
110
111
112
113
114
115
116
117
118
119
rm(list=ls())
source("dmc/dmc.R")
load_model("LBA", "lbaN_B.R")
# Load
print(load("samples/sTPPM_B_sdvS.RData"))
samples <- h.run.dmc(samples,
cores = 24,
report = 10,
p.migrate = 0.05)
save(samples, file = "samples/sTPPM_B_sdvS.RData")
# # Check for convergence
# gelman.diag.dmc(samples)
# check.recovery.dmc(samples)
samples1 <- h.RUN.dmc(h.samples.dmc(samples = samples,
nmc = 100,
thin = 10,
n.chains = 78),
cores = 24,
report = 10,
verbose = TRUE)
# Save samples
save(samples, samples1, file = "samples/sTPPM_B_sdvS.RData")
# samples2 <- h.RUN.dmc(h.samples.dmc(samples = samples1,
# nmc = 100,
# thin = 10,
# n.chains = 78),
# cores = 16,
# report = 10)
# # Save samples
# save(samples, samples1, samples2, file = "samples/sTPPM_B_sdvS.RData")
# Sampling plots ----------------------------------------------------------
# Load samples
load("samples/sTPPM_B_sdvS.RData")
# samples1 <- samples2
# Chains
pdf("plots/chains_B_sdvS.pdf", height = 6, width = 8)
par(mfrow = c(3, 3))
for(i in 1:length(samples1)) {
plot.dmc(samples1, pll.chain = TRUE, subject = i, start = 1)
title(paste0('Subject ', names(samples1)[i]))
}
dev.off()
# Chains by subject
pdf("plots/chains_by_sub_B_sdvS.pdf", height = 6, width = 8)
for(sub in names(samples1)) {
par(mfrow = c(1, 1))
plot.dmc(samples1, hyper = FALSE, subject = sub, pll.chain = TRUE)
mtext(paste0('Participant ', sub, ' posterior log likelihoods of all chains'),
outer = TRUE, line = -1.5)
plot.dmc(samples1, hyper = FALSE, subject = sub, layout = c(2, 2), density = TRUE)
mtext(paste0('Participant ', sub), outer = TRUE, line = -1.5)
}
dev.off()
# Posterior predictive fits -----------------------------------------------
# Sample posterior predictives
post_preds <- h.post.predict.dmc(samples1)
save(post_preds, file = "deriv/post_preds_B_sdvS.RData")
# Load posterior predictives
load("deriv/post_preds_B_sdvS.RData")
# Plot fits to CDFs
pdf("plots/cdfs_B_sdvS.pdf", height = 6, width = 8)
plot.pp.dmc(post_preds, "cdf", layout = c(2, 2), model.legend = FALSE)
dev.off()
# Subject-level fits
pdf("plots/cdfs_by_sub_B_sdvS.pdf", height = 6, width = 8)
for (i in seq_len(length(post_preds))) {
plot.pp.dmc(post_preds[[i]], "cdf", layout = c(2, 2), model.legend = FALSE)
}
dev.off()
# Fit indices
h.IC.dmc(samples1, DIC = TRUE)
# Save posterior predictive sims
post_pred_sims <- h.post.predict.dmc(samples1, save.simulation = TRUE)
save(post_pred_sims, file = "deriv/post_pred_sims_B_sdvS.RData")
# Parameters --------------------------------------------------------------
# Summary
parms <- summary.dmc(samples1)
parms <- do.call(rbind, lapply(parms, function(x) x$statistics[, 1]))
# head(parms)
save(parms, file = "deriv/map_parms_B_sdvS.RData")
# Recovery
h.check.recovery.dmc(samples1, digits = 2)
# Effective sample size
ess <- effectiveSize.dmc(samples1)
do.call(rbind, ess)
min(do.call(rbind, ess))
# # Correlations
# pdf("plots/pairs_B_sdvS.pdf", height = 6, width = 8)
# pairs.dmc(samples1, start = 50)
# dev.off()