diff --git a/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl b/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl
index ca0f71b0dd5e02932a18c0f5c8447918e430dcfb..93886bdb93965756add734d02a9665e05ffa94e2 100644
--- a/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl
+++ b/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl
@@ -12,7 +12,6 @@ import StatsBase:mean
 using Dates
 using ThreadsX
 using DelimitedFiles
-
 using NamedTupleTools
 # using CUDA
 using NetworkLayout:Stress
@@ -30,7 +29,7 @@ include("graphs.jl")
 const RNG = Xoroshiro128Star()
 const population = 14.57e6 #population of ontario
 
-const color_palette = palette(:seaborn_bright) #color theme for the plots
+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()
@@ -43,18 +42,20 @@ using BenchmarkTools
 function main()
     # display(size.(agent_model.demographic_index_vectors))
 
-    agent_model = AgentModel(5000)
+    # Random.seed!(RNG,1)
+    VectorizedRNG.seed!(1)
+    agent_model = AgentModel(1000)
     u_0 = get_u_0(length(agent_model.demographics))
     steps = 300
-    Random.seed!(RNG,1)
 
     @btime solve!($u_0,$get_parameters(),$steps,$agent_model,$vaccinate_uniformly!);    
     # agent_model = AgentModel(10_000)
-    @btime solve!($u_0,$get_parameters(),$steps,$agent_model,$vaccinate_uniformly!);
+    # @btime solve!($u_0,$get_parameters(),$steps,$agent_model,$vaccinate_uniformly!);
+    println("done")
     # sol,graphs = solve!(u_0,get_parameters(),steps,agent_model,vaccinate_uniformly!);
     # return aggregate_timeseries(sol)
     # plot_model_spatial_gif(agent_model.base_network,graphs,sol1)
 end
 
 
-end
\ No newline at end of file
+end
diff --git a/CovidAlertVaccinationModel/src/graphs.jl b/CovidAlertVaccinationModel/src/graphs.jl
index e7b95cd2c39461b506ec7836d00f168ca54bfbd1..b0a52082706d7b8a35c9cd90bc3d2e0e0aa305cf 100644
--- a/CovidAlertVaccinationModel/src/graphs.jl
+++ b/CovidAlertVaccinationModel/src/graphs.jl
@@ -39,8 +39,8 @@ end
 function advance_mixing_graph!(t,mixing_graph,index_vectors)
     mixing_graph.graph .= 0
     mixing_graph.graph .= mixing_graph.base_graph
-    generate_mixing_graph!(mixing_graph.graph, index_vectors,mixing_graph.contact_vector_rest[t])
-    generate_mixing_graph!(mixing_graph.graph, index_vectors,mixing_graph.contact_vector_ws[t])
+    generate_mixing_graph!(mixing_graph.graph, index_vectors,mixing_graph.contact_vector_rest[t],contact_time_distribution_matrix)
+    generate_mixing_graph!(mixing_graph.graph, index_vectors,mixing_graph.contact_vector_ws[t],contact_time_distribution_matrix)
     return nothing
 end
 
@@ -85,7 +85,7 @@ function generate_contact_vectors!(ij_dist,ji_dist,i_to_j_contacts, j_to_i_conta
     return nothing
 end
 #add a bipartite graph derived from mixing matrices onto g
-function generate_mixing_graph!(g,index_vectors,contacts)
+function generate_mixing_graph!(g,index_vectors,contacts,edge_weight_dists)
     for i in 1:length(index_vectors), j in 1:length(index_vectors)  
         random_bipartite_graph_fast_CL!(g,index_vectors[i],index_vectors[j],contacts.contact_array[i,j],contacts.contact_array[j,i])
     end
@@ -117,8 +117,8 @@ function random_bipartite_graph_fast_CL!(g::T,anodes,bnodes,aseq,bseq) where T<:
     @assert sum(aseq) == sum(bseq) "degree sequences must have equal sum"
     astubs = sample(RNG,anodes,StatsBase.weights(aseq./m), m)
     bstubs = sample(RNG,bnodes,StatsBase.weights(bseq./m), m)
-    for k in 1:m
-        g[astubs[k],bstubs[k]] = 1
+    for i in 1:m
+        g[astubs[i],bstubs[i]] = 1
     end
     return g
 end
diff --git a/CovidAlertVaccinationModel/src/model.jl b/CovidAlertVaccinationModel/src/model.jl
index e81b7c860536f2cf512be1851ec6954e6c802a22..89945cf6b4d172c6d0d14bbc3688e87aadecb8d9 100644
--- a/CovidAlertVaccinationModel/src/model.jl
+++ b/CovidAlertVaccinationModel/src/model.jl
@@ -25,13 +25,13 @@ function agents_step!(t,u_next,u,population_list,graph,params,index_vectors,vacc
     if t >= vaccination_start_day
         vaccination_algorithm!(Int(floor(vaccines_per_day*length(u))),u_next,u,graph,population_list,index_vectors) 
     end
-    for i in 1:length(u)
+    @inbounds for i in 1:length(u)
         agent_status = u[i]
         agent_demo = population_list[i]
         if agent_status == Susceptible
-            for j in 1:length(u) #LightGraphs.neighbors(graph,i)
-                if graph[i,j] == 1 && u[j] == Infected && rand(RNG) < contact_weight(p, agent_demo,population_list[j])
-                    u_next[i] = Infected
+            @inbounds for j in 1:length(u) #LightGraphs.neighbors(graph,i)
+                if graph[i,j] == 1 && u[j] == Infected && rand(RNG) < contact_weight(p, agent_demo,population_list[j]) #only need contact weights for infected subnetwork
+                    (u_next[i] = Infected)
                 end
             end
         elseif agent_status == Infected
diff --git a/IntervalsModel/Manifest.toml b/IntervalsModel/Manifest.toml
index 9289de637ad611adcb11d0433ab0083d2626e46d..bf4fca4d9e5fbbf6d3f7a0c2ba45ee277c15499d 100644
--- a/IntervalsModel/Manifest.toml
+++ b/IntervalsModel/Manifest.toml
@@ -115,6 +115,11 @@ git-tree-sha1 = "ac5f2213e56ed8a34a3dd2f681f4df1166b34929"
 uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
 version = "0.12.6"
 
+[[Combinatorics]]
+git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860"
+uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+version = "1.0.2"
+
 [[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"
@@ -359,6 +364,12 @@ git-tree-sha1 = "1179250d910c99810d8a7ff55c50c4ed68c77a58"
 uuid = "e33a78d0-f292-5ffc-b300-72abe9b543c8"
 version = "2.4.0+0"
 
+[[HypothesisTests]]
+deps = ["Combinatorics", "Distributions", "LinearAlgebra", "Random", "Rmath", "Roots", "Statistics", "StatsBase"]
+git-tree-sha1 = "552892528991c3e17eb60e623f1ac94c0663eb7d"
+uuid = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
+version = "0.10.2"
+
 [[IfElse]]
 git-tree-sha1 = "28e837ff3e7a6c3cdb252ce49fb412c8eb3caeef"
 uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
@@ -833,6 +844,12 @@ git-tree-sha1 = "d76185aa1f421306dec73c057aa384bad74188f0"
 uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
 version = "0.2.2+1"
 
+[[Roots]]
+deps = ["Printf"]
+git-tree-sha1 = "369e25546984dff5df351bc056fccc30de615080"
+uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
+version = "1.0.8"
+
 [[SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 
diff --git a/IntervalsModel/Project.toml b/IntervalsModel/Project.toml
index a0a0b53522660c7867280c851f3835ab70cbfb72..8ba621b310fe1181f9315673017013ac2395c8aa 100644
--- a/IntervalsModel/Project.toml
+++ b/IntervalsModel/Project.toml
@@ -10,6 +10,7 @@ CovidAlertVaccinationModel = "9260c4ec-b5cf-4bc2-ad29-e1d23bf2bd6f"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
 JLSO = "9da8a3cd-07a3-59c0-a743-3fdc52c30d11"
 KissABC = "9c9dad79-530a-4643-a18b-2704674d4108"
diff --git a/IntervalsModel/plots/Distributions.Normal_hh.pdf b/IntervalsModel/plots/Distributions.Normal_hh.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..6757ba382a52fb4786cdd1b4c736cc363ba29afc
Binary files /dev/null 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
new file mode 100644
index 0000000000000000000000000000000000000000..c85459d516e568c7b2a06e7ea5bc047a2175c872
Binary files /dev/null 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
new file mode 100644
index 0000000000000000000000000000000000000000..9dbf8ec6f609f582ed8c9797a97588d5f6b32b6c
Binary files /dev/null and b/IntervalsModel/plots/Distributions.Weibull_hh.pdf 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
new file mode 100644
index 0000000000000000000000000000000000000000..10fb973badeed1480cc7ee022c8a988b3c5c5177
Binary files /dev/null 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
new file mode 100644
index 0000000000000000000000000000000000000000..03d29483aa66e1561bf528f35ec764a22057e55e
Binary files /dev/null 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
new file mode 100644
index 0000000000000000000000000000000000000000..590d3cf5ae8c04be68ad6378acd512a0dce79858
Binary files /dev/null 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
new file mode 100644
index 0000000000000000000000000000000000000000..72ef2b4b548c167659ba79fbb7e886cdfd4df021
Binary files /dev/null and b/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Poisson, T} where T_ws.pdf differ
diff --git a/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Weibull, T} where T_rest.pdf b/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Weibull, T} where T_rest.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..05844bdcf18c09491b42eacc43ff7055ec21ba6c
Binary files /dev/null and b/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Weibull, T} where T_rest.pdf differ
diff --git a/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Weibull, T} where T_ws.pdf b/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Weibull, T} where T_ws.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..a4bd200925d467679e037579a1af857fd726a2c9
Binary files /dev/null and b/IntervalsModel/plots/ZeroWeightedDistributions.ZWDist{Distributions.Weibull, T} where T_ws.pdf differ
diff --git a/IntervalsModel/plots/hh.pdf b/IntervalsModel/plots/hh.pdf
index 38123c8048b694cdb472e8481da98192dadfe564..6d1520f30443388c556edf18cb49f2a242a08138 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 d206a27805ca5948a0e22fea0391defbf1985b61..7b0b1497c81930dabd94aabd149b32c5cf270fa5 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 472c2d44417e4e27cef9621fa7687b11d7f61a1f..2390c675a8dc57aefd07b6934bce38945f125ad1 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 c08687ddd0d038d7dcfa0584408a2eb9cbb85051..e60971face4281f8eab838dd03a19238dbe68faa 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 ed9fae52b48b1c7c870aa9a0062a38ccddc8479e..3e687ae2e68205c56786e7733a02253ae1a84cc5 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 99445fa56b8a39ae77c31625ad7770a8395c722f..0f1731c22bd459912895e1ccee17c0ff42caa9ec 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 7e557e70ccdbb713564e46ebbc7e064bd98a3249..7153f85d9e34144f53a5b29b6f26a7577c9452c2 100644
--- a/IntervalsModel/src/IntervalsModel.jl
+++ b/IntervalsModel/src/IntervalsModel.jl
@@ -12,9 +12,11 @@ const PACKAGE_FOLDER = dirname(dirname(pathof(IntervalsModel)))
 using DataStructures:OrderedDict
 using Serialization
 
+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
 
 include("interval_overlap_sampling.jl")
 include("utils.jl")
@@ -29,10 +31,10 @@ using JLSO
 using Plots
 
 
-const μ_bounds = (6,12*6) 
-const σ_bounds = (1,48)
+const μ_bounds = (1,144) 
+const σ_bounds = (1,144)
 const α_bounds = (0.0,0.8) 
-
+pgfplotsx()
 function do_hh(particles)
     dists = [
         Normal,
@@ -42,7 +44,7 @@ function do_hh(particles)
     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]),
+        ([(0.1,144)  for i = 1:6], [(0.1,144.0) for i = 1:6]),
     ])
     bayesian_estimate("hh",err_hh,dists,bounds_list,particles)
 end
@@ -51,14 +53,14 @@ function do_ws(particles)
     dists = [
         ZWDist{Normal},
         ZWDist{Poisson},
-        ZWDist{Weibull}
+        # 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,12*6)  for i = 1:6], [(0.1,48.0) 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
     bayesian_estimate("ws",err_ws,dists,bounds_list,particles)
@@ -75,7 +77,7 @@ function do_rest(particles)
     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]],
+        [[α_bounds for i = 1:6],[(0.1,144)  for i = 1:6], [(0.1,144.0) for i = 1:6]],
     ])
     @show bounds_list
     bayesian_estimate("rest",err_rest,dists,bounds_list,particles)
@@ -89,7 +91,7 @@ function bayesian_estimate(fname,err_func,dists,bounds_list,particles)
 
         @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.90,  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 = 2,  parallel = true) #apply sequential monte carlo with 200 particles
         
         return dist => out
     end |> OrderedDict
@@ -100,14 +102,21 @@ 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)))
-
+    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)
+    end
 end
 
 
diff --git a/IntervalsModel/src/plotting_functions.jl b/IntervalsModel/src/plotting_functions.jl
index 6d25eb3070042389699a853c8f0f7d3b30f6f178..69c97e74390c4b37c56aabf2390e22234d340745 100644
--- a/IntervalsModel/src/plotting_functions.jl
+++ b/IntervalsModel/src/plotting_functions.jl
@@ -8,16 +8,20 @@ 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])
     ymo = ["Y","M","O"]
+    x_range = 0.0:144.0
     for i in YOUNG:OLD, j in YOUNG:OLD
-        for (dist_constructor,p_estimate) in zip(dist_constructors,p_estimate_as_arrays)
+        for (k,(dist_constructor,p_estimate)) in enumerate(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) 
-            dist_pts = [pdf.(dists,i) for i in 0.0:144.0]
+        
+            dist_pts = [pdf.(dists,i) for i in x_range]
             mean_dat = median.(dist_pts)
-            err_down = quantile.(dist_pts,0.05)
-            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))
+            display(typeof(dist_pts))
+            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])
         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
@@ -27,25 +31,28 @@ function plot_dists(fname,dist_constructors,data)
 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)
+function plot_posteriors(fname,data)
+    p_list = map(x -> plot(),1:length(data.P))
     for i in 1:length(data.P)
-        a = stephist(
-            data.P[i].particles;
-            normalize = true,
-            title = "$(parameter_names[div(i,6)])_(i % 6)"
-        )
-        push!(p_list,a)
+        # display(data.P[i].particles)
+        hist = fit(Histogram,data.P[i].particles; nbins = 50)
+        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)
+        # display(kernel_data)
+        # vline!(p_list[i],[argmax(kernel_data)]; seriescolor = color_palette[2])
+
+    end   
+    p = plot(p_list...;layout=(length(p_list) ÷ 6,6), size = (1000,(length(p_list) ÷ 6)*300), seriescolor = color_palette[1])
+    savefig(p,joinpath(PACKAGE_FOLDER,"plots","$fname.pdf"))
+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
-    p = plot(p_list...)
-    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 3b16688679bd4e17c0f5ceec10682d4596bb9d3f..1512128e33c2589391a4aa3857a0078ba2f0818e 100644
--- a/IntervalsModel/src/ws_durations_model.jl
+++ b/IntervalsModel/src/ws_durations_model.jl
@@ -16,7 +16,7 @@ const rest_data = (
     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 = 2000
+const comparison_samples = 1000
 
 ws_distributions = CovidAlertVaccinationModel.initial_workschool_mixing_matrix
 
@@ -40,12 +40,12 @@ 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)
-                # kde_est = kde(sample_list)
+                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
-                errsum += (mean(sample_list) - mean(ws_samples[age_sample]))^2
+                errsum += mapreduce(+,0:2:durmax) do i
+                    return (pdf(kde_est,i) - pdf(kerneldensity_data[age_sample],i))^2
+                end
+                # errsum += err#(mean(sample_list) - mean(ws_samples[age_sample]))^2
             end
         end
     end
diff --git a/Manifest.toml b/Manifest.toml
new file mode 100644
index 0000000000000000000000000000000000000000..976d564f40ea7f42d0f9812717f91ea87ce5d799
--- /dev/null
+++ b/Manifest.toml
@@ -0,0 +1,169 @@
+# This file is machine-generated - editing it directly is not advised
+
+[[ArgTools]]
+uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+
+[[ArrayInterface]]
+deps = ["IfElse", "LinearAlgebra", "Requires", "SparseArrays"]
+git-tree-sha1 = "ee07ae00e3cc277dcfa5507ce25be522313ecc3e"
+uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
+version = "3.1.1"
+
+[[Artifacts]]
+uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+
+[[Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[Downloads]]
+deps = ["ArgTools", "LibCURL", "NetworkOptions"]
+uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+
+[[Hwloc]]
+deps = ["Hwloc_jll"]
+git-tree-sha1 = "2e3d1d4ab0e7296354539b2be081f71f4b694c0b"
+uuid = "0e44f5e4-bd66-52a0-8798-143a42290a1d"
+version = "1.2.0"
+
+[[Hwloc_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "1179250d910c99810d8a7ff55c50c4ed68c77a58"
+uuid = "e33a78d0-f292-5ffc-b300-72abe9b543c8"
+version = "2.4.0+0"
+
+[[IfElse]]
+git-tree-sha1 = "28e837ff3e7a6c3cdb252ce49fb412c8eb3caeef"
+uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
+version = "0.1.0"
+
+[[InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[JLLWrappers]]
+git-tree-sha1 = "a431f5f2ca3f4feef3bd7a5e94b8b8d4f2f647a0"
+uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
+version = "1.2.0"
+
+[[LibCURL]]
+deps = ["LibCURL_jll", "MozillaCACerts_jll"]
+uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
+
+[[LibCURL_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
+uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
+
+[[LibGit2]]
+deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+
+[[LibSSH2_jll]]
+deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
+uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+
+[[Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[LinearAlgebra]]
+deps = ["Libdl"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[MbedTLS_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+
+[[MozillaCACerts_jll]]
+uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+
+[[NetworkOptions]]
+uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+
+[[Pkg]]
+deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
+[[Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
+[[Random]]
+deps = ["Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[Requires]]
+deps = ["UUIDs"]
+git-tree-sha1 = "cfbac6c1ed70c002ec6361e7fd334f02820d6419"
+uuid = "ae029012-a4dd-5104-9daa-d747884805df"
+version = "1.1.2"
+
+[[SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+
+[[Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[TOML]]
+deps = ["Dates"]
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+
+[[Tar]]
+deps = ["ArgTools", "SHA"]
+uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+
+[[UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+
+[[UnPack]]
+git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
+uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
+version = "1.0.2"
+
+[[Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+
+[[VectorizationBase]]
+deps = ["ArrayInterface", "Hwloc", "IfElse", "Libdl", "LinearAlgebra"]
+git-tree-sha1 = "338930400f561a120b9b317a456c7c1cd62eac13"
+uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
+version = "0.18.10"
+
+[[VectorizedRNG]]
+deps = ["Distributed", "Random", "UnPack", "VectorizationBase"]
+git-tree-sha1 = "59c95a188efd11c6ed762154397ed5ea94a95e30"
+uuid = "33b4df10-0173-11e9-2a0c-851a7edac40e"
+version = "0.2.7"
+
+[[Zlib_jll]]
+deps = ["Libdl"]
+uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+
+[[nghttp2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
diff --git a/Project.toml b/Project.toml
new file mode 100644
index 0000000000000000000000000000000000000000..835d8e2f14612183dae65db2052406b4b78232bb
--- /dev/null
+++ b/Project.toml
@@ -0,0 +1,2 @@
+[deps]
+VectorizedRNG = "33b4df10-0173-11e9-2a0c-851a7edac40e"
diff --git a/ZeroWeightedDistributions/src/ZeroWeightedDistributions.jl b/ZeroWeightedDistributions/src/ZeroWeightedDistributions.jl
index 8b5b044ecc95e5883dac1a8657fbf0f715e8cb7f..9a6c8f554cc9f0a29d8913a9f25dc4ab41248e43 100644
--- a/ZeroWeightedDistributions/src/ZeroWeightedDistributions.jl
+++ b/ZeroWeightedDistributions/src/ZeroWeightedDistributions.jl
@@ -25,7 +25,7 @@ function Distributions.pdf(d::ZWDist, x)
     if x == 0 
         return d.α + (1-d.α)*pdf(d.base_dist,0)
     else
-        return pdf(d.base_dist,0)
+        return pdf(d.base_dist,x)
     end
 end