diff --git a/IntervalsModel/plots/Distributions.Poisson_hh.pdf b/IntervalsModel/plots/Distributions.Poisson_hh.pdf
index 76f28596b32d2dff18b708363c1bbd952704ec3d..b9fce03e61511d16e5f2fbb75458c6718726dd7c 100644
Binary files a/IntervalsModel/plots/Distributions.Poisson_hh.pdf and b/IntervalsModel/plots/Distributions.Poisson_hh.pdf differ
diff --git a/IntervalsModel/plots/Distributions.Poisson_rest.pdf b/IntervalsModel/plots/Distributions.Poisson_rest.pdf
index 04bcddc7c3dcc3c07b4880734b7dfed3d27c013c..a5f91163fd581ff87018b85000917701c1b76e34 100644
Binary files a/IntervalsModel/plots/Distributions.Poisson_rest.pdf and b/IntervalsModel/plots/Distributions.Poisson_rest.pdf differ
diff --git a/IntervalsModel/plots/Distributions.Poisson_ws.pdf b/IntervalsModel/plots/Distributions.Poisson_ws.pdf
index 8e09fe4b6c6c16b1182c5ae8cc4db8a63b00831c..a6cd3b525241c59754d86e746ec6fa51b95f2761 100644
Binary files a/IntervalsModel/plots/Distributions.Poisson_ws.pdf and b/IntervalsModel/plots/Distributions.Poisson_ws.pdf differ
diff --git a/IntervalsModel/plots/hh.pdf b/IntervalsModel/plots/hh.pdf
index fa2572fa0a1315ff43a946aadaa07b71c4ca7b8a..a10a3fe0fcbcc978e9aac893526077aa9822f23b 100644
Binary files a/IntervalsModel/plots/hh.pdf and b/IntervalsModel/plots/hh.pdf differ
diff --git a/IntervalsModel/plots/rest.pdf b/IntervalsModel/plots/rest.pdf
index 9f8728a9d14cf1cc15f8dce54ad848a05f4ed09a..575d1b48b5ed83553a2023d0d41bc6539d16b387 100644
Binary files a/IntervalsModel/plots/rest.pdf and b/IntervalsModel/plots/rest.pdf differ
diff --git a/IntervalsModel/plots/ws.pdf b/IntervalsModel/plots/ws.pdf
index 1e36f2f154e101e73e18e7d0d4b744ed260b4bce..db7380cb4c5e13b476f055d40db3f05c8e050c02 100644
Binary files a/IntervalsModel/plots/ws.pdf and b/IntervalsModel/plots/ws.pdf differ
diff --git a/IntervalsModel/simulation_data/hh.dat b/IntervalsModel/simulation_data/hh.dat
new file mode 100644
index 0000000000000000000000000000000000000000..1ef311949c975c6f24703ef4e2dea20b79699043
Binary files /dev/null and b/IntervalsModel/simulation_data/hh.dat differ
diff --git a/IntervalsModel/simulation_data/ws.dat b/IntervalsModel/simulation_data/ws.dat
index a20eeb8fe50afeeac37ba5c0332ec3777938896b..9017a988d8706c87779aa79bed4776b30a3b53af 100644
Binary files a/IntervalsModel/simulation_data/ws.dat and b/IntervalsModel/simulation_data/ws.dat differ
diff --git a/IntervalsModel/src/IntervalsModel.jl b/IntervalsModel/src/IntervalsModel.jl
index 51275b5cda728d45e2a1abc755c7cd4231f06eb8..7d0ec6b8a3ec8b47313c3482e24405e3fb3ffd96 100644
--- a/IntervalsModel/src/IntervalsModel.jl
+++ b/IntervalsModel/src/IntervalsModel.jl
@@ -25,7 +25,7 @@ const rng = Xoroshiro128Plus()
#consts that let give us nicer names for the indices
const YOUNG, MIDDLE,OLD = 1,2,3
const durmax = 144
-const sample_repeat = 10
+const sample_repeat = 100
"""
Number of start times to sample for a given set of distribution parameters.
"""
@@ -44,9 +44,9 @@ include("plotting_functions.jl")
Runs parameter estimation for the three scenarios, hh, ws, and rest.
"""
function main()
- # do_hh(400)
- do_ws(400)
- do_rest(400)
+ do_hh(500)
+ do_ws(500)
+ do_rest(500)
plot_all()
end
@@ -88,7 +88,7 @@ Number of particles to obtain representing posterior distribution of the distrib
Threshold adaptive rate, see `KissABC.smc` for more details.
"""
-function bayesian_estimation(fname, err_func, priors_list, dists, particles; alpha = 0.99)
+function bayesian_estimation(fname, err_func, priors_list, dists, particles; alpha = 0.98)
data_pairs = map(zip(dists,priors_list)) do (dist,priors)
init = rand(priors)
@@ -112,11 +112,12 @@ function do_hh(particles)
]
poisson_priors = filter_priors("home")
+ display(poisson_priors)
# Set parameter priors for fitting
priors_list = [
Factored(poisson_priors...),
]
- bayesian_estimation("hh",err_hh,priors_list,dists,particles; alpha = 0.99)
+ bayesian_estimation("hh",err_hh,priors_list,dists,particles; alpha = 0.98)
end
@@ -136,7 +137,7 @@ function do_ws(particles)
priors_list = [
Factored(poisson_priors...)
]
- bayesian_estimation("ws",err_ws,priors_list,dists,particles; alpha = 0.99)
+ bayesian_estimation("ws",err_ws,priors_list,dists,particles; alpha = 0.98)
end
"""
diff --git a/IntervalsModel/src/hh_durations_model.jl b/IntervalsModel/src/hh_durations_model.jl
index 18f65bffb5aaa4a2a91eb0ded9ef691c6038c37a..79a9d17814bf2a4067750aa37fe118125c43809d 100644
--- a/IntervalsModel/src/hh_durations_model.jl
+++ b/IntervalsModel/src/hh_durations_model.jl
@@ -38,7 +38,7 @@ function err_hh(p,dist)
# display(p_matrix)
AGERESP = dat.AGERESP #age of the respondent
errsum = 0
- @inbounds for i = 1:length(duration_subarray) #loop through entire file
+ @inbounds for i = 1:length(AGERESP) #loop through entire file
age_sample = AGERESP[i]
@inbounds for age_j in YOUNG:OLD #for a given age_sample loop over possible contact ages
diff --git a/IntervalsModel/src/plotting_functions.jl b/IntervalsModel/src/plotting_functions.jl
index 1e85e408c16044a3ed65de8f2a73e10c81c8861f..9691bd085a3af45c90033ff0d59090f3be7adc26 100644
--- a/IntervalsModel/src/plotting_functions.jl
+++ b/IntervalsModel/src/plotting_functions.jl
@@ -1,6 +1,7 @@
using LaTeXStrings
using Plots
-const color_palette = palette(:seaborn_bright) #color theme for the plots
+const color_palette = palette(:seaborn_pastel) #color theme for the plots
+const color_palette_bright = palette(:seaborn_bright) #bright color theme for the plots
default(dpi = 300)
default(framestyle = :box)
@@ -34,7 +35,12 @@ function plot_dists(fname,dist_constructors,data)
hasnans = any(any.(map(l -> isnan.(l),dist_pts)))
err_down = hasnans ? 0 : quantile.(dist_pts,0.05)
err_up = hasnans ? 0 : quantile.(dist_pts,0.95)
- plot!(p_matrix[i,j] ,x_range,mean_dat; ribbon = ( mean_dat .- err_down,err_up .- mean_dat),legend = false,label = string(dist_constructor),seriescolor = color_palette[k])
+ plot!(p_matrix[i,j] ,x_range,mean_dat;
+ ribbon = ( mean_dat .- err_down,err_up .- mean_dat),
+ legend = false,
+ label = string(dist_constructor),
+ seriescolor = color_palette_bright[k]
+ )
end
annotate!(p_matrix[i,j],compute_x_pos(p_matrix[i,j]),compute_y_pos(p_matrix[i,j]), Plots.text("$(ymo[i])→$(ymo[j])", :left, 10))
end
@@ -54,7 +60,7 @@ function plot_posteriors(fname,data)
hist = fit(Histogram,data.P[i].particles; nbins = 40)
kde_est = kde(data.P[i].particles)
kernel_data = [pdf(kde_est,x) for x in minimum(data.P[i].particles):maximum(data.P[i].particles)]
- plot!(p_list[i],hist;legend = false, xlabel = L"\lambda_%$i")
+ plot!(p_list[i],hist;legend = false, xlabel = L"\lambda_%$i", color = color_palette[1] )
# display(kernel_data)
# vline!(p_list[i],[argmax(kernel_data)]; seriescolor = color_palette[2])