diff --git a/IntervalsModel/plots/Distributions.Poisson_hh.pdf b/IntervalsModel/plots/Distributions.Poisson_hh.pdf
index 1ba31a7eefe870e7e4072ddd7e266fa00e022d8d..76f28596b32d2dff18b708363c1bbd952704ec3d 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 59cd700004b933f910fd4effd1523590f4e0fdda..04bcddc7c3dcc3c07b4880734b7dfed3d27c013c 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 5cd34c7f527494712afc603cc3ccf0ebeed8f2d4..8e09fe4b6c6c16b1182c5ae8cc4db8a63b00831c 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 faef2d1ad0c79eae6626b480ba089d20f93ae287..fa2572fa0a1315ff43a946aadaa07b71c4ca7b8a 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 89decf0a94ea8a69fdb091e2efe540b70359061b..9f8728a9d14cf1cc15f8dce54ad848a05f4ed09a 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 243733f078430548209c9e597c56d8cbfa32fb7f..1e36f2f154e101e73e18e7d0d4b744ed260b4bce 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
deleted file mode 100644
index fea7ea8ac0d3dd37a0ab30ff5ce899e8d56515ce..0000000000000000000000000000000000000000
Binary files a/IntervalsModel/simulation_data/hh.dat and /dev/null differ
diff --git a/IntervalsModel/simulation_data/rest.dat b/IntervalsModel/simulation_data/rest.dat
index e43cb1058e2f3a3ba27ad57bc4518efe012f4b57..41b8f3dcb9ec866d5d1743c12a56e3d88ed46de6 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
index ee07dd956d8f31c12c7c6005b403b532cdc55b5b..a20eeb8fe50afeeac37ba5c0332ec3777938896b 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 5d730605994b10c4473c7a13e065a05eae9d483e..51275b5cda728d45e2a1abc755c7cd4231f06eb8 100644
--- a/IntervalsModel/src/IntervalsModel.jl
+++ b/IntervalsModel/src/IntervalsModel.jl
@@ -25,11 +25,11 @@ 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
"""
Number of start times to sample for a given set of distribution parameters.
"""
-const comparison_samples = 70
+const comparison_samples = 100
const sparam = [60, 12]
# Swap age brackets for numbers
const swap_dict = OrderedDict("Y" => YOUNG, "M" => MIDDLE, "O" => OLD)
@@ -44,7 +44,7 @@ include("plotting_functions.jl")
Runs parameter estimation for the three scenarios, hh, ws, and rest.
"""
function main()
- do_hh(400)
+ # do_hh(400)
do_ws(400)
do_rest(400)
plot_all()
@@ -131,6 +131,7 @@ function do_ws(particles)
Poisson,
]
poisson_priors = filter_priors("workschool")
+ # display(poisson_priors)
# Set parameter priors for fitting
priors_list = [
Factored(poisson_priors...)
diff --git a/IntervalsModel/src/hh_durations_model.jl b/IntervalsModel/src/hh_durations_model.jl
index fe80396c424e430bbef1cb03f577b171a5f36180..18f65bffb5aaa4a2a91eb0ded9ef691c6038c37a 100644
--- a/IntervalsModel/src/hh_durations_model.jl
+++ b/IntervalsModel/src/hh_durations_model.jl
@@ -13,7 +13,7 @@ function make_dat_array()
Int.(HHYMO[!,"MNUM"]),
Int.(HHYMO[!,"ONUM"]),
)
- WGHT = Weights(HHYMO[!,"WGHT_PER"]./cnst_hh.Wghttotal)
+ WGHT = Weights(HHYMO[!,"WGHT_PER"]./Wghttotal)
AGERESP = map(r -> swap_dict[r],HHYMO[!,"AGERESP"])
return (;
nums,
@@ -31,10 +31,7 @@ const dat = make_dat_array() #assign a constant data array
Calculate error between observed total contact durations in HH data and simulated total contact durations based on `dist(p[i]...)`. This one works a bit differently than the WS or rest error function, since we do not have a distribution of total durations to compare to, but rather one average total duration, and we do have the explicit samples of numbers of contacts.
"""
function err_hh(p,dist)
- p_matrix = zeros(3,3)
- for i in 1:9
- p_matrix[i] = p[i]
- end
+ p_matrix = get_params(p)
age_dists = [dist(p_matrix[i,j]...) for i in YOUNG:OLD, j in YOUNG:OLD]
duration_subarray = dat.durs #durations, these are call subarray because we used to only sample a subset
num_contacts_subarray = dat.nums #numbers of contacts for each duration
diff --git a/IntervalsModel/src/interval_overlap_sampling.jl b/IntervalsModel/src/interval_overlap_sampling.jl
index c64d83f47cd94223ba51fa1bad083b9cfaeb5ed0..3684fc5520a1b3f14f9c42eec8b1196983ed80de 100644
--- a/IntervalsModel/src/interval_overlap_sampling.jl
+++ b/IntervalsModel/src/interval_overlap_sampling.jl
@@ -71,4 +71,4 @@ function tot_dur_sample!(sample_list,durlist)
union!(int_list)
sample_list[i] = mapreduce(Intervals.span,+,int_list)
end
-end
+end
\ No newline at end of file
diff --git a/IntervalsModel/src/plotting_functions.jl b/IntervalsModel/src/plotting_functions.jl
index f03414f3d936738a01135ca86c29ede1dd3dd2b9..1e85e408c16044a3ed65de8f2a73e10c81c8861f 100644
--- a/IntervalsModel/src/plotting_functions.jl
+++ b/IntervalsModel/src/plotting_functions.jl
@@ -17,14 +17,8 @@ function plot_estimate(fname)
end
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]
- end
- return p_matrix
- end
- p_estimate_as_arrays = map(d -> get_non_symmetric_params(d.P),data)
+
+ p_estimate_as_arrays = map(d -> get_params(d.P),data)
p_matrix = map(x -> plot(),p_estimate_as_arrays[1])
ymo = ["Y","M","O"]
x_range = 0.0:144.0
@@ -67,7 +61,6 @@ function plot_posteriors(fname,data)
end
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
diff --git a/IntervalsModel/src/rest_durations_model.jl b/IntervalsModel/src/rest_durations_model.jl
index 53678dfeca5afda3b305a757e2574d7f7b19ce7b..37ce4e9b2ea0f82bd8462757d34c293b60df2073 100644
--- a/IntervalsModel/src/rest_durations_model.jl
+++ b/IntervalsModel/src/rest_durations_model.jl
@@ -12,15 +12,9 @@ const kerneldensity_data_rest = (
Calculate error between observed contact durations in rest data and simulated contact durations based on `dist(p[i]...)`.
"""
function err_rest(p,dist)
- #we get the vector p as a 9 element vector from KissABC, but more convenient as a 3x3 matrix
- p_matrix = zeros(3,3)
- for i in 1:9
- p_matrix[i] = p[i]
- end
+ #we get the vector p as a 6n element vector from KissABC, but more convenient as a 3x3 matrix of tuples
+ p_matrix = get_params(p)
- #sample a matrix of neighourhood sizes for distributions in rest_distributions
- #this is actually not correct? neighbourhood sizes should be symmetric!
- neighourhoods = rand.(rng,rest_distributions)
#create a matrix of distributions of type dist, from the parameters in p_matrix
age_dists = [dist(p_matrix[i,j]...) for i in YOUNG:OLD, j in YOUNG:OLD]
@@ -29,24 +23,30 @@ function err_rest(p,dist)
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
- #get durations from our candidate distribtions for each of the contacts in neighbourhood
- durs = trunc.(Int,rand(rng,age_dists[age_sample,age_j],neighourhoods[age_sample,age_j])) .% durmax
-
- #this MODIFIES sample_list to contain samples from the distribution of total_durations, given intervals of length dur.
- tot_dur_sample!(sample_list,durs)
-
- #compute a kernel density from the list of samples
- kde_est = kde(sample_list)
+ for _ in 1:sample_repeat
+
+ #sample a matrix of neighourhood sizes for distributions in rest_distributions
+ #this is actually not correct? neighbourhood sizes should be symmetric!
+ neighourhoods = rand.(rng,rest_distributions)
+ 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
+ #get durations from our candidate distribtions for each of the contacts in neighbourhood
+ durs = trunc.(Int,rand(rng,age_dists[age_sample,age_j],neighourhoods[age_sample,age_j])) .% durmax
+
+ #this MODIFIES sample_list to contain samples from the distribution of total_durations, given intervals of length dur.
+ tot_dur_sample!(sample_list,durs)
+
+ #compute a kernel density from the list of samples
+ kde_est = kde(sample_list)
- #compute L1 distance between observed distribution of total durations from survey, and observed distribution from simulation
- errsum += mapreduce(+,0:1:durmax) do i
- return abs(pdf(kde_est,i) - pdf(kerneldensity_data_ws[age_sample],i))
+ #compute L1 distance between observed distribution of total durations from survey, and observed distribution from simulation
+ errsum += mapreduce(+,0:1:durmax) do i
+ return abs(pdf(kde_est,i) - pdf(kerneldensity_data_ws[age_sample],i))
+ end
end
end
end
end
- return errsum
+ return errsum./sample_repeat
end
\ No newline at end of file
diff --git a/IntervalsModel/src/utils.jl b/IntervalsModel/src/utils.jl
index a02d85dad890a3a0ebc7df138be723d79477c2c5..843730bebff44c5998352a409af181342b354cce 100644
--- a/IntervalsModel/src/utils.jl
+++ b/IntervalsModel/src/utils.jl
@@ -24,6 +24,14 @@ function as_symmetric_matrix(l)
]
end
+"""
+ Turn a symmetric 3x3 matrix, l, into a vector of length 6, probably a nicer way to do this exists but meh.
+"""
+function symmetric_matrix_as_vector(A)
+ return [A[1,1],A[1,2],A[1,3],A[2,2],A[2,3],A[3,3]]
+end
+
+
"""
Load rest data from `network-data/Timeuse/Rest/RData`.
"""
@@ -82,13 +90,14 @@ function filter_priors(location_key)
end
priors_probs = filter(row-> row["location"] == location_key, eachrow(priors)) |>
- df -> map(r -> (r[:index_out],r[:index_in]) => convert(Vector,r[string.(0:143)]),df) |>
+ df -> map(r -> (r[:index_out],r[:index_in]) => convert(Vector,r[string.(0:94)]),df) |> #this is pretty inflexible to changing support unfortunately
Dict |>
df -> Dict(map(v -> v => Categorical(df[v]), collect(keys(df))))
output_array = Array{eltype(values(priors_probs)),2}(undef,(3,3))
for (key,value) in priors_probs
output_array[key...] = value
+ output_array[reverse(key)...] = value #corny way to get symmetry assuming key is a 2-tuple
end
- return output_array
+ return symmetric_matrix_as_vector(output_array)
end
\ No newline at end of file
diff --git a/IntervalsModel/src/ws_durations_model.jl b/IntervalsModel/src/ws_durations_model.jl
index 91aadca465abdd4cec66f14883598903cf57c5ff..1f6b81362225b1fe273c68f866be217736cea58e 100644
--- a/IntervalsModel/src/ws_durations_model.jl
+++ b/IntervalsModel/src/ws_durations_model.jl
@@ -16,10 +16,8 @@ Calculate error between observed contact durations in workschool data and simula
"""
function err_ws(p,dist)
#we get the vector p as a 9 element vector from KissABC, but more convenient as a 3x3 matrix
- p_matrix = zeros(3,3)
- for i in 1:9
- p_matrix[i] = p[i] #this will fill p_matrix with the elements of p in column-major order
- end
+
+ p_matrix = get_params(p)
#create a matrix of distributions of type dist, from the parameters in p_matrix
age_dists = [dist(p_matrix[i,j]...) for i in YOUNG:OLD, j in YOUNG:OLD]
@@ -27,28 +25,29 @@ function err_ws(p,dist)
#initialize the list of samples and the total error
sample_list = zeros(comparison_samples)
errsum = 0
-
- for age_sample in YOUNG:OLD
- #Sample neighbourhood sizes. Note that this includes the weekend/weekday coinflip, implemented in CovidAlertVaccinationModel.
- 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
- #get durations from our candidate distribtions for each of the contacts in neighbourhood
- durs = trunc.(Int,rand(rng,age_dists[age_sample,age_j],neighourhoods[age_j])) .% durmax
-
- #this MODIFIES sample_list to contain samples from the distribution of total_durations, given intervals of length dur.
- tot_dur_sample!(sample_list,durs)
-
- #compute a kernel density from the list of samples
- kde_est = kde(sample_list)
-
- #compute L1 distance between observed distribution of total durations from survey, and observed distribution from simulation
- errsum += mapreduce(+,0:1:durmax) do i
- return abs(pdf(kde_est,i) - pdf(kerneldensity_data_ws[age_sample],i))
+ for _ in 1:sample_repeat
+ for age_sample in YOUNG:OLD
+ #Sample neighbourhood sizes. Note that this includes the weekend/weekday coinflip, implemented in CovidAlertVaccinationModel.
+ 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
+ #get durations from our candidate distribtions for each of the contacts in neighbourhood
+ durs = trunc.(Int,rand(rng,age_dists[age_sample,age_j],neighourhoods[age_j])) .% durmax
+
+ #this MODIFIES sample_list to contain samples from the distribution of total_durations, given intervals of length dur.
+ tot_dur_sample!(sample_list,durs)
+
+ #compute a kernel density from the list of samples
+ kde_est = kde(sample_list)
+
+ #compute L1 distance between observed distribution of total durations from survey, and observed distribution from simulation
+ errsum += mapreduce(+,0:1:durmax) do i
+ return abs(pdf(kde_est,i) - pdf(kerneldensity_data_ws[age_sample],i))
+ end
end
end
end
end
- return errsum
+ return errsum./sample_repeat
end