diff --git a/IntervalsModel/Manifest.toml b/IntervalsModel/Manifest.toml
index b631eb5ad75dd8b834f83f549385f62ae064f0b9..aae98665fa39cfd0a2853329a0ec1cc328d3efbd 100644
--- a/IntervalsModel/Manifest.toml
+++ b/IntervalsModel/Manifest.toml
@@ -184,6 +184,12 @@ version = "1.0.0"
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+[[DefaultApplication]]
+deps = ["InteractiveUtils"]
+git-tree-sha1 = "fc2b7122761b22c87fec8bf2ea4dc4563d9f8c24"
+uuid = "3f0dd361-4fe0-5fc6-8523-80b14ec94d85"
+version = "1.0.0"
+
[[DefineSingletons]]
git-tree-sha1 = "77b4ca280084423b728662fe040e5ff8819347c5"
uuid = "244e2a9f-e319-4986-a169-4d1fe445cd52"
@@ -203,6 +209,12 @@ git-tree-sha1 = "f0e06a5b5ccda38e2fb8f59d91316e657b67047d"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.24.12"
+[[DocStringExtensions]]
+deps = ["LibGit2", "Markdown", "Pkg", "Test"]
+git-tree-sha1 = "50ddf44c53698f5e784bbebb3f4b21c5807401b1"
+uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+version = "0.8.3"
+
[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
@@ -684,6 +696,18 @@ git-tree-sha1 = "95a4038d1011dfdbde7cecd2ad0ac411e53ab1bc"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.10.1"
+[[PGFPlotsX]]
+deps = ["ArgCheck", "DataStructures", "Dates", "DefaultApplication", "DocStringExtensions", "MacroTools", "Parameters", "Requires", "Tables"]
+git-tree-sha1 = "1adde3d07cce96b6a3bb88572612db4bd9d6153b"
+uuid = "8314cec4-20b6-5062-9cdb-752b83310925"
+version = "1.2.10"
+
+[[Parameters]]
+deps = ["OrderedCollections", "UnPack"]
+git-tree-sha1 = "2276ac65f1e236e0a6ea70baff3f62ad4c625345"
+uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+version = "0.12.2"
+
[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "50c9a9ed8c714945e01cd53a21007ed3865ed714"
diff --git a/IntervalsModel/Project.toml b/IntervalsModel/Project.toml
index 58ade50c253cbfd3ccc4dc94a0ce5c11a6aca850..da17ab876bfd66563f04a8ebba678cc559be09d9 100644
--- a/IntervalsModel/Project.toml
+++ b/IntervalsModel/Project.toml
@@ -12,7 +12,9 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
JLSO = "9da8a3cd-07a3-59c0-a743-3fdc52c30d11"
KissABC = "9c9dad79-530a-4643-a18b-2704674d4108"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
ModularIntervals = "1e07f51d-1de3-4a1c-a782-ac336877d585"
+PGFPlotsX = "8314cec4-20b6-5062-9cdb-752b83310925"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
RandomNumbers = "e6cf234a-135c-5ec9-84dd-332b85af5143"
Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
diff --git a/IntervalsModel/plots/hh.pdf b/IntervalsModel/plots/hh.pdf
index ff07d86f1768d78ebcb85ba23d876fd21400036e..e4475483c8bb20c1286e6008921c7c46baff0758 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
new file mode 100644
index 0000000000000000000000000000000000000000..1b5c6eb890568ffb59e5963a4c9d1ab6f4c8bd42
Binary files /dev/null and b/IntervalsModel/plots/rest.pdf differ
diff --git a/IntervalsModel/plots/ws.pdf b/IntervalsModel/plots/ws.pdf
index ae202cdaa2684da2bc800d06d34a630bd163ee7b..d46cb378594dd87fb76401b8ec447b228576f970 100644
Binary files a/IntervalsModel/plots/ws.pdf and b/IntervalsModel/plots/ws.pdf differ
diff --git a/IntervalsModel/simulation_data/rest.dat b/IntervalsModel/simulation_data/rest.dat
new file mode 100644
index 0000000000000000000000000000000000000000..ed9fae52b48b1c7c870aa9a0062a38ccddc8479e
Binary files /dev/null and b/IntervalsModel/simulation_data/rest.dat differ
diff --git a/IntervalsModel/simulation_data/ws.dat b/IntervalsModel/simulation_data/ws.dat
index 7c51dccc9f356d9700754b1d2ec1da35dabfb87a..6bb21aa2aa083ff4534e8c7ec019f774982de332 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 929c5c12de8fd17d70b87fa056638cea5a828b06..c9e6b9a555e7dbbd06b4683e5d2f3c86f7dbc0fa 100644
--- a/IntervalsModel/src/IntervalsModel.jl
+++ b/IntervalsModel/src/IntervalsModel.jl
@@ -9,7 +9,7 @@ using Distributions
using CovidAlertVaccinationModel
using ZeroWeightedDistributions
const PACKAGE_FOLDER = dirname(dirname(pathof(IntervalsModel)))
-
+using DataStructures:OrderedDict
using Serialization
const rng = Xoroshiro128Plus(1)
@@ -31,7 +31,7 @@ using Plots
const μ_bounds = (6,12*6)
const σ_bounds = (1,48)
-const α_bounds = (0.0,1.0)
+const α_bounds = (0.0,0.1)
function do_hh(particles)
dists = [
@@ -63,6 +63,23 @@ function do_ws(particles)
@show bounds_list
bayesian_estimate("ws",err_ws,dists,bounds_list,particles)
+end
+function do_rest(particles)
+ dists = [
+ ZWDist{Normal},
+ ZWDist{Poisson},
+ ZWDist{Weibull}
+ ]
+
+ # Set parameter bounds for fitting
+ bounds_list = map(l -> vcat(l...),[
+ [[α_bounds for i = 1:6],[μ_bounds for i = 1:6], [σ_bounds for i = 1:6]],
+ [[α_bounds for i = 1:6],[μ_bounds for i = 1:6]],
+ [[α_bounds for i = 1:6],[μ_bounds for i = 1:6], [σ_bounds for i = 1:6]],
+ ])
+ @show bounds_list
+ bayesian_estimate("rest",err_rest,dists,bounds_list,particles)
+
end
function bayesian_estimate(fname,err_func,dists,bounds_list,particles)
@@ -70,22 +87,27 @@ function bayesian_estimate(fname,err_func,dists,bounds_list,particles)
init = [b[1] for b in bounds]
priors = Factored([Uniform(l,u) for (l,u) in bounds]...) #assume uniform priors
- # @btime err_ws($init,$dist) #compute benchmark of the error function, not rly necessary
- # @btime err_ws($init,$dist) #compute benchmark of the error function, not rly necessary
+ @btime err_ws($init,$dist) #compute benchmark of the error function, not rly necessary
- out = smc(priors,p -> err_func(p, dist), verbose=true, nparticles=particles, alpha=0.95, parallel = true) #apply sequential monte carlo with 200 particles
+ out = smc(priors,p -> err_func(p, dist), verbose=true, nparticles=particles, alpha=0.99, M = 1,epstol = 0.1, parallel = true) #apply sequential monte carlo with 200 particles
return dist => out
- end |> Dict
+ end |> OrderedDict
+ display(data_pairs)
serialize(joinpath(PACKAGE_FOLDER,"simulation_data","$fname.dat"),data_pairs)
end
function plot_estimates()
data = deserialize(joinpath(PACKAGE_FOLDER,"simulation_data","ws.dat"))
plot_dists("ws",collect(keys(data)),collect(values(data)))
+ display(collect(values(data)))
data = deserialize(joinpath(PACKAGE_FOLDER,"simulation_data","hh.dat"))
plot_dists("hh",collect(keys(data)),collect(values(data)))
+
+
+ data = deserialize(joinpath(PACKAGE_FOLDER,"simulation_data","rest.dat"))
+ plot_dists("rest",collect(keys(data)),collect(values(data)))
end
diff --git a/IntervalsModel/src/plotting_functions.jl b/IntervalsModel/src/plotting_functions.jl
index c60cba2e7b913d731aa34086bc41332c0b6772a6..6d25eb3070042389699a853c8f0f7d3b30f6f178 100644
--- a/IntervalsModel/src/plotting_functions.jl
+++ b/IntervalsModel/src/plotting_functions.jl
@@ -1,9 +1,15 @@
+using LaTeXStrings
+using Plots
+
+default(dpi = 300)
+default(framestyle = :box)
+pgfplotsx()
function plot_dists(fname,dist_constructors,data)
p_estimate_as_arrays = map(d -> get_params(d.P),data)
p_matrix = map(x -> plot(),p_estimate_as_arrays[1])
- for (dist_constructor,p_estimate) in zip(dist_constructors,p_estimate_as_arrays)
- for i in YOUNG:OLD, j in YOUNG:OLD
- @show p_estimate[i,j]
+ ymo = ["Y","M","O"]
+ for i in YOUNG:OLD, j in YOUNG:OLD
+ for (dist_constructor,p_estimate) in zip(dist_constructors,p_estimate_as_arrays)
dists = map(p -> p.particles, p_estimate[i,j]) |>
t -> zip(t...) |>
l -> map(t -> dist_constructor(t...),l)
@@ -13,20 +19,33 @@ function plot_dists(fname,dist_constructors,data)
err_up = quantile.(dist_pts,0.95)
plot!(p_matrix[i,j] ,0:144,mean_dat; ribbon = ( mean_dat .- err_down,err_up .- mean_dat),legend = false,label = string(dist_constructor))
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
plot!(p_matrix[end,1]; legend = true)
- p = plot(p_matrix..., size = (600,400))
+ p = plot(p_matrix..., size = (800,600))
savefig(p,joinpath(PACKAGE_FOLDER,"plots", "$fname.pdf"))
end
+compute_x_pos(p) = xlims(p)[1] + 0.02*((xlims(p)[2] - xlims(p)[1]))
+compute_y_pos(p) = ylims(p)[2] - 0.11*((ylims(p)[2] - ylims(p)[1]))
function plot_posteriors(fname,parameter_names,data)
for i in 1:length(data.P)
a = stephist(
data.P[i].particles;
normalize = true,
- title = i <=6 ? "μ_$i" : "σ_$i"
+ title = "$(parameter_names[div(i,6)])_(i % 6)"
)
push!(p_list,a)
end
p = plot(p_list...)
- savefig(p,"$norm.png")
-end
\ No newline at end of file
+ savefig(p,"$fname.png")
+end
+# function add_subplot_letters!(plot_list; pos = :top)
+# for (i,sp) in enumerate(plot_list)
+# letter = string(Char(i+96))
+# if pos == :top
+# annotate!(sp,xlims(sp)[1] + 0.02*((xlims(sp)[2] - xlims(sp)[1])),ylims(sp)[2] - 0.11*((ylims(sp)[2] - ylims(sp)[1])), Plots.text("$letter)", :left, 18))
+# elseif pos == :bottom
+# annotate!(sp,xlims(sp)[1] + 0.02*((xlims(sp)[2] - xlims(sp)[1])),ylims(sp)[1] + 0.11*((ylims(sp)[2] - ylims(sp)[1])), Plots.text("$letter)", :left, 18))
+# end
+# end
+# end
\ No newline at end of file
diff --git a/IntervalsModel/src/ws_durations_model.jl b/IntervalsModel/src/ws_durations_model.jl
index 735a8f2a7e2c109ad1ec4110ef41e4ed5810fa82..e27869e9c7b83827526c1619930f740da32458a9 100644
--- a/IntervalsModel/src/ws_durations_model.jl
+++ b/IntervalsModel/src/ws_durations_model.jl
@@ -5,10 +5,23 @@ const ws_data = (
M = CSV.File("$PACKAGE_FOLDER/network-data/Timeuse/WS/WorkschoolDataM.csv") |> Tables.matrix |> x -> dropdims(x;dims = 2),
O = CSV.File("$PACKAGE_FOLDER/network-data/Timeuse/WS/WorkschoolDataO.csv") |> Tables.matrix |> x -> dropdims(x;dims = 2),
)
+using KernelDensity
+const kerneldensity_data = (
+ Y = kde(IntervalsModel.ws_data.Y),
+ M = kde(IntervalsModel.ws_data.M),
+ O = kde(IntervalsModel.ws_data.O),
+)
+const rest_data = (
+ Y = CSV.File("$PACKAGE_FOLDER/network-data/Timeuse/Rest/RDataY.csv") |> Tables.matrix |> x -> dropdims(x;dims = 2),
+ M = CSV.File("$PACKAGE_FOLDER/network-data/Timeuse/Rest/RDataM.csv") |> Tables.matrix |> x -> dropdims(x;dims = 2),
+ O = CSV.File("$PACKAGE_FOLDER/network-data/Timeuse/Rest/RDataO.csv") |> Tables.matrix |> x -> dropdims(x;dims = 2),
+)
const comparison_samples = 10_000
ws_distributions = CovidAlertVaccinationModel.initial_workschool_mixing_matrix
+rest_distributions = CovidAlertVaccinationModel.initial_rest_mixing_matrix
+
#error function
function err_ws(p,dist)
params = get_params(p)
@@ -16,9 +29,37 @@ function err_ws(p,dist)
age_dists = [dist(params[i,j]...) for i in YOUNG:OLD, j in YOUNG:OLD]
sample_list = zeros(comparison_samples)
ws_samples = (
- Y = sample(rng, ws_data.Y, comparison_samples),
- M = sample(rng, ws_data.M, comparison_samples),
- O = sample(rng, ws_data.O, comparison_samples),
+ Y = sample(rng, rest_data.Y, comparison_samples),
+ M = sample(rng, rest_data.M, comparison_samples),
+ O = sample(rng, rest_data.O, comparison_samples),
+ )
+ errsum = 0
+ for age_sample in YOUNG:OLD
+ for age_j in YOUNG:OLD #for a given age_sample loop over possible contact ages
+ if neighourhoods[age_sample,age_j] > 0
+ durs = trunc.(Int,rand(rng,age_dists[age_sample,age_j],neighourhoods[age_sample,age_j])) .% durmax
+ # display(durs)
+ tot_dur_sample!(sample_list,cnst_hh.Sparam,durs)
+ kde_est = kde(sample_list)
+ # err = (1 - pvalue(KSampleADTest(ws_samples[age_sample],sample_list)))^2 #need to maximize probability of null hypothesis, not rly valid but everyone does it so idk
+ errsum += mapreduce(+,0:0.05:durmax) do i
+ return (pdf(kde_est,i) - pdf(kerneldensity_data[age_sample],i))^2
+ end
+ end
+ end
+ end
+ return errsum
+end
+
+function err_rest(p,dist)
+ params = get_params(p)
+ neighourhoods = rand.(rng,rest_distributions)
+ age_dists = [dist(params[i,j]...) for i in YOUNG:OLD, j in YOUNG:OLD]
+ sample_list = zeros(comparison_samples)
+ rest_samples = (
+ Y = sample(rng, rest_data.Y, comparison_samples),
+ M = sample(rng, rest_data.M, comparison_samples),
+ O = sample(rng, rest_data.O, comparison_samples),
)
errsum = 0
for age_sample in YOUNG:OLD
@@ -27,7 +68,7 @@ function err_ws(p,dist)
durs = trunc.(Int,rand(rng,age_dists[age_sample,age_j],neighourhoods[age_sample,age_j])) .% durmax
# display(durs)
tot_dur_sample!(sample_list,cnst_hh.Sparam,durs)
- err = 1 - pvalue(KSampleADTest(ws_samples[age_sample],sample_list)) #need to maximize probability of null hypothesis, not rly valid but everyone does it so idk
+ err = 1 - pvalue(KSampleADTest(rest_samples[age_sample],sample_list)) #need to maximize probability of null hypothesis, not rly valid but everyone does it so idk
errsum += err
end
end
diff --git a/ws.pdf b/ws.pdf
deleted file mode 100644
index 0273814a404871a726d8147ad6b4745d5d570619..0000000000000000000000000000000000000000
Binary files a/ws.pdf and /dev/null differ