/
parameter_dist.R
59 lines (48 loc) · 1.6 KB
/
parameter_dist.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
load("35563325713.Rdat")
require(socialR)
require(warningsignals)
plot_m <- function(pow, X=NULL){
test <- pow$test_par_dist
shade_plot <- function(par, lower=-Inf, upper=Inf){
# disgard nonsense values
par <- par[par>lower & par < upper]
dist <- density(par)
plot(dist, type="n")
polygon(dist$x, dist$y, col="gray", border="darkgray")
}
upper <- 1
lower <- -1
if(!is.null(X)){
lower <- -10*length(X)/max(time(X))
upper <- -lower
}
shade_plot(test[2,], lower, upper)
}
par(mfrow=c(1,2))
plot(density(deterior_mc$test_par_dist[2,], bw=1), xlim=c(-20,1))
plot(density(constant_mc$test_par_dist[2,c(-173,-1717)]))
## parameter distributions
par_names <- c("Ro", "m", "theta", "sigma")
plot_parameters <- function(pow){
test <- pow$test_par_dist
null <- pow$null_par_dist
shade_plot <- function(par, lower=-Inf, upper=Inf, ...){
par <- par[par>lower & par < upper]
dist <- density(par)
plot(dist, type="n", ...)
polygon(dist$x, dist$y, col="gray", border="darkgray")
}
n <- dim(test)[1]
par(mfrow=c(2,n))
for(i in 1:n ){
shade_plot(test[i,], main=par_names[i])
}
m <- dim(null)[1]
for(i in 1:m ){
if(i==2) plot(test[2,], type="n")
shade_plot(null[i,], main="")
}
}
social_plot(plot_parameters(deterior_mc), tag="sim deterior warningsignals stochpop parameters")
social_plot(plot_parameters(constant_mc), tag="sim constant warningsignals stochpop parameters")
social_plot(plot_parameters(deut3_mc), tag="sim deut3 warningsignals stochpop parameters")