diff --git a/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl b/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl
index b4de1c4291da542acc829966a4f7dcd0a3becfdb..85623176d324248c70017bdd81f3a68d8c0e7749 100644
--- a/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl
+++ b/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl
@@ -19,19 +19,21 @@ using NetworkLayout:SFDP
import LightGraphs.add_edge!
+const PACKAGE_FOLDER = dirname(dirname(pathof(CovidAlertVaccinationModel)))
+
+const RNG = Xoroshiro128Star()
+
include("agents.jl")
-include("mixing_distributions.jl")
include("data.jl")
+include("mixing_distributions.jl")
include("model.jl")
include("plotting.jl")
include("graphs.jl")
-const RNG = Xoroshiro128Star()
const population = 14.57e6 #population of ontario
const color_palette = palette(:seaborn_pastel) #color theme for the plots
const age_bins = [(0.0, 25.0),(25.0,65.0),(65.0,Inf)]
-const PACKAGE_FOLDER = dirname(dirname(pathof(CovidAlertVaccinationModel)))
const infection_data = parse_cases_data()
const demographic_distribution = get_canada_demographic_distribution()
diff --git a/CovidAlertVaccinationModel/src/data.jl b/CovidAlertVaccinationModel/src/data.jl
index f125cffb81346d8c46416e5c6480be2d2ac90485..65ff2d8935d722d8103b56be299f2eadfd04af9c 100644
--- a/CovidAlertVaccinationModel/src/data.jl
+++ b/CovidAlertVaccinationModel/src/data.jl
@@ -28,9 +28,9 @@ function get_canada_case_fatality()::Tuple{Vector{Tuple{Float64,Float64}},Vector
end
function find_household_composition(df_row)
age_resp_to_bin = Dict(
+ "Y" => 1,
"M" => 2,
"O" => 3,
- "Y" => 1
)
u25_bins = [:U15CHILD,:O15CHILD,:YSPOUSE]
m_bins = [:MPAR, :MCHILD,:MHHADULT]
@@ -51,6 +51,16 @@ function sample_household_data(n)
return sample(RNG,households, Weights(weight_vector),n)
# https://www.publichealthontario.ca/-/media/documents/ncov/epi/covid-19-severe-outcomes-ontario-epi-summary.pdf?la=en
end
+function get_household_data_proportions()
+ f = readdlm(joinpath(PACKAGE_FOLDER,"data/csv/home_compositions.csv"), ',')
+ df = DataFrame([f[1,i] => f[2:end, i] for i = 1:length(f[1,:])])
+ weight_vector::Vector{Float64} = df[!,:WGHT_PER]/sum(df[!,:WGHT_PER])
+ households = map(find_household_composition,eachrow(df))
+ households_by_demographic_sum = sum.([map(l-> l[i], households) for i in 1:3])
+ return households_by_demographic_sum./sum(households_by_demographic_sum)
+ # https://www.publichealthontario.ca/-/media/documents/ncov/epi/covid-19-severe-outcomes-ontario-epi-summary.pdf?la=en
+end
+
function parse_string_as_float_pair(s)::Tuple{Float64,Float64}
parsed = strip(s, ['(',')']) |> s -> split(s, ",")
diff --git a/CovidAlertVaccinationModel/src/mixing_distributions.jl b/CovidAlertVaccinationModel/src/mixing_distributions.jl
index 398209c9d35f9feac13f0f128d36bd1e3e5439b2..f054ef68750769ef6538538bcc34f7181eeb7da1 100644
--- a/CovidAlertVaccinationModel/src/mixing_distributions.jl
+++ b/CovidAlertVaccinationModel/src/mixing_distributions.jl
@@ -63,24 +63,41 @@ function alpha_matrix(alphas)
return M
end
-const alphas = alpha_matrix((
- Y = 0.11,
- M = 0.14,
- O = 0.87
-))
-
-
-const initial_workschool_mixing_matrix = map(t->from_mean(t...),[
- (ZWDist{Geometric{Float64},Discrete}, 0.433835,4.104848) (ZWDist{Geometric{Float64},Discrete},0.406326,2.568782) (ZWDist{Geometric{Float64},Discrete},0.888015,0.017729)
- (ZWDist{Geometric{Float64},Discrete}, 0.600966,0.975688) (ZWDist{Geometric{Float64},Discrete},0.4306,5.057572) (ZWDist{Geometric{Float64},Discrete}, 0.84513,0.021307)
- (ZWDist{Geometric{Float64},Discrete},0.887392,0.001937) (ZWDist{Geometric{Float64},Discrete},0.793004,0.00722) (ZWDist{Geometric{Float64},Discrete},0.940473, 0.022134)
-])
+const unemployment_matrix = alpha_matrix(
+ (
+ Y = 0.11,
+ M = 0.14,
+ O = 0.87
+ )
+)
+
+function make_workschool_mixing_matrix()
+ ws_mixing = map(t->from_mean(t...),[
+ (Geometric{Float64}, 4.104848) (Geometric{Float64},2.568782) (Geometric{Float64},0.017729)
+ (Geometric{Float64}, 0.975688) (Geometric{Float64},5.057572) (Geometric{Float64},0.021307)
+ (Geometric{Float64},0.001937) (Geometric{Float64},0.00722) (Geometric{Float64}, 0.022134)
+ ])
+
+ ws_mixing_w_unemployment_symmetrized = symmetrize_means(get_household_data_proportions(),ws_mixing)
+ ws_adjust_mean(W) = (5/7) .* (1 .- unemployment_matrix) ./ ( mean.(W) + (5/7) .* (1 .- unemployment_matrix))
+ ws_mixing_w_unemployment_symmetrized_weekday_adjusted = ZWDist.(unemployment_matrix,Geometric.(ws_adjust_mean(ws_mixing_w_unemployment_symmetrized)))
+ return ws_mixing_w_unemployment_symmetrized_weekday_adjusted
+end
-const initial_rest_mixing_matrix = map(t->from_mean(t...),[
+const initial_rest_mixing_matrix = symmetrize_means(get_household_data_proportions(),map(t->from_mean(t...),[
(Geometric{Float64},2.728177) (Geometric{Float64},1.382557) (Geometric{Float64},0.206362)
(Geometric{Float64},1.139072) (Geometric{Float64},3.245594) (Geometric{Float64},0.785297)
(Geometric{Float64},0.264822) (Geometric{Float64},0.734856) (Geometric{Float64},0.667099)
-])
+]))
+const initial_workschool_mixing_matrix = make_workschool_mixing_matrix()
+@views function ws_sample(age)
+ return rand.(initial_workschool_mixing_matrix[age,:]) * (rand(RNG) < (5/7))
+end
+
+# ws_adjust_mean(M) = from_mean.(typeof.(M),(7/5) .* mean.(M) ./ (1 .- alphas)) (5/7) * (1-u[i,j]) / { E[X[i,j]] + (5/7) * (1-u[i,j]) }
+# const initial_workschool_mixing_whole_population_symmetrized = ws_adjust_mean(symmetrize_means(get_household_data_proportions(),initial_workschool_mixing_matrix))
+
+# const initial_rest_mixing_matrix_whole_population_symmetrized = symmetrize_means(get_household_data_proportions(),initial_rest_mixing_matrix)
const contact_time_distribution_matrix = [Geometric() for i in 1:(AgentDemographic.size-1), j in 1:(AgentDemographic.size-1)]
diff --git a/IntervalsModel/plots/Distributions.Poisson_hh.pdf b/IntervalsModel/plots/Distributions.Poisson_hh.pdf
index 4249051de9b3a7f0d668193bf4d4925764c255f3..1ba31a7eefe870e7e4072ddd7e266fa00e022d8d 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 b4c5605e624f6801ef17e596452d7ce6e50841c8..59cd700004b933f910fd4effd1523590f4e0fdda 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
new file mode 100644
index 0000000000000000000000000000000000000000..5cd34c7f527494712afc603cc3ccf0ebeed8f2d4
Binary files /dev/null and b/IntervalsModel/plots/Distributions.Poisson_ws.pdf differ
diff --git a/IntervalsModel/plots/hh.pdf b/IntervalsModel/plots/hh.pdf
index 549d2e5f7d39268a26b416ba71874e37b74268b8..faef2d1ad0c79eae6626b480ba089d20f93ae287 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 dd59bed8c51c36322556f8afbc893e568af39eb3..89decf0a94ea8a69fdb091e2efe540b70359061b 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
new file mode 100644
index 0000000000000000000000000000000000000000..243733f078430548209c9e597c56d8cbfa32fb7f
Binary files /dev/null and b/IntervalsModel/plots/ws.pdf differ
diff --git a/IntervalsModel/simulation_data/hh.dat b/IntervalsModel/simulation_data/hh.dat
index f2f9d6e3e86458f084ef5018c3c799492c0f9b4e..eb57911038e6bf2d99523665eb316df80e6f8ec2 100644
Binary files a/IntervalsModel/simulation_data/hh.dat and b/IntervalsModel/simulation_data/hh.dat differ
diff --git a/IntervalsModel/simulation_data/rest.dat b/IntervalsModel/simulation_data/rest.dat
index 0c2050dd5e991dd5f29f85e0e4d3124bd0a802d7..e43cb1058e2f3a3ba27ad57bc4518efe012f4b57 100644
Binary files a/IntervalsModel/simulation_data/rest.dat and b/IntervalsModel/simulation_data/rest.dat differ
diff --git a/IntervalsModel/simulation_data/ws.dat b/IntervalsModel/simulation_data/ws.dat
new file mode 100644
index 0000000000000000000000000000000000000000..687a6b2c3245be85104c1a36ed856e2dcbce8fbe
Binary files /dev/null and b/IntervalsModel/simulation_data/ws.dat differ
diff --git a/IntervalsModel/src/IntervalsModel.jl b/IntervalsModel/src/IntervalsModel.jl
index 6c8a963e7d4b063dc5d2c584c94cdf6c75ac8ef3..e79d487f5afa7320a6bfad7cec5473a5d7f96bc5 100644
--- a/IntervalsModel/src/IntervalsModel.jl
+++ b/IntervalsModel/src/IntervalsModel.jl
@@ -7,6 +7,7 @@ using CSV
using DataFrames
using RandomNumbers.Xorshifts
using Random
+using KernelDensity
using StatsBase
using Distributions
using CovidAlertVaccinationModel
@@ -19,49 +20,75 @@ using BenchmarkTools
using Plots
-const rng = Xoroshiro128Plus(1)
+const rng = Xoroshiro128Plus()
+
+#lets us use nicer names for the indices
const YOUNG, MIDDLE,OLD = 1,2,3
const durmax = 144
-const color_palette = palette(:seaborn_bright) #color theme for the plots
+const comparison_samples = 100
const sparam = [60, 12]
- # Swap age brackets for numbers
+# Swap age brackets for numbers
const swap_dict = OrderedDict("Y" => YOUNG, "M" => MIDDLE, "O" => OLD)
include("utils.jl")
include("interval_overlap_sampling.jl")
-
include("hh_durations_model.jl")
-include("ws_rest_durations_model.jl")
+include("ws_durations_model.jl")
+include("rest_durations_model.jl")
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)
+ plot_all()
+end
-const priors = get_priors()
-const α_priors_mean_ws = (
- Y = 0.364286,
- M = 0.385714,
- O = 0.907143
-)
+"""
+ bayesian_estimation(
+ fname::String,
+ err_func::Function,
+ priors_list::Vector{Distribution},
+ dists::Vector{Distribution},
+ particles::Integer;
+ alpha = 0.95
+ )
-const α_priors_mean_rest = (
- Y = 0.635714,
- M = 0.614286,
- O = 0.092857
-)
+# Arguments
+- fname
+Filename to save output particles. Will be saved in `simulations/` as `fname.dat`. The format is just straight serialized Julia objects (using Serialization.jl), which are neither small nor fast to load, but very convient.
-pgfplotsx()
-function main()
+- err_func(p,distribution)
- do_hh(400)
- # do_ws(400)
- do_rest(1000)
+A function that takes two arguments, a vector of particles with (currently since nonsymmetric) 9n elements for some integer n, and a distribution to fit that takes n parameters.
-end
-using BenchmarkTools
-function bayesian_estimation(fname, err_func, priors_list, dists, particles; alpha = 0.95)
+- priors_list
+
+A list of 9n distributions, for some integer n, that describe the prior distribution of parameters for a distribution with n parameters.
+
+- dists
+
+A list of univariate distributions to fit to the distribution of contact times.
+
+- particles
+
+Number of particles to obtain representing posterior distribution of the distribution parameters.
+
+# Optional Arguments
+
+- alpha = 0.95
+
+Threshold adaptive rate, see `KissABC.smc` for more details.
+"""
+function bayesian_estimation(fname, err_func, priors_list, dists, particles; alpha = 0.99)
+
data_pairs = map(zip(dists,priors_list)) do (dist,priors)
init = rand(priors)
- @btime $err_func($init,$dist)
+ # @btime $err_func($init,$dist)
out = smc(priors,p -> err_func(p, dist), verbose=true, nparticles=particles, alpha=alpha, parallel = true)
return dist => out
end |> OrderedDict
@@ -70,5 +97,52 @@ function bayesian_estimation(fname, err_func, priors_list, dists, particles; alp
serialize(joinpath(PACKAGE_FOLDER,"simulation_data","$fname.dat"),data_pairs)
end
+"""
+Fit HH durations. This function specifies a list of distributions to fit, and loads the priors for "hh", and then calls `bayesian_estimation` with `err_hh`. This function passes parameters to `bayesian_estimation` depending on what produces the best fit for the HH data.
+"""
+function do_hh(particles)
+ dists = [
+ Poisson,
+ ]
+
+ poisson_priors = filter_priors("home")
+ # Set parameter priors for fitting
+ priors_list = [
+ Factored(poisson_priors...),
+ ]
+ bayesian_estimation("hh",err_hh,priors_list,dists,particles; alpha = 0.99)
+end
+
+
+"""
+Fit WS durations. This function specifies a list of distributions to fit, and loads the priors for "workschool", and then calls `bayesian_estimation` with `err_ws`. This function passes parameters to `bayesian_estimation` depending on what produces the best fit for the WS data.
+"""
+function do_ws(particles)
+ dists = [
+ Poisson,
+ ]
+ poisson_priors = filter_priors("workschool")
+ # Set parameter priors for fitting
+ priors_list = [
+ Factored(poisson_priors...)
+ ]
+ bayesian_estimation("ws",err_ws,priors_list,dists,particles; alpha = 0.99)
+end
+
+"""
+Fit rest durations. This function specifies a list of distributions to fit, and loads the priors for "rest", and then calls `bayesian_estimation` with `err_rest`. This function passes parameters to `bayesian_estimation` depending on what produces the best fit for the Rest data.
+"""
+function do_rest(particles)
+ dists = [
+ Poisson,
+ ]
+ poisson_priors = filter_priors("rest")
+ # Set parameter priors for fitting
+ priors_list = [
+ Factored(poisson_priors...)
+ ]
+ bayesian_estimation("rest",err_rest, priors_list, dists, particles, alpha = 0.99)
+end
+
end # module
\ No newline at end of file
diff --git a/IntervalsModel/src/hh_durations_model.jl b/IntervalsModel/src/hh_durations_model.jl
index 124dd11c5dc1787b6bf8f18740b579443b7b5d3f..2f2eebc038a3421939f7942279c68c12fb9cd1f8 100644
--- a/IntervalsModel/src/hh_durations_model.jl
+++ b/IntervalsModel/src/hh_durations_model.jl
@@ -1,30 +1,10 @@
const HHYMO = DataFrame(CSV.File("$PACKAGE_FOLDER/network-data/Timeuse/HH/HHYMO.csv"))
-
const cnst_hh = (
- # Set parameters for intervals sample and subpopulation size
- numsamples = 100,
subsize = size(HHYMO)[1],
- # Total weight in survey
Wghttotal = sum(HHYMO[:,"WGHT_PER"]),
)
-
-function do_hh(particles)
- dists = [
- Poisson,
- ]
-
-
- poisson_priors = filter_priors("rest")
-
- # # Set parameter priors for fitting
- priors_list = [
- Factored(poisson_priors...)
- ]
- fname = "hh"
- bayesian_estimation(fname,err_hh,priors_list,dists,particles; alpha = 0.95)
-end
function make_dat_array()
durs = hcat(
Int.(HHYMO[!,"YDUR"*string(sparam[2])]),
@@ -36,7 +16,6 @@ function make_dat_array()
Int.(HHYMO[!,"MNUM"]),
Int.(HHYMO[!,"ONUM"]),
)
-
WGHT = Weights(HHYMO[!,"WGHT_PER"]./cnst_hh.Wghttotal)
AGERESP = map(r -> swap_dict[r],HHYMO[!,"AGERESP"])
return (;
@@ -65,8 +44,8 @@ function err_hh(p,dist)
age_sample = AGERESP[i]
@inbounds for age_j in YOUNG:OLD #for a given age_sample loop over possible contact ages
durs = trunc.(Int,rand(rng,age_dists[age_sample,age_j],num_contacts_subarray[i,age_j])) .% durmax
- expdur = tot_dur_sample(cnst_hh.numsamples,durs)
- errsum += (expdur/cnst_hh.numsamples - duration_subarray[i,age_j])^2 #compute total
+ expdur = tot_dur_sample(comparison_samples,durs)
+ errsum += (expdur/comparison_samples - duration_subarray[i,age_j])^2 #compute total
end
end
# display("2")
diff --git a/IntervalsModel/src/plotting_functions.jl b/IntervalsModel/src/plotting_functions.jl
index f6b8b1b0fc53213584f3fd71c2e481ee41efde32..f03414f3d936738a01135ca86c29ede1dd3dd2b9 100644
--- a/IntervalsModel/src/plotting_functions.jl
+++ b/IntervalsModel/src/plotting_functions.jl
@@ -1,5 +1,6 @@
using LaTeXStrings
using Plots
+const color_palette = palette(:seaborn_bright) #color theme for the plots
default(dpi = 300)
default(framestyle = :box)
@@ -17,7 +18,6 @@ function plot_estimate(fname)
end
function plot_dists(fname,dist_constructors,data)
function get_non_symmetric_params(p)
-
p_matrix = [[zero(eltype(p))] for i in 1:3, j in 1:3]
for i in 1:9
p_matrix[i][1] = p[i]
@@ -40,7 +40,7 @@ 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 = true,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[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
@@ -57,15 +57,18 @@ function plot_posteriors(fname,data)
p_list = map(x -> plot(),1:length(data.P))
for i in 1:length(data.P)
# display(data.P[i].particles)
- hist = fit(Histogram,data.P[i].particles; nbins = 50)
+ 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)
+ plot!(p_list[i],hist;legend = false, xlabel = L"\lambda_%$i")
# display(kernel_data)
# vline!(p_list[i],[argmax(kernel_data)]; seriescolor = color_palette[2])
end
- p = plot(p_list...;layout=(length(p_list) ÷ 9,9), size = (1500,(length(p_list) ÷ 9)*300), seriescolor = color_palette[1])
+ plot!(p_list[1]; ylabel = "No. of particles")
+ plot!(p_list[4]; ylabel = "No. of particles")
+ plot!(p_list[7]; ylabel = "No. of particles")
+ p = plot(p_list...; size = (800,600), seriescolor = color_palette[1])
savefig(p,joinpath(PACKAGE_FOLDER,"plots","$fname.pdf"))
end
function add_subplot_letters!(plot_list; pos = :top)
diff --git a/IntervalsModel/src/rest_durations_model.jl b/IntervalsModel/src/rest_durations_model.jl
new file mode 100644
index 0000000000000000000000000000000000000000..e40e217e62e6793b16d1f96b66db77f0acc815be
--- /dev/null
+++ b/IntervalsModel/src/rest_durations_model.jl
@@ -0,0 +1,35 @@
+
+const rest_distributions = CovidAlertVaccinationModel.initial_rest_mixing_matrix
+const rest_data = get_rest_data()
+const kerneldensity_data_rest = (
+ Y = kde(rest_data.Y.durs,weights = rest_data.Y.weights),
+ M = kde(rest_data.M.durs,weights = rest_data.M.weights),
+ O = kde(rest_data.O.durs,weights = rest_data.O.weights),
+)
+
+#error function for rest
+function err_rest(p,dist)
+ p_matrix = zeros(3,3)
+ for i in 1:9
+ p_matrix[i] = p[i]
+ end
+
+ neighourhoods = rand.(rng,rest_distributions)
+ age_dists = [dist(p_matrix[i,j]...) for i in YOUNG:OLD, j in YOUNG:OLD]
+ sample_list = zeros(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
+ tot_dur_sample!(sample_list,durs)
+ kde_est = kde(sample_list)
+ errsum += mapreduce(+,0:1:durmax) do i
+ return abs(pdf(kde_est,i) - pdf(kerneldensity_data_ws[age_sample],i))
+ end
+ end
+ end
+ end
+ return errsum
+end
\ No newline at end of file
diff --git a/IntervalsModel/src/utils.jl b/IntervalsModel/src/utils.jl
index 6a8c7c58669dad642b157b13bce358827c4f737e..5215ce85bfafa4ad5e5d8b1fa83a7fb72bd1acff 100644
--- a/IntervalsModel/src/utils.jl
+++ b/IntervalsModel/src/utils.jl
@@ -1,19 +1,39 @@
+import Pandas: read_csv
+using DataFrames
-nparams(d::Type{T}) where T<:Sampleable = length(fieldnames(d))
-nparams(::Type{ZWDist{T,S}}) where {T<:Sampleable,S} = 1+nparams(T)
-
function get_params(params)
p_list = [as_symmetric_matrix(params[i:i+5]) for i = 1:6:length(params)]
return zip(p_list...) |> collect
end
-import Pandas: read_csv
-using DataFrames
+
+function get_rest_data()
+ path = "$PACKAGE_FOLDER/network-data/Timeuse/Rest/RData"
+ data_table_by_age = map(collect(keys(swap_dict))) do age
+ data_table = CSV.File(path*"$age.csv") |> Tables.matrix
+ weights = Weights(data_table[:,1])
+ durs = data_table[:,2]
+ Symbol(age) => (;durs, weights)
+ end
+ return (;data_table_by_age...)
+end
+
+function get_ws_data()
+ path = "$PACKAGE_FOLDER/network-data/Timeuse/WS/WorkschoolData"
+ data_table_by_age = map(collect(keys(swap_dict))) do age
+ data_table = CSV.File(path*"$age.csv") |> Tables.matrix
+ weights = Weights(data_table[:,1])
+ durs = data_table[:,2]
+ Symbol(age) => (;durs, weights)
+ end
+ return (;data_table_by_age...)
+end
+
function get_priors()
df = DataFrame(read_csv("IntervalsModel/network-data/POLYMOD/AALPoisPriors.csv"))
df_w_indicies = transform(
df,
- :Age_in => e->map(x-> swap_dict[x],e),
+ :Age_in => e -> map(x-> swap_dict[x],e),
:Age_out => e -> map(x -> swap_dict[x],e),
:Age_in,:Age_out
)
@@ -22,6 +42,10 @@ function get_priors()
end
function filter_priors(key)
+ priors = get_priors()
+ if !(key in priors[!,:location])
+ throw(ArgumentError("location key not in priors file!"))
+ end
priors_probs = filter(row-> row["location"] == key, eachrow(priors)) |>
df -> map(r -> (r[:index_out],r[:index_in]) => convert(Vector,r[string.(0:143)]),df) |>
Dict |>
diff --git a/IntervalsModel/src/ws_durations_model.jl b/IntervalsModel/src/ws_durations_model.jl
new file mode 100644
index 0000000000000000000000000000000000000000..cafdb6cb82ed50d1ff60fa011df3e9a4236d7d70
--- /dev/null
+++ b/IntervalsModel/src/ws_durations_model.jl
@@ -0,0 +1,34 @@
+
+const ws_distributions = CovidAlertVaccinationModel.initial_workschool_mixing_matrix
+const ws_data = get_ws_data()
+const kerneldensity_data_ws = (
+ Y = kde(ws_data.Y.durs;weights = ws_data.Y.weights),
+ M = kde(ws_data.M.durs;weights = ws_data.M.weights),
+ O = kde(ws_data.O.durs;weights = ws_data.O.weights),
+)
+
+#error function for ws
+function err_ws(p,dist)
+ p_matrix = zeros(3,3)
+ for i in 1:9
+ p_matrix[i] = p[i]
+ end
+ age_dists = [dist(p_matrix[i,j]...) for i in YOUNG:OLD, j in YOUNG:OLD]
+ sample_list = zeros(comparison_samples)
+ errsum = 0
+ for age_sample in YOUNG:OLD
+ neighourhoods = CovidAlertVaccinationModel.ws_sample(age_sample)
+ for age_j in YOUNG:OLD #for a given age_sample loop over possible contact ages
+ if neighourhoods[age_j] > 0
+ durs = trunc.(Int,rand(rng,age_dists[age_sample,age_j],neighourhoods[age_j])) .% durmax
+ tot_dur_sample!(sample_list,durs)
+ kde_est = kde(sample_list)
+
+ errsum += mapreduce(+,0:1:durmax) do i
+ return abs(pdf(kde_est,i) - pdf(kerneldensity_data_ws[age_sample],i))
+ end
+ end
+ end
+ end
+ return errsum
+end
diff --git a/IntervalsModel/src/ws_rest_durations_model.jl b/IntervalsModel/src/ws_rest_durations_model.jl
deleted file mode 100644
index b75164c8c84d2c4992132348b804ffb9788a6ae0..0000000000000000000000000000000000000000
--- a/IntervalsModel/src/ws_rest_durations_model.jl
+++ /dev/null
@@ -1,120 +0,0 @@
-using KernelDensity
-function get_ws_rest_data(path)
- data_table_by_age = map(collect(keys(swap_dict))) do age
- data_table = CSV.File(path*"$age.csv") |> Tables.matrix
- weights = Weights(data_table[:,1])
- durs = data_table[:,2]
- Symbol(age) => (;durs, weights)
- end
- return (;data_table_by_age...)
-end
-
-
-const ws_data = get_ws_rest_data("$PACKAGE_FOLDER/network-data/Timeuse/WS/WorkschoolData")
-const kerneldensity_data_ws = (
- Y = kde(ws_data.Y.durs;weights = ws_data.Y.weights),
- M = kde(ws_data.M.durs;weights = ws_data.M.weights),
- O = kde(ws_data.O.durs;weights = ws_data.O.weights),
-)
-
-const rest_data = get_ws_rest_data("$PACKAGE_FOLDER/network-data/Timeuse/Rest/RData")
-const kerneldensity_data_rest = (
- Y = kde(rest_data.Y.durs,weights = rest_data.Y.weights),
- M = kde(rest_data.M.durs,weights = rest_data.M.weights),
- O = kde(rest_data.O.durs,weights = rest_data.O.weights),
-)
-
-
-const comparison_samples = 5000
-
-ws_distributions = CovidAlertVaccinationModel.initial_workschool_mixing_matrix
-rest_distributions = CovidAlertVaccinationModel.initial_rest_mixing_matrix
-
-
-function do_ws(particles)
- dists = [
- Poisson,
- ]
-
- # Set parameter priors for fitting
- # alpha_priors = [TriangularDist(α*(0.8),min(1,α*(1.2),α)) for α in alpha_matrix(α_priors_mean_ws)]
- # priors_list = map(l -> Factored(vcat(l...)...),[
- # [alpha_priors,[Uniform(μ_bounds...) for i=1:6], [Uniform(σ_bounds...) for i = 1:6]],
- # [alpha_priors,[Uniform(μ_bounds...) for i=1:6]],
- # [alpha_priors,[Uniform(μ_bounds...) for i=1:6], [Uniform(σ_bounds...) for i = 1:6]],
- # ])
- # priors_list = map(l -> Factored(vcat(l...)...),[
- # filter(priors,
- # ])
- bayesian_estimation("ws",err_ws,priors_list,dists,particles; alpha = 0.99)
-end
-
-function do_rest(particles)
- dists = [
- Poisson,
- ]
-
-
- poisson_priors = filter_priors("rest")
-
- # # Set parameter priors for fitting
- priors_list = [
- Factored(poisson_priors...)
- ]
- display(poisson_priors)
- bayesian_estimation("rest",err_rest, priors_list, dists, particles, alpha = 0.98)
-end
-
-
-#error function for ws
-function err_ws(p,dist)
- p_matrix = zeros(3,3)
- for i in 1:9
- p_matrix[i] = p[i]
- end
- neighourhoods = rand.(rng,ws_distributions)
- age_dists = [dist(p_matrix[i,j]...) for i in YOUNG:OLD, j in YOUNG:OLD]
- sample_list = zeros(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
- tot_dur_sample!(sample_list,durs)
- kde_est = kde(sample_list)
-
- errsum += mapreduce(+,0:1:durmax) do i
- return abs(pdf(kde_est,i) - pdf(kerneldensity_data_ws[age_sample],i))
- end
- end
- end
- end
- return errsum
-end
-
-#error function for rest
-function err_rest(p,dist)
- p_matrix = zeros(3,3)
- for i in 1:9
- p_matrix[i] = p[i]
- end
-
- neighourhoods = rand.(rng,rest_distributions)
- age_dists = [dist(p_matrix[i,j]...) for i in YOUNG:OLD, j in YOUNG:OLD]
- sample_list = zeros(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
- tot_dur_sample!(sample_list,durs)
- kde_est = kde(sample_list)
- errsum += mapreduce(+,0:1:durmax) do i
- return abs(pdf(kde_est,i) - pdf(kerneldensity_data_ws[age_sample],i))
- end
- end
- end
- end
- return errsum
-end
\ No newline at end of file
diff --git a/ZeroWeightedDistributions/src/ZeroWeightedDistributions.jl b/ZeroWeightedDistributions/src/ZeroWeightedDistributions.jl
index 81cccddb2218a8752a2dc57ba72a2b812e61c33b..4d2c1abf423280f7a394661f1f29d7d0c5c7abd7 100644
--- a/ZeroWeightedDistributions/src/ZeroWeightedDistributions.jl
+++ b/ZeroWeightedDistributions/src/ZeroWeightedDistributions.jl
@@ -32,13 +32,9 @@ function Distributions.pdf(d::ZWDist, x)
end
end
-# function Distributions.mean(d::ZWDist, x)
-# if x == 0
-# return d.α + (1-d.α)*pdf(d.base_dist,0)
-# else
-# return pdf(d.base_dist,x)
-# end
-# end
+function Distributions.mean(d::ZWDist, x)
+ return (1-d.α)*StatsBase.mean(d.base_dist,0)
+end
function Base.rand(rng::AbstractRNG, s::ZWDist{DType,S}, n::Int) where {DType, S}