-
Notifications
You must be signed in to change notification settings - Fork 0
/
template.R
130 lines (119 loc) · 5.8 KB
/
template.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
120
121
122
123
124
125
126
127
128
129
td <- 12
inits <- pilot.inits[[m]][1:reps]
recompile <- FALSE # recomile .tpl files?
set.seed(seed)
seeds <- sample(1:1e4, size=reps)
d <- m
## Run model with hbf=1 to get right covariance matrix and MLEs for NUTS
setwd(d)
if(recompile)
system(paste('admb', m, '-f')) # -f does optimized mode
## Add the -mcmc so SS turns off bias adjustment and we get a covariance
## that matches this to use for sampling. Otherwise there is a mismatch.
if(m!='snowcrab2'){
system(paste(m, '-hbf 1 -mcmc 15 -iprint 100 -nox'))
} else {
system(paste(m, '-hbf 1 -iprint 100 -nox -ainp snowcrab2.par -phase 50'))
}
setwd('..')
### Run NUTS for different mass matrices
## Mass matrix is MLE covariance
fit.nuts.mle <-
sample_admb(m, iter=iter, init=inits, algorithm='NUTS', seeds=seeds,
parallel=TRUE, chains=reps, warmup=warmup, path=d, cores=reps,
control=list(max_treedepth=td, metric="mle", adapt_delta=ad))
fit.nuts.mle <- add.monitor(fit.nuts.mle, FALSE)
## Mass matrix is dense one estimated from previous run.
fit.nuts.dense <-
sample_admb(m, iter=iter, init=inits, algorithm='NUTS', seeds=seeds,
mceval=TRUE,
parallel=TRUE, chains=reps, warmup=warmup, path=d, cores=reps,
control=list(max_treedepth=td, metric=fit.nuts.mle$covar.est, adapt_delta=ad))
## These are presumably the best posterior samples so use them for
## management quantities by running mceval and then saving the results onto
## the output for later use in plotting and analysis.
fit.nuts.dense <- add.monitor(fit.nuts.dense, FALSE, metric='dense')
fit.nuts.dense$dq.post <- get.dq(m, TRUE)
## Now run RWM but using a thinning rate similar to NUTS so the time is
## roughly equivalent.
## Rerun model with hbf=0
setwd(m)
if(m!='snowcrab2'){
system(paste(m, '-hbf 0 -mcmc 15 -nox -iprint 100'))
} else {
system(paste0(m, ' -nox -phase 50 -iprint 100 -ainp ', m,'.par'))
# system(paste(m, '-hbf 0 -ainp snowcrab2.par -phase 50 -nox'))
}
setwd('..')
tt <- floor(4*mean(extract_sampler_params( fit.nuts.mle)$n_leapfrog__))
fit.rwm.mle <-
sample_admb(m, iter=tt*iter, init=inits, thin=tt, seeds=seeds,
parallel=TRUE, chains=reps, warmup=tt*warmup,
path=d, cores=reps, control=list(metric=NULL),
algorithm='RWM')
fit.rwm.mle <- add.monitor(fit.rwm.mle, FALSE)
### This just didn't really work so took it out
## tt <- floor(4*mean(extract_sampler_params( fit.nuts.dense)$n_leapfrog__))
## fit.rwm.dense <-
## sample_admb(m, iter=tt*iter, init=inits, thin=tt, seeds=seeds,
## parallel=TRUE, chains=reps, warmup=tt*warmup,
## path=d, cores=reps, control=list(metric=fit.rwm.mle$covar.est),
## algorithm='RWM')
fit.rwm.dense <- NULL
## Save fits to file
message(paste("Saving results to file for model",m))
saveRDS(list(fit.nuts.mle=fit.nuts.mle, fit.nuts.dense=fit.nuts.dense,
fit.rwm.mle=fit.rwm.mle, fit.rwm.dense=fit.rwm.dense),
file=paste0('results/', d,'_fits.RDS'))
message(paste("Calculating metrics for model",m))
### Gather adaptation and performance metrics
ff <- function(labels, ...){
fits <- list(...)
if(length(fits)!=length(labels)) stop("bad labels")
do.call(rbind, lapply(1:length(fits), function(i){
x <- fits[[i]]$sampler_params
data.frame(m=labels[i], chain=1:length(x),
eps=as.numeric(do.call(rbind, lapply(x, function(l) tail(l[,2],1)))),
divergences=as.numeric(do.call(rbind, lapply(x, function(l) sum(l[-(1:fits[[i]]$warmup),5])))),
accept_prob=as.numeric(do.call(rbind, lapply(x,
function(l) mean(l[-(1:fits[[i]]$warmup),1])))),
nsteps=as.numeric(do.call(rbind, lapply(x,
function(l) mean(l[-(1:fits[[i]]$warmup),4])))))
}))}
adaptation <- ff(c("mle", "dense"), fit.nuts.mle, fit.nuts.dense)
stats.nuts.mle <- with(fit.nuts.mle,
data.frame(alg='nuts', m='mle', time.total=sum(time.total), monitor))
perf.nuts.mle <- data.frame(alg='nuts', m='mle',
efficiency=min(stats.nuts.mle$n_eff)/sum(fit.nuts.mle$time.total))
stats.nuts.dense <- with(fit.nuts.dense,
data.frame(alg='nuts', m='dense', time.total=sum(time.total), monitor))
perf.nuts.dense <- data.frame(alg='nuts', m='dense',
efficiency=min(stats.nuts.dense$n_eff)/sum(fit.nuts.dense$time.total))
stats.rwm.mle <- with(fit.rwm.mle,
data.frame(alg='rwm', m='mle', time.total=sum(time.total), monitor))
perf.rwm.mle <- data.frame(alg='rwm', m='mle',
efficiency=min(stats.rwm.mle$n_eff)/sum(fit.rwm.mle$time.total))
## stats.rwm.dense <- with(fit.rwm.dense, data.frame(alg='rwm', m='dense', time.total=sum(time.total), rstan::monitor(samples, warmup=warmup, probs=.5, print=FALSE)))
## perf.rwm.dense <- data.frame(alg='rwm', m='dense', efficiency=min(stats.rwm.dense$n_eff)/sum(fit.rwm.dense$time.total))
stats.all <- rbind(stats.nuts.mle, stats.nuts.dense, stats.rwm.mle)
stats.all[,c('mean', 'se_mean', 'sd', 'X50.', 'model', 'metric', 'version',
'alg.1', 'pilot')] <- NULL
stats.all <- ddply(stats.all, .(alg, m), mutate, perf=(n_eff)/time.total)
stats.long <- reshape2::melt(stats.all, c('alg', 'm'))
perf.all <- rbind(perf.nuts.mle, perf.nuts.dense, perf.rwm.mle)
adaptation.long <- reshape2::melt(adaptation, c('m', 'chain'))
## Quick plots
message(paste("Making plots for model",m))
ggwidth <- 7
ggheight <- 5
pdf(paste0('plots/', d, '_comparison.pdf'), width=7, height=5)
g <- ggplot(adaptation.long, aes(y=value, x=m)) + geom_point(alpha=.5) +
facet_wrap('variable', scales='free')
print(g)
g <- ggplot(stats.long, aes(y=value, x=m, color=alg)) + geom_jitter(alpha=.5) +
facet_wrap('variable', scales='free')
print(g)
g <- ggplot(perf.all, aes(m, efficiency, color=alg)) + geom_point()
print(g)
plot.ess(rwm=fit.rwm.mle, nuts=fit.nuts.mle)
dev.off()