new intervals model sims done, behavioural fitting still not working..

CovidAlertVaccinationModel/abm_fit_parameters.jl

@@ -22,7 +22,7 @@ const parameters = (
     β = 10.0,
     notification_parameter = 0.001,
     vaccinator_prob = 0.2,
-    app_user_fraction = 0.2,
+    app_user_fraction = 0.0,
     notification_threshold = 17,
     immunizing = true,
     immunization_delay = 14,

CovidAlertVaccinationModel/abm_timeseries.jl

@@ -9,7 +9,7 @@ const parameters = (
     base_transmission_probability = 0.001,
     recovery_rate = 1/7,
     immunization_loss_prob = 0.0055, #mean time of 6 months
-    π_base = -1.0,
+    π_base = -0.8,
     η = 0.0,
     κ = 0.0,
     ω = 0.00005,
@@ -19,7 +19,7 @@ const parameters = (
     ω_en = 0.00,
     ρ_en =  [0.0,0.0,0.0],
     γ = 0.0,
-    β = 5.0,
+    β = 10.0,
     notification_parameter = 0.001,
     vaccinator_prob = 0.2,
     app_user_fraction = 0.2,

CovidAlertVaccinationModel/src/ABM/parameter_optimization.jl

@@ -48,55 +48,53 @@ end
 
 function fit_behavioural_parameters(p_tuple, target_growth_rate)
     samples = 1
-    p_names = (:π_base,:κ, :ρ_y,:ρ_m,:ρ_o)
-    priors = Factored(Uniform(-2.0,0.0),Uniform(0.0,1.0),Uniform(0.0,1.0),Uniform(0.0,1.0),Uniform(0.0,1.0))
+    p_names = (:π_base,:ρ_y,:ρ_m,:ρ_o)
+    priors = Factored(Uniform(-2.0,0.0),Uniform(0.0,1.0),Uniform(0.0,1.0),Uniform(0.0,1.0))
     #simulation begins in august
     #30 days for opinion dynamics to stabilize, then immunization begins in september,
     #infection is introduced at the end of november 
-    sim_length = 90
+    sim_length = 120
     p_tuple_adjust = merge(p_tuple,
         (
             sim_length = sim_length,
             I_0_fraction = 0.000,
-            immunization_begin_day =30, 
+            immunization_begin_day =60, 
             infection_introduction_day = 91,
             immunizing = true,
         )
     )
     target_cumulative_vac_proportion = 0.33
-    vaccination_data = @SVector [0.0,0.043,0.115,0.03,0.005] #by month starting in august
+    vaccination_data = @SVector [0.0,0.043,0.385,0.424,0.115,0.03,0.005] #by month starting in august
     ymo_vac = @SVector [0.255,0.278,0.602]
     function cost(p)
         new_params = merge(p_tuple_adjust, NamedTuple{p_names}(ntuple(i -> p[i],length(p_names))))
         out = DebugRecorder(0,sim_length)
         model = abm(new_params,out)
-        vaccination_ts = mean.(out.recorded_status_totals.V)
-        # vaccination_ts = mean.(out.total_vaccinators)
-        # display(out.new_ymo_immunization)
-        ymo_vaccination_ts = out.new_ymo_immunization
-        total_ymo_vaccination = sum.(eachrow(ymo_vaccination_ts))
-        display(total_ymo_vaccination)
+
         target_cumulative_vaccinations = target_cumulative_vac_proportion*model.nodes
-        total_err = 0
-        target_ymo_vac = ymo_vac .* target_montly_vac .* target_cumulative_vaccinations
-        return sum((monthly_ymo_vac .- target_monthly_ymo_vac).^2)    
+        target_ymo_vac = ymo_vac .* sum(vaccination_data[1:4]) .* target_cumulative_vaccinations
+        ymo_vaccination_ts = out.new_ymo_immunization
+        total_preinfection_vaccination = sum.(eachrow(ymo_vaccination_ts))
+        # display(target_ymo_vac)
+        # display(total_preinfection_vaccination)
+        return sum((total_preinfection_vaccination .- target_ymo_vac).^2)    
     end
 #smc(priors,cost; verbose = true, nparticles = 100, parallel = true)#
-    out = ABCDE(priors,cost,1e6; verbose=true, nparticles=100,generations=50,  parallel = true) #this one has better NaN handling
+    out = ABCDE(priors,cost,1e6; verbose=true, nparticles=300,generations=100,  parallel = true) #this one has better NaN handling
     return out
 end
 
 
 function plot_behavioural_fit(particles,p_tuple)
-    p_names = (:π_base, :κ, :ρ_y,:ρ_m,:ρ_o)
+    p_names = (:π_base, :ρ_y,:ρ_m,:ρ_o)
     sim_length = 210
     samples = 1
     p_tuple_adjust = merge(p_tuple,
-        (
+    (
             sim_length = sim_length,
-            I_0_fraction = 0.005,
-            immunization_begin_day =30, 
-            infection_introduction_day = 90,
+            I_0_fraction = 0.000,
+            immunization_begin_day =60, 
+            infection_introduction_day = 91,
             immunizing = true,
         )
     )
@@ -104,7 +102,12 @@ function plot_behavioural_fit(particles,p_tuple)
     display(p)
     new_params = merge(p_tuple_adjust, NamedTuple{p_names}(ntuple(i -> p[i],length(p_names))))
     out = mean_solve(samples, new_params ,DebugRecorder)
-    display(new_params)
+    target_cumulative_vac_proportion = 0.33
+    vaccination_data = @SVector [0.0,0.043,0.385,0.424,0.115,0.03,0.005] #by month starting in august
+    ymo_vac = @SVector [0.255,0.278,0.602]
+    ymo_vaccination_ts = mean.(out.new_ymo_immunization)
+    total_preinfection_vaccination = sum.(eachrow(ymo_vaccination_ts))
+    display(total_preinfection_vaccination)
     p = plot_model(nothing,[nothing],[out],new_params.infection_introduction_day,new_params.immunization_begin_day)
     savefig(p,"behaviour_fit.pdf")
     return out

CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl

@@ -56,11 +56,11 @@ include("ABM/output.jl")
 include("ABM/solve.jl")
 include("ABM/abm.jl")
 
-
 include("IntervalsModel/intervals_model.jl")
 include("IntervalsModel/interval_overlap_sampling.jl")
 include("IntervalsModel/hh_durations_model.jl")
 include("IntervalsModel/ws_durations_model.jl")
 include("IntervalsModel/rest_durations_model.jl")
 include("IntervalsModel/plotting_functions.jl")
+
 end

CovidAlertVaccinationModel/src/IntervalsModel/plotting_functions.jl

@@ -48,7 +48,7 @@ function plot_dists(fname,dist_constructors,data)
     p = plot(p_matrix..., size = (800,600))
     savefig(p,joinpath(PACKAGE_FOLDER,"plots", "$fname.pdf"))
 end
-
+using Statistics
 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]))
 
@@ -59,7 +59,7 @@ function plot_posteriors(fname,data)
     ymo = ["Y","M","O"]
     for  i in YOUNG:OLD, j in YOUNG:OLD
         # display(data.P[i].particles)
-        hist = fit(Histogram,sym_data[i,j].particles; nbins = 40)
+        hist = StatsBase.fit(Histogram,sym_data[i,j].particles; nbins = 40)
         plot!(p_matrix[i,j],hist;legend = false, xlabel = L"\lambda_{%$i,%$j}", color = color_palette[1]  )
         # display(kernel_data)
         # vline!(p_list[i],[argmax(kernel_data)]; seriescolor = color_palette[2]