first attempt with O distribution shift.. did not work

CovidAlertVaccinationModel/Manifest.toml

@@ -520,8 +520,8 @@ uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
 version = "0.6.3"
 
 [[KissABC]]
-deps = ["AbstractMCMC", "Distributions", "MonteCarloMeasurements", "Random"]
-git-tree-sha1 = "693d4bd4cbb9ad6515414e0935bbd646b15fc16a"
+deps = ["AbstractMCMC", "Distributed", "Distributions", "MonteCarloMeasurements", "Random", "ThreadsX"]
+path = "/home/peterj/Projects/PRs/KissABC.jl/"
 uuid = "9c9dad79-530a-4643-a18b-2704674d4108"
 version = "3.0.1"
 

CovidAlertVaccinationModel/src/ABM/parameter_optimization.jl

@@ -69,17 +69,19 @@ function fit_all_parameters(p_tuple)
     return NamedTuple{p_names}(ntuple(i -> out.P[i].particles,length(p_names)))
 end
 function fit_distribution_parameters(p_tuple)
-    p_names = (:ω,:β_y,:π_base_y,:π_base_m,:π_base_o,:α_y,:α_m,:α_o,:O_distribution_shift)
+    p_names = (:ω,:β_y,:β_m,:β_o,:π_base_y,:π_base_m,:π_base_o,:α_y,:α_m,:α_o,:O_distribution_shift)
     priors = Factored(
         Uniform(0.0,0.01),
         Uniform(0.0,0.5),  
+        Uniform(0.0,0.5),  
+        Uniform(0.0,1.0),  
         Uniform(-5.0,3.0),
         Uniform(-5.0,3.0),
         Uniform(-5.0,3.0),
         TriangularDist(23.0, 44.0, 34.0),
         TriangularDist(22.0, 54.0, 40.0),
         TriangularDist(9.0, 59.0, 39.0),   
-        Uniform(0.0,0.2),   
+        Uniform(0.0,1.0),   
     )
     #simulation begins in july
     #60 days for opinion dynamics to stabilize, then immunization begins in september,
@@ -113,7 +115,7 @@ function fit_distribution_parameters(p_tuple)
     
         return sum((normalize_by_pop(total_preinfection_vaccination .- target_ymo_vac)).^2)  
         + sum(normalize_by_pop((total_postinf_vaccination .- target_ymo_vac)).^2)  
-        + sum(normalize_by_pop((final_size .- target_final_size)).^2)
+        + 5*sum(normalize_by_pop((final_size .- target_final_size)).^2)
         # display(final_size)
         # display( sum.(eachrow(ymo_vaccination_ts[:,1:end])))
         # display( length.(model.index_vectors))
@@ -121,7 +123,7 @@ function fit_distribution_parameters(p_tuple)
     end
     # display(cost((0.0,0.0005,-3.0,-3.0,-0.2,0.25,0.25,0.0,0.1)))
     #
-    out = smc(priors,cost; verbose = true, nparticles = 300, parallel = true)
+    out = smc(priors,cost; verbose = true, nparticles = 800, parallel = true)
     return NamedTuple{p_names}(ntuple(i -> out.P[i].particles,length(p_names)))
 end
 function fit_parameters(default_parameters)
@@ -131,17 +133,17 @@ function fit_parameters(default_parameters)
     fit_dist_parameters_path = joinpath(PACKAGE_FOLDER,"abm_parameter_fits","fit_dist_parameters.dat")
 
 
-    output = fit_distribution_parameters(default_parameters)
+    # output = fit_distribution_parameters(default_parameters)
     # serialize(fit_all_parameters_path,output)
 
     
-    # fitted_parameter_tuple = deserialize(fit_all_parameters_path)
-    # fitted_sol,avg_populations = plot_fitting_posteriors("post_inf_fitting",fitted_parameter_tuple,default_parameters)
+    fitted_parameter_tuple = deserialize(fit_all_parameters_path)
+    fitted_sol,avg_populations = plot_fitting_posteriors("post_inf_fitting",fitted_parameter_tuple,default_parameters)
 
-    # final_vac = sum.(eachrow(mean.(fitted_sol.daily_immunized_by_age)))
-    # final_size = sum.(eachrow(mean.(fitted_sol.daily_unvac_cases_by_age)))
-    # display(final_vac./avg_populations)
-    # display(final_size./avg_populations)
+    final_vac = sum.(eachrow(mean.(fitted_sol.daily_immunized_by_age)))
+    final_size = sum.(eachrow(mean.(fitted_sol.daily_unvac_cases_by_age)))
+    display(final_vac./avg_populations)
+    display(final_size./avg_populations)
     return fitted_sol
 end 
 
@@ -152,13 +154,13 @@ function plot_fitting_posteriors(fname,particles_tuple,parameters)
     avg_populations = [0.0,0.0,0.0]
     names = keys(particles_tuple)
     samples = collect(zip(particles_tuple...))
-    for p_set in sample(samples,10)
+    for p_set in sample(samples,20)
         p_set_named = NamedTuple{names}(p_set)
         sol = abm(merge(parameters,p_set_named),output_recorder)
         avg_populations .+= length.(sol.index_vectors)
         fit!(stat_recorder,output_recorder)
     end
-    avg_populations ./= 10
+    avg_populations ./= 20
     p = plot_model(nothing,[nothing],[stat_recorder],parameters.infection_introduction_day,parameters.immunization_begin_day)
     savefig(p, "$fname.pdf")