Skip to contents

Bayesian hyperparameter optimisation

hyperparameter_opt.Rmd

Overview

This vignette demonstrates the optimisation of algorithm hyperparameters using Bayesian optimisation

Load the required libraries:

Set up the simulation and run it with our initial parameters:

Initialise the cluster for parallel computing:

cores <- 2
cl <- makeCluster(cores)

clusterSetRNGStream(cl, iseed = 20250501)
Epi_pars <- data.frame(
  Pathogen = c("COVID-19", "Ebola"),
  R0 = c(2.2, 2.5),
  gen_time = c(6.5, 15.0),
  gen_time_var = c(2.1, 2.1),
  CFR = c(0.0132, 0.5),
  mortality_mean = c(10.0, 10.0),
  mortality_var = c(1.1, 1.1)
)

Noise_pars <- data.frame(
  repd_mean = 10.5,  # Reporting delay mean
  del_disp = 5.0,    # Reporting delay variance
  ur_mean = 0.3,     # Under-reporting mean
  ur_beta_a = 50.0   # Beta distribution alpha for under-reporting
)

Action_space <- data.frame(
  NPI = c("No restrictions", "Lockdown"),
  R_coeff = c(1.0, 0.3),
  R_beta_a = c(0.0, 5.0),
  cost_of_NPI = c(0.0, 0.15)
)

C_target <- 5000
D_target <- 12
r_trans_len <- 7

sim_settings <- list(
  ndays = 40 * 7, #simulation length
  start_day = 1,
  N = 1e7, # population size
  I0 = 10, # initial infections
  C_target = C_target, #target cases
  C_target_pen = C_target*1.5, #overshoot penalty threshold
  R_target = 1.0,
  D_target = D_target, #one way to get peaks at 400 is to increase this to 15
  D_target_pen = 50, #max death
  alpha = 1.3/C_target, #~proportional gain (regulates error in cases) covid
  #alpha = 3.25/C_target #~proportional gain (regulates error in cases) ebola
  alpha_d = 0*1.3/D_target,
  ovp = 5.0, #overshoot penalty
  dovp = 0*10.0, #death overshoot penalty
  gamma = 0.95, #discounting factor
  n_ens = 20L, #MC assembly size for 4
  sim_ens = 10L, #assembly size for full simulation
  rf = 14L, #days 14
  R_est_wind = 5L, #rf-2 #window for R estimation
  susceptibles = 0,
  pred_days = 28L,
  r_trans_steep = 1.5,  # Growth rate
  r_trans_len = r_trans_len,  # Number of days for the transition
  t0 = r_trans_len  / 2, # Midpoint of the transition
  pathogen = 1,
  susceptibles = 0,
  delay = 1,
  ur = 1,
  r_dir = 2,
  LD_on = 14, #on threshold
  LD_off = 7, #off threshold
  v_max_rate = 0.8,
  vac_scale = 100,
  vac_start = 370,
  delta_scale = 40,
  delta_start = 550,
  delta_multiplier = 1.75,
  v_protection_delta = (58+85)/200,
  v_protection_alpha = 0.83
)

episettings <- list(
  sim_function = Epi_MPC_run_wd,
  reward_function = reward_fun_wd,
  R_estimator = R_epiestim,
  noise_par = Noise_pars,
  epi_par = Epi_pars,
  actions = Action_space,
  sim_settings = sim_settings,
  parallel = FALSE
)

column_names <- c("days", "sim_id", "I", "Lambda", "C", "Lambda_C", "S", "Deaths", "Re", "Rew", "Rest", "R0est", "policy", "R_coeff")
episim_data <- data.frame(matrix(0, nrow = (sim_settings$ndays), ncol = length(column_names)))
colnames(episim_data) <- column_names

episim_data$policy <- rep(1, sim_settings$ndays)
episim_data$days <- 1:(sim_settings$ndays)
episim_data[1, ] <- c(1, 1, sim_settings$I0, sim_settings$I0, Noise_pars$ur_mean * sim_settings$I0, Noise_pars$ur_mean * sim_settings$I0, sim_settings$N - sim_settings$I0, 0, Epi_pars[1, "R0"], Epi_pars[1, "R0"], 1, 1, 1, 1)

episim_data_ens <- replicate(sim_settings$sim_ens, episim_data, simplify = FALSE)

for (ii in 1:sim_settings$sim_ens) {
  episim_data_ens[[ii]]$sim_id <- rep(ii, sim_settings$ndays)
}


# Ensure actions$cost_of_NPI is numeric
Action_space$cost_of_NPI <- as.numeric(as.character(Action_space$cost_of_NPI))

# Verify numeric parameters in sim_settings
numeric_params <- c("alpha", "alpha_d", "ovp", "dovp", "C_target", "C_target_pen", "D_target", "D_target_pen", "pred_days", "sim_ens", "ndays", "N")
sim_settings[numeric_params] <- lapply(sim_settings[numeric_params], as.numeric)

# Don't forget to stop the cluster at the end

# Export required functions and objects to cluster
clusterExport(cl, ls())
clusterExport(cl, c("reward_fun_wd", "Epi_pred_wd", "Epi_MPC_run_wd", "episim_data_ens", "episettings", "Action_space", "Epi_pars", "Noise_pars", "R_estim"))

parallel::clusterEvalQ(cl, {
  library(pbapply)
  library(VGAM)
})
#> [[1]]
#>  [1] "VGAM"      "splines"   "stats4"    "pbapply"   "stats"     "graphics" 
#>  [7] "grDevices" "utils"     "datasets"  "methods"   "base"     
#> 
#> [[2]]
#>  [1] "VGAM"      "splines"   "stats4"    "pbapply"   "stats"     "graphics" 
#>  [7] "grDevices" "utils"     "datasets"  "methods"   "base"

# Add the cluster object to your settings
episettings$parallel <- TRUE
episettings$cl <- cl

# Run the epicontrol function
results <- epicontrol(episim_data_ens, episettings)
#> [1] 0
#>   |                                                          |                                                  |   0%

Defining our objective function

We want to minimise the combined costs of being off target and implementation of lockdowns. We want to optimise the projection horizon and discount factor when calculating costs at decision points.


a<- 1.3/5000

objective <- function(horizon, discount) {

  sim_settings$pred_days <- horizon
  sim_settings$gamma <- discount

  episettings <- list(
    sim_function = Epi_MPC_run_wd,
    reward_function = reward_fun_wd,
    R_estimator = R_estim,
    noise_par = Noise_pars,
    epi_par = Epi_pars,
    actions = Action_space,
    sim_settings = sim_settings,
    parallel = FALSE
  )

  results <- epicontrol(episim_data_ens, episettings)

  rewards <- numeric(length(results))
  for (i in seq_along(results)) {
    res <- results[[i]]
    cases_v <- res$C
    interventions <- res$policy

    # single bracket: keeps length matching 'interventions'
    interv_costs <- Action_space$cost_of_NPI[interventions]

    case_err <- abs(cases_v - sim_settings$C_target)

    rewards[i] <- sum(-a * case_err - interv_costs)
  }

  mean_reward <- mean(rewards, na.rm = TRUE)
  list(Score = mean_reward,
       Pred = 0)
}

Let’s do the optimisation:

OPT_Res <- BayesianOptimization(objective,
                                bounds = list(horizon = c(14L, 35L),
                                              discount = c(0.7, 0.99)),
                                init_points = 5, n_iter = 30,
                                acq = "ucb", kappa = 2.576, eps = 0.0,
                                verbose = TRUE)
#> elapsed = 25.748 Round = 1   horizon = 20.0000   discount = 0.960533 Value = -289.4258 
#> elapsed = 29.116 Round = 2   horizon = 22.0000   discount = 0.9739558    Value = -287.4896 
#> elapsed = 32.574 Round = 3   horizon = 26.0000   discount = 0.8916314    Value = -300.3452 
#> elapsed = 40.853 Round = 4   horizon = 33.0000   discount = 0.8824431    Value = -340.5870 
#> elapsed = 22.722 Round = 5   horizon = 18.0000   discount = 0.717918 Value = -309.5872 
#> elapsed = 20.264 Round = 6   horizon = 14.0000   discount = 0.9900   Value = -388.0072 
#> elapsed = 27.276 Round = 7   horizon = 21.0000   discount = 0.7019426    Value = -288.4061 
#> elapsed = 28.103 Round = 8   horizon = 21.0000   discount = 0.9900   Value = -289.2448 
#> elapsed = 36.017 Round = 9   horizon = 29.0000   discount = 0.7050185    Value = -319.2199 
#> elapsed = 29.606 Round = 10  horizon = 24.0000   discount = 0.7311516    Value = -291.2297 
#> elapsed = 40.818 Round = 11  horizon = 35.0000   discount = 0.7000   Value = -349.5131 
#> elapsed = 27.498 Round = 12  horizon = 22.0000   discount = 0.8379493    Value = -290.3592 
#> elapsed = 29.419 Round = 13  horizon = 24.0000   discount = 0.9900   Value = -288.6685 
#> elapsed = 25.714 Round = 14  horizon = 22.0000   discount = 0.7987672    Value = -286.1254 
#> elapsed = 27.365 Round = 15  horizon = 22.0000   discount = 0.8193503    Value = -288.8146 
#> elapsed = 27.287 Round = 16  horizon = 22.0000   discount = 0.9043957    Value = -284.1406 
#> elapsed = 28.713 Round = 17  horizon = 22.0000   discount = 0.9900   Value = -289.6519 
#> elapsed = 29.03  Round = 18  horizon = 22.0000   discount = 0.9877141    Value = -292.0235 
#> elapsed = 26.081 Round = 19  horizon = 22.0000   discount = 0.7002133    Value = -283.5007 
#> elapsed = 28.229 Round = 20  horizon = 22.0000   discount = 0.936678 Value = -283.4220 
#> elapsed = 27.899 Round = 21  horizon = 22.0000   discount = 0.8680492    Value = -282.9043 
#> elapsed = 26.538 Round = 22  horizon = 22.0000   discount = 0.9176258    Value = -290.0770 
#> elapsed = 26.605 Round = 23  horizon = 22.0000   discount = 0.7225638    Value = -284.3034 
#> elapsed = 26.023 Round = 24  horizon = 22.0000   discount = 0.8049441    Value = -286.8484 
#> elapsed = 26.595 Round = 25  horizon = 22.0000   discount = 0.8184396    Value = -293.2562 
#> elapsed = 25.782 Round = 26  horizon = 22.0000   discount = 0.7796007    Value = -288.4966 
#> elapsed = 26.603 Round = 27  horizon = 22.0000   discount = 0.9867218    Value = -287.4619 
#> elapsed = 26.012 Round = 28  horizon = 22.0000   discount = 0.8082374    Value = -293.9573 
#> elapsed = 26.213 Round = 29  horizon = 22.0000   discount = 0.9673567    Value = -287.8284 
#> elapsed = 26.366 Round = 30  horizon = 22.0000   discount = 0.8244174    Value = -293.3437 
#> elapsed = 25.968 Round = 31  horizon = 22.0000   discount = 0.822453 Value = -286.4732 
#> elapsed = 26.499 Round = 32  horizon = 22.0000   discount = 0.807405 Value = -291.9611 
#> elapsed = 27.056 Round = 33  horizon = 22.0000   discount = 0.8030033    Value = -288.3032 
#> elapsed = 27.057 Round = 34  horizon = 22.0000   discount = 0.9433525    Value = -287.4136 
#> elapsed = 26.847 Round = 35  horizon = 22.0000   discount = 0.8788171    Value = -286.3122 
#> 
#>  Best Parameters Found: 
#> Round = 21   horizon = 22.0000   discount = 0.8680492    Value = -282.9043