diff --git a/IntervalsModel/plots/Distributions.Normal_hh.pdf b/IntervalsModel/plots/Distributions.Normal_hh.pdf
index 03834507d1a5193d06e61b642bbaa6cb379b1e47..9a0bd8d323e91d704de87f9a254dd5075680dc0b 100644
Binary files a/IntervalsModel/plots/Distributions.Normal_hh.pdf and b/IntervalsModel/plots/Distributions.Normal_hh.pdf differ
diff --git a/IntervalsModel/plots/Distributions.Poisson_hh.pdf b/IntervalsModel/plots/Distributions.Poisson_hh.pdf
index 6536eaa2168a25add30348ea4edfd1f8a309854c..b482737f8ebd79e9cfc6e25558072001a4f479d0 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.Weibull_hh.pdf b/IntervalsModel/plots/Distributions.Weibull_hh.pdf
deleted file mode 100644
index ef42347bce28a7f3b52c89e142fad3b234ce91dc..0000000000000000000000000000000000000000
Binary files a/IntervalsModel/plots/Distributions.Weibull_hh.pdf and /dev/null differ
diff --git a/IntervalsModel/plots/Normal_hh.pdf b/IntervalsModel/plots/Normal_hh.pdf
deleted file mode 100644
index 6849eff0dd197a2cc7228e68465da3152df7fa54..0000000000000000000000000000000000000000
Binary files a/IntervalsModel/plots/Normal_hh.pdf and /dev/null differ
diff --git a/IntervalsModel/plots/Poisson_hh.pdf b/IntervalsModel/plots/Poisson_hh.pdf
deleted file mode 100644
index b7b6ad9d8b74484f5f179ff317a8930e5d32d919..0000000000000000000000000000000000000000
Binary files a/IntervalsModel/plots/Poisson_hh.pdf and /dev/null differ
diff --git a/IntervalsModel/plots/Weibull_hh.pdf b/IntervalsModel/plots/Weibull_hh.pdf
deleted file mode 100644
index d62a6de9361a715b460deb88d7d760b1323b9078..0000000000000000000000000000000000000000
Binary files a/IntervalsModel/plots/Weibull_hh.pdf and /dev/null differ
diff --git a/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Normal, T} where T_rest.pdf b/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Normal, T} where T_rest.pdf
index b77d4c633631d67b7ecd8aad0fc8b47105a00cee..21dea6b8c769ca8f3495f044c78c472399126a9a 100644
Binary files a/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Normal, T} where T_rest.pdf and b/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Normal, T} where T_rest.pdf differ
diff --git a/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Normal, T} where T_ws.pdf b/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Normal, T} where T_ws.pdf
index 065bcc665f35364c93efcb273c037daf4f7228e0..4f8481018141d1c8eef4a919f33b789a8318d6c3 100644
Binary files a/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Normal, T} where T_ws.pdf and b/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Normal, T} where T_ws.pdf differ
diff --git a/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Poisson, T} where T_rest.pdf b/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Poisson, T} where T_rest.pdf
index baf11d19f9a3df84cdce6e28795ad48eb12779d8..2cd1e4c3ad4c2fab5789bdde6c1b9a518314b67a 100644
Binary files a/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Poisson, T} where T_rest.pdf and b/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Poisson, T} where T_rest.pdf differ
diff --git a/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Poisson, T} where T_ws.pdf b/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Poisson, T} where T_ws.pdf
index 635e3cf3bfed0c613b4b4b13123f76f4c9423d1d..5b3d72c516bae8fa679fa36fbc4e2425aa2a7c62 100644
Binary files a/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Poisson, T} where T_ws.pdf and b/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Poisson, T} where T_ws.pdf differ
diff --git a/IntervalsModel/plots/hh.pdf b/IntervalsModel/plots/hh.pdf
index 1c4db69e65a7e839320190ee9ce6e6bd983efafe..15e3ec080369d64a46f916bb4d96e4c9fb1c2310 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 efafae5c8c8ff313f186fbe4260dc9f3ccd410ef..b4d005643cbafc81ec35b9bd461b7179329ce2c6 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 0aac651f1f01fd3763b0e5480bc229e8a86d313d..f032c78bcb18419046082e01705878913a20ab84 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
index 8f24246423c0ecfc005dbc78fe77ef6188a5de09..138e35aabd71c8d03026b4d69a0d91605f0fe3df 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 072fb9a31dccec542b8d238ae79f9b2e253efeb7..bf71a13ce03e5fe4bee757535808ec23d8e3fba5 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 e098055573cf81c847b94ebd2e523e193716eebf..759b4125ffb46ca276dbf645a80ade2c0bcf34a8 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 09df7cd8a317724fe66a69515276d159c9660550..3fb843da0724df325977096dea40cb5292ddd8cb 100644
--- a/IntervalsModel/src/IntervalsModel.jl
+++ b/IntervalsModel/src/IntervalsModel.jl
@@ -1,4 +1,6 @@
 module IntervalsModel
+export main
+
 
 using Intervals
 using CSV
@@ -11,25 +13,26 @@ using ZeroWeightedDistributions
 const PACKAGE_FOLDER = dirname(dirname(pathof(IntervalsModel)))
 using DataStructures:OrderedDict
 using Serialization
-
+using KissABC
+using BenchmarkTools
 using Plots
+
+
 const rng = Xoroshiro128Plus(1)
 const YOUNG, MIDDLE,OLD = 1,2,3
 const durmax = 144
 const color_palette = palette(:seaborn_bright) #color theme for the plots
+const sparam = [60, 12]
 
+    # Swap age brackets for numbers
+const swap_dict = OrderedDict("Y" => YOUNG, "M" => MIDDLE, "O" => OLD)
 include("interval_overlap_sampling.jl")
 include("utils.jl")
 include("hh_durations_model.jl")
-include("ws_durations_model.jl")
+include("ws_rest_durations_model.jl")
 include("plotting_functions.jl")
 
 
-using KissABC
-using BenchmarkTools
-using JLSO
-using Plots
-
 
 const μ_bounds = (1,144) 
 const σ_bounds = (1,144)
@@ -46,119 +49,32 @@ const α_priors_mean_rest = (
     O = 0.092857
 )
 
-function alpha_matrix(alphas)
-    M = zeros(length(alphas),length(alphas))
-    for i in 1:length(alphas), j in 1:length(alphas)
-        M[i,j] = alphas[i] + alphas[j] - alphas[j]*alphas[i]
-    end
-    return [M[1,1], M[1,2],M[1,3],M[2,2],M[2,3],M[3,3]] #lol
-end
 pgfplotsx()
-function do_hh(particles)
-    dists = [
-        Normal,
-        Poisson,
-        Weibull
-    ]
-    bounds_list =  map(l -> vcat(l...),[
-        ([μ_bounds for i = 1:6], [σ_bounds for i = 1:6]),
-        ([μ_bounds for i = 1:6],),
-        ([(1.0,144)  for i = 1:6], [(1.0,144.0) for i = 1:6]),
-    ])
-    @show bounds_list
-    fname = "hh"
-    data_pairs = map(zip(dists,bounds_list)) do (dist,bounds)
-        priors = Factored([Uniform(l,u) for (l,u) in bounds[1:end]]...) #assume uniform priors
-
-        # @btime err_ws($init,$dist) #compute benchmark of the error function, not rly necessary
-        
-        out = smc(priors,p -> err_hh(p, dist), verbose=true, nparticles=particles,  alpha=0.98,parallel = true)
-        
-        return dist => out
-    end |> OrderedDict
-
-    display(data_pairs)
-    serialize(joinpath(PACKAGE_FOLDER,"simulation_data","$fname.dat"),data_pairs)
+function main()
+    do_hh(400)
+    do_ws(400)
+    do_rest(400)
 
 end
 
-function do_ws(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],[(0.1,144)  for i = 1:6], [(0.1,144.0) for i = 1:6]],
-    ])
-    fname = "ws"
-    data_pairs = map(zip(dists,bounds_list)) do (dist,bounds)
-        init = [b[1] for b in bounds]
-        priors = Factored([TriangularDist(μ*(0.8),min(1,μ*(1.2),μ)) for μ in alpha_matrix(α_priors_mean_ws)]...,[Uniform(l,u) for (l,u) in bounds[7:end]]...) #assume uniform priors
-
-        # @btime err_ws($init,$dist) #compute benchmark of the error function, not rly necessary
-        
-        out = smc(priors,p -> err_ws(p, dist), verbose=true, nparticles=particles,  alpha=0.995,  parallel = true)
+function bayesian_estimation(fname, err_func, priors_list, dists, particles; alpha = 0.995)
+    data_pairs = map(zip(dists,priors_list)) do (dist,priors)
+        out = smc(priors,p -> err_ws(p, dist), verbose=true, nparticles=particles,  alpha=alpha,  parallel = true)
         
         return dist => out
     end |> OrderedDict
 
     display(data_pairs)
     serialize(joinpath(PACKAGE_FOLDER,"simulation_data","$fname.dat"),data_pairs)
-
 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],[(0.1,144)  for i = 1:6], [(0.1,144.0) for i = 1:6]],
-    ])
-    # @show bounds_list
-    fname = "rest"
-    data_pairs = map(zip(dists,bounds_list)) do (dist,bounds)
-        init = [b[1] for b in bounds]
-        priors = Factored([TriangularDist(μ*(0.8),min(1,μ*(1.2),μ)) for μ in alpha_matrix(α_priors_mean_rest)]...,[Uniform(l,u) for (l,u) in bounds[7:end]]...) #assume uniform priors
-
-        # @btime err_ws($init,$dist) #compute benchmark of the error function, not rly necessary
-        
-        out = smc(priors,p -> err_rest(p, dist), verbose=true, nparticles=particles,  alpha=0.995,   parallel = true) #apply sequential monte carlo with 200 particles
-        
-        return dist => out
-    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(keys(data)))
-    for (k,v) in zip(keys(data),values(data)) 
-    plot_posteriors("$(k)_ws",v)
-    end
-    data = deserialize(joinpath(PACKAGE_FOLDER,"simulation_data","hh.dat"))
-    plot_dists("hh",collect(keys(data)),collect(values(data)))
-    for (k,v) in zip(keys(data),values(data)) 
-        plot_posteriors("$(k)_hh",v)
-    end
 
-    data = deserialize(joinpath(PACKAGE_FOLDER,"simulation_data","rest.dat"))
-    plot_dists("rest",collect(keys(data)),collect(values(data)))
-    for (k,v) in zip(keys(data),values(data)) 
-        plot_posteriors("$(k)_rest",v)
+function alpha_matrix(alphas)
+    M = zeros(length(alphas),length(alphas))
+    for i in 1:length(alphas), j in 1:length(alphas)
+        M[i,j] = alphas[i] + alphas[j] - alphas[j]*alphas[i]
     end
+    return [M[1,1], M[1,2],M[1,3],M[2,2],M[2,3],M[3,3]] #lol
 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 4b44e1c2d7c946ea572c6c7f9d100984bbe840e9..dfa9d07f77639b575512e1c7b256e6d68e2ed2b0 100644
--- a/IntervalsModel/src/hh_durations_model.jl
+++ b/IntervalsModel/src/hh_durations_model.jl
@@ -5,21 +5,31 @@ const HHYMO = DataFrame(CSV.File("$PACKAGE_FOLDER/network-data/Timeuse/HH/HHYMO.
 # In particular, we avoid having to modify any strings in the error function.
 
 const cnst_hh = (
-    # Set the underlying parameters for the intervals model
-    Sparam = [60,12],
     # Set parameters for intervals sample and subpopulation size
-    numsamples = 10,
+    numsamples = 10_000,
     subsize = size(HHYMO)[1],
-    # Swap age brackets for numbers
-    swap = OrderedDict("Y" => YOUNG, "M" => MIDDLE, "O" => OLD),
     # Total weight in survey
     Wghttotal = sum(HHYMO[:,"WGHT_PER"]),
 )
+
+
+function do_hh(particles)
+    dists = [
+        Normal,
+        Poisson,
+    ]
+    priors_list = map(l -> Factored(vcat(l...)...),[
+        [[Uniform(μ_bounds...) for i=1:6], [Uniform(σ_bounds...)  for i = 1:6]],
+        [[Uniform(μ_bounds...) for i=1:6]],
+    ])
+    fname = "hh"
+    bayesian_estimation(fname,err_hh,priors_list,dists,particles)
+end
 function make_dat_array()
     durs = hcat(
-        Int.(HHYMO[!,"YDUR"*string(cnst_hh.Sparam[2])]),
-        Int.(HHYMO[!,"MDUR"*string(cnst_hh.Sparam[2])]),
-        Int.(HHYMO[!,"ODUR"*string(cnst_hh.Sparam[2])]),
+        Int.(HHYMO[!,"YDUR"*string(sparam[2])]),
+        Int.(HHYMO[!,"MDUR"*string(sparam[2])]),
+        Int.(HHYMO[!,"ODUR"*string(sparam[2])]),
     )
     nums = hcat(
         Int.(HHYMO[!,"YNUM"]),
@@ -28,7 +38,7 @@ function make_dat_array()
     )
 
     WGHT = Weights(HHYMO[!,"WGHT_PER"]./cnst_hh.Wghttotal)
-    AGERESP = map(r -> cnst_hh.swap[r],HHYMO[!,"AGERESP"])
+    AGERESP = map(r -> swap_dict[r],HHYMO[!,"AGERESP"])
     return (;
         nums,
         durs,
@@ -53,7 +63,7 @@ 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,cnst_hh.Sparam,durs)
+            expdur = tot_dur_sample(cnst_hh.numsamples,durs)
             errsum += (expdur/cnst_hh.numsamples - duration_subarray[i,age_j])^2 #compute total 
         end
     end
diff --git a/IntervalsModel/src/interval_overlap_sampling.jl b/IntervalsModel/src/interval_overlap_sampling.jl
index 18227477f2a414c0b7a57121c200fca8ef99f425..8918004c31bec9cece1d5febe62aa49ac453e239 100644
--- a/IntervalsModel/src/interval_overlap_sampling.jl
+++ b/IntervalsModel/src/interval_overlap_sampling.jl
@@ -1,3 +1,6 @@
+# Set the underlying parameters for the intervals model
+const startdist = Normal(sparam...)
+
 
 function coverage!(cov,S_j,E_j)
     if E_j < S_j
@@ -8,7 +11,7 @@ function coverage!(cov,S_j,E_j)
     end
 end
 #compute the total duration of a sample of intervals
-function tot_dur_sample(n, dist,durlist)
+function tot_dur_sample(n,durlist)
     if isempty(durlist)
         return 0
     end
@@ -18,7 +21,7 @@ function tot_dur_sample(n, dist,durlist)
     int_list = Vector{Interval{Int,Closed,Closed}}()
     sizehint!(int_list,numcontact*2)
 
-    start_matrix = trunc.(Int,(rand(rng,dist,(numcontact,n))))
+    start_matrix = trunc.(Int,rand(rng,startdist,(numcontact,n)))
     @inbounds for i in 1:n  
         empty!(int_list)
         @inbounds for j in 1:numcontact
@@ -31,7 +34,9 @@ function tot_dur_sample(n, dist,durlist)
     end
     return total_dur
 end
-function tot_dur_sample!(sample_list, dist,durlist)
+
+
+function tot_dur_sample!(sample_list,durlist)
     if isempty(durlist)
         sample_list .= 0.0
         return
@@ -40,7 +45,8 @@ function tot_dur_sample!(sample_list, dist,durlist)
     n = length(sample_list)
     int_list = Vector{Interval{Int,Closed,Closed}}()
     sizehint!(int_list,numcontact*2)
-    start_matrix = trunc.(Int,(rand(rng,dist,(numcontact,n))))
+    # @show rand(rng,startdist,(numcontact,n))
+    start_matrix = trunc.(Int,rand(rng,startdist,(numcontact,n)))
     for i in 1:n  
         empty!(int_list)
         for j in 1:numcontact
diff --git a/IntervalsModel/src/plotting_functions.jl b/IntervalsModel/src/plotting_functions.jl
index cd2d4f67564b11c9c4fced92362506ffcdbead6f..a1408eb97838a65f4bb2b7761c2052222feb77dd 100644
--- a/IntervalsModel/src/plotting_functions.jl
+++ b/IntervalsModel/src/plotting_functions.jl
@@ -4,6 +4,17 @@ using Plots
 default(dpi = 300)
 default(framestyle = :box)
 pgfplotsx()
+function plot_all()
+    fnames = ["hh","ws","rest"]
+    map(plot_estimate,fnames)
+end
+function plot_estimate(fname)
+    data = deserialize(joinpath(PACKAGE_FOLDER,"simulation_data","$fname.dat"))
+    plot_dists("$fname",collect(keys(data)),collect(values(data)))
+    for (k,v) in zip(keys(data),values(data)) 
+    plot_posteriors("$(k)_$fname",v)
+    end
+end
 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])
diff --git a/IntervalsModel/src/ws_durations_model.jl b/IntervalsModel/src/ws_rest_durations_model.jl
similarity index 66%
rename from IntervalsModel/src/ws_durations_model.jl
rename to IntervalsModel/src/ws_rest_durations_model.jl
index a1f8722192973d69882782b8dac4daf957b62fc4..71b929853da2b231d8381c9239b7bcb04823c071 100644
--- a/IntervalsModel/src/ws_durations_model.jl
+++ b/IntervalsModel/src/ws_rest_durations_model.jl
@@ -1,18 +1,16 @@
-
-using HypothesisTests
-
+using KernelDensity
 function get_ws_rest_data(path)
-    data_table_by_age = map(collect(keys(cnst_hh.swap))) do age
+    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")
-using KernelDensity
 const kerneldensity_data_ws = (
     Y = kde(ws_data.Y.durs),
     M = kde(ws_data.M.durs),
@@ -30,10 +28,43 @@ const kerneldensity_data_rest = (
 const comparison_samples = 5000
 
 ws_distributions = CovidAlertVaccinationModel.initial_workschool_mixing_matrix
-
 rest_distributions = CovidAlertVaccinationModel.initial_rest_mixing_matrix
 
-#error function
+
+function do_ws(particles)
+    dists = [
+        ZWDist{Normal},
+        ZWDist{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]],
+    ])
+    bayesian_estimation("ws",err_ws,priors_list,dists,particles; alpha = 0.99)
+end
+
+function do_rest(particles)
+    dists = [
+        ZWDist{Normal},
+        ZWDist{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]],
+    ])
+
+    bayesian_estimation("rest",err_rest, priors_list, dists, particles, alpha = 0.99)
+end
+
+
+#error function for ws
 function err_ws(p,dist)
     params = get_params(p)
     neighourhoods = rand.(rng,ws_distributions)
@@ -49,20 +80,18 @@ function err_ws(p,dist)
         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)
+                tot_dur_sample!(sample_list,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:1:durmax) do i
                     return abs(pdf(kde_est,i) - pdf(kerneldensity_data_ws[age_sample],i))
                 end
-                # errsum += err#(mean(sample_list) - mean(ws_samples[age_sample]))^2
             end
         end
     end
     return errsum
 end
 
+#error function for rest
 function err_rest(p,dist)
     params = get_params(p)
     neighourhoods = rand.(rng,rest_distributions)
@@ -78,12 +107,8 @@ function err_rest(p,dist)
         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)
-
+                tot_dur_sample!(sample_list,durs)
                 kde_est = kde(sample_list)
-                # 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
                 errsum += mapreduce(+,0:1:durmax) do i
                     return abs(pdf(kde_est,i) - pdf(kerneldensity_data_ws[age_sample],i))
                 end
diff --git a/Manifest.toml b/Manifest.toml
index 976d564f40ea7f42d0f9812717f91ea87ce5d799..327d5bd1fdf0159ffed57efae15b1676c7fd9d19 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -1,5 +1,17 @@
 # This file is machine-generated - editing it directly is not advised
 
+[[AbstractFFTs]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "8ed9de2f1b1a9b1dee48582ad477c6e67b83eb2c"
+uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
+version = "1.0.0"
+
+[[Adapt]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "ffcfa2d345aaee0ef3d8346a073d5dd03c983ebe"
+uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+version = "3.2.0"
+
 [[ArgTools]]
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
 
@@ -12,21 +24,88 @@ version = "3.1.1"
 [[Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 
+[[AxisAlgorithms]]
+deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"]
+git-tree-sha1 = "a4d07a1c313392a77042855df46c5f534076fab9"
+uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950"
+version = "1.0.0"
+
 [[Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
+[[ChainRulesCore]]
+deps = ["Compat", "LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "de4f08843c332d355852721adb1592bce7924da3"
+uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+version = "0.9.29"
+
+[[Compat]]
+deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
+git-tree-sha1 = "919c7f3151e79ff196add81d7f4e45d91bbf420b"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "3.25.0"
+
+[[CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+
+[[DataAPI]]
+git-tree-sha1 = "dfb3b7e89e395be1e25c2ad6d7690dc29cc53b1d"
+uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
+version = "1.6.0"
+
+[[DataStructures]]
+deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
+git-tree-sha1 = "4437b64df1e0adccc3e5d1adbc3ac741095e4677"
+uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+version = "0.18.9"
+
 [[Dates]]
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
+[[DelimitedFiles]]
+deps = ["Mmap"]
+uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+
 [[Distributed]]
 deps = ["Random", "Serialization", "Sockets"]
 uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 
+[[Distributions]]
+deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"]
+git-tree-sha1 = "0fc424e725eaec6ea3e9fa8df773bee18a1ab503"
+uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
+version = "0.24.14"
+
+[[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"
 
+[[FFTW]]
+deps = ["AbstractFFTs", "FFTW_jll", "IntelOpenMP_jll", "Libdl", "LinearAlgebra", "MKL_jll", "Reexport"]
+git-tree-sha1 = "1b48dbde42f307e48685fa9213d8b9f8c0d87594"
+uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+version = "1.3.2"
+
+[[FFTW_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "5a0d4b6a22a34d17d53543bd124f4b08ed78e8b0"
+uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
+version = "3.3.9+7"
+
+[[FillArrays]]
+deps = ["LinearAlgebra", "Random", "SparseArrays"]
+git-tree-sha1 = "4705cc4e212c3c978c60b1b18118ec49b4d731fd"
+uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
+version = "0.11.5"
+
 [[Hwloc]]
 deps = ["Hwloc_jll"]
 git-tree-sha1 = "2e3d1d4ab0e7296354539b2be081f71f4b694c0b"
@@ -44,15 +123,37 @@ git-tree-sha1 = "28e837ff3e7a6c3cdb252ce49fb412c8eb3caeef"
 uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 version = "0.1.0"
 
+[[IntelOpenMP_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "d979e54b71da82f3a65b62553da4fc3d18c9004c"
+uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
+version = "2018.0.3+2"
+
 [[InteractiveUtils]]
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
+[[Interpolations]]
+deps = ["AxisAlgorithms", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"]
+git-tree-sha1 = "eb1dd6d5b2275faaaa18533e0fc5f9171cec25fa"
+uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+version = "0.13.1"
+
 [[JLLWrappers]]
 git-tree-sha1 = "a431f5f2ca3f4feef3bd7a5e94b8b8d4f2f647a0"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
 version = "1.2.0"
 
+[[KernelDensity]]
+deps = ["Distributions", "DocStringExtensions", "FFTW", "Interpolations", "StatsBase"]
+git-tree-sha1 = "09aeec87bdc9c1fa70d0b508dfa94a21acd280d9"
+uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
+version = "0.6.2"
+
+[[LazyArtifacts]]
+deps = ["Artifacts", "Pkg"]
+uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
+
 [[LibCURL]]
 deps = ["LibCURL_jll", "MozillaCACerts_jll"]
 uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
@@ -79,6 +180,12 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 [[Logging]]
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
+[[MKL_jll]]
+deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
+git-tree-sha1 = "c253236b0ed414624b083e6b72bfe891fbd2c7af"
+uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
+version = "2021.1.1+1"
+
 [[Markdown]]
 deps = ["Base64"]
 uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
@@ -87,12 +194,44 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 deps = ["Artifacts", "Libdl"]
 uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
 
+[[Missings]]
+deps = ["DataAPI"]
+git-tree-sha1 = "f8c673ccc215eb50fcadb285f522420e29e69e1c"
+uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
+version = "0.4.5"
+
+[[Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
 [[MozillaCACerts_jll]]
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
 
 [[NetworkOptions]]
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
 
+[[OffsetArrays]]
+deps = ["Adapt"]
+git-tree-sha1 = "f64fbf703e7f5af5d82ea7b7385b443bb7765ce9"
+uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+version = "1.6.1"
+
+[[OpenSpecFun_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "9db77584158d0ab52307f8c04f8e7c08ca76b5b3"
+uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
+version = "0.5.3+4"
+
+[[OrderedCollections]]
+git-tree-sha1 = "4fa2ba51070ec13fcc7517db714445b4ab986bdf"
+uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+version = "1.4.0"
+
+[[PDMats]]
+deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse", "Test"]
+git-tree-sha1 = "95a4038d1011dfdbde7cecd2ad0ac411e53ab1bc"
+uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
+version = "0.10.1"
+
 [[Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
@@ -101,6 +240,12 @@ uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 deps = ["Unicode"]
 uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
+[[QuadGK]]
+deps = ["DataStructures", "LinearAlgebra"]
+git-tree-sha1 = "12fbe86da16df6679be7521dfb39fbc861e1dc7b"
+uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+version = "2.4.1"
+
 [[REPL]]
 deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
 uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
@@ -109,25 +254,89 @@ uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 deps = ["Serialization"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
+[[Ratios]]
+git-tree-sha1 = "37d210f612d70f3f7d57d488cb3b6eff56ad4e41"
+uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
+version = "0.4.0"
+
+[[Reexport]]
+git-tree-sha1 = "57d8440b0c7d98fc4f889e478e80f268d534c9d5"
+uuid = "189a3867-3050-52da-a836-e630ba90ab69"
+version = "1.0.0"
+
 [[Requires]]
 deps = ["UUIDs"]
 git-tree-sha1 = "cfbac6c1ed70c002ec6361e7fd334f02820d6419"
 uuid = "ae029012-a4dd-5104-9daa-d747884805df"
 version = "1.1.2"
 
+[[Rmath]]
+deps = ["Random", "Rmath_jll"]
+git-tree-sha1 = "86c5647b565873641538d8f812c04e4c9dbeb370"
+uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
+version = "0.6.1"
+
+[[Rmath_jll]]
+deps = ["Libdl", "Pkg"]
+git-tree-sha1 = "d76185aa1f421306dec73c057aa384bad74188f0"
+uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
+version = "0.2.2+1"
+
 [[SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 
 [[Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
+[[SharedArrays]]
+deps = ["Distributed", "Mmap", "Random", "Serialization"]
+uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+
 [[Sockets]]
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
+[[SortingAlgorithms]]
+deps = ["DataStructures", "Random", "Test"]
+git-tree-sha1 = "03f5898c9959f8115e30bc7226ada7d0df554ddd"
+uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
+version = "0.3.1"
+
 [[SparseArrays]]
 deps = ["LinearAlgebra", "Random"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
+[[SpecialFunctions]]
+deps = ["ChainRulesCore", "OpenSpecFun_jll"]
+git-tree-sha1 = "5919936c0e92cff40e57d0ddf0ceb667d42e5902"
+uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
+version = "1.3.0"
+
+[[StaticArrays]]
+deps = ["LinearAlgebra", "Random", "Statistics"]
+git-tree-sha1 = "9da72ed50e94dbff92036da395275ed114e04d49"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "1.0.1"
+
+[[Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[StatsBase]]
+deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics"]
+git-tree-sha1 = "400aa43f7de43aeccc5b2e39a76a79d262202b76"
+uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+version = "0.33.3"
+
+[[StatsFuns]]
+deps = ["Rmath", "SpecialFunctions"]
+git-tree-sha1 = "3b9f665c70712af3264b61c27a7e1d62055dafd1"
+uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
+version = "0.9.6"
+
+[[SuiteSparse]]
+deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
+uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
+
 [[TOML]]
 deps = ["Dates"]
 uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
@@ -136,6 +345,10 @@ uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 deps = ["ArgTools", "SHA"]
 uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
 
+[[Test]]
+deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
 [[UUIDs]]
 deps = ["Random", "SHA"]
 uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
@@ -160,6 +373,12 @@ git-tree-sha1 = "59c95a188efd11c6ed762154397ed5ea94a95e30"
 uuid = "33b4df10-0173-11e9-2a0c-851a7edac40e"
 version = "0.2.7"
 
+[[WoodburyMatrices]]
+deps = ["LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "59e2ad8fd1591ea019a5259bd012d7aee15f995c"
+uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6"
+version = "0.5.3"
+
 [[Zlib_jll]]
 deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
diff --git a/Project.toml b/Project.toml
index 835d8e2f14612183dae65db2052406b4b78232bb..9befd605486c79001b4c04461c462384f2f324c3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,2 +1,3 @@
 [deps]
+KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
 VectorizedRNG = "33b4df10-0173-11e9-2a0c-851a7edac40e"
diff --git a/ZeroWeightedDistributions/test/runtests.jl b/ZeroWeightedDistributions/test/runtests.jl
index 20fe1feeba58dc55d1ef02922d0594c1e3aaaa6c..3a2ef2aee5c205dcb6144a505d4af35d5b1a0231 100644
--- a/ZeroWeightedDistributions/test/runtests.jl
+++ b/ZeroWeightedDistributions/test/runtests.jl
@@ -11,19 +11,19 @@ using Plots
 # const RNG = Xoroshiro128Star(1)
 # const tol = 1e-1
 function test_dist(d)
-    samples_vec = rand(RNG,d, 100)
-    display(samples_vec)
-    samples_map = map(x->rand(RNG,d), 1:10_000)
-    epdf = kde(samples_vec)
-    vec_err = all([abs(pdf(epdf,k) - pdf(d,k)) for k in 0.0:0.1:10] .< tol)
-    epdf = kde(samples_map)
-    map_err= all([abs(pdf(epdf,k) - pdf(d,k)) for k in 0.0:0.1:10] .< tol)
-    # displlay((vec_err,map_err))
+    # samples_vec = rand(RNG,d, 100)
+    # display(samples_vec)
+    # samples_map = map(x->rand(RNG,d), 1:10_000)
+    # epdf = kde(samples_vec)
+    # vec_err = all([abs(pdf(epdf,k) - pdf(d,k)) for k in 0.0:0.1:10] .< tol)
+    # epdf = kde(samples_map)
+    # map_err= all([abs(pdf(epdf,k) - pdf(d,k)) for k in 0.0:0.1:10] .< tol)
+    # # displlay((vec_err,map_err))
     
-    p = scatter([pdf(epdf,k) for k in 0.0:0.1:10])
-    scatter!(p,[pdf(d,k) for k in 0.0:0.1:10])
-    display(p)
-    return vec_err && map_err
+    # p = scatter([pdf(epdf,k) for k in 0.0:0.1:10])
+    # scatter!(p,[pdf(d,k) for k in 0.0:0.1:10])
+    # display(p)
+    # return vec_err && map_err
 end