Likelihood Free Markhov Chain Monte Carlo (LFMCMC)

Andrew Pulsipher

2024-11-20

Introduction

The purpose of the “lfmcmc” function is to perform a Likelihood-Free Markhov Chain Monte Carlo simulation.

Example: Using LFMCMC to calibrate an SIR Model

Setup and Running the Model

Create an SIR Model and add a small world population. Then, run the model.

library(epiworldR)

model_seed <- 122

model_sir <- ModelSIR(
  name = "COVID-19",
  prevalence = .1,
  transmission_rate = .3,
  recovery_rate = .3
)

agents_smallworld(
  model_sir,
  n = 1000,
  k = 5,
  d = FALSE,
  p = 0.01
)

verbose_off(model_sir)

run(
  model_sir,
  ndays = 50,
  seed = model_seed
)

summary(model_sir)
#> ________________________________________________________________________________
#> ________________________________________________________________________________
#> SIMULATION STUDY
#> 
#> Name of the model   : Susceptible-Infected-Recovered (SIR)
#> Population size     : 1000
#> Agents' data        : (none)
#> Number of entities  : 0
#> Days (duration)     : 50 (of 50)
#> Number of viruses   : 1
#> Last run elapsed t  : 737.00µs
#> Last run speed      : 67.84 million agents x day / second
#> Rewiring            : off
#> 
#> Global events:
#>  (none)
#> 
#> Virus(es):
#>  - COVID-19
#> 
#> Tool(s):
#>  (none)
#> 
#> Model parameters:
#>  - Recovery rate     : 0.3000
#>  - Transmission rate : 0.3000
#> 
#> Distribution of the population at time 50:
#>   - (0) Susceptible :  900 -> 285
#>   - (1) Infected    :  100 -> 0
#>   - (2) Recovered   :    0 -> 715
#> 
#> Transition Probabilities:
#>  - Susceptible  0.98  0.02  0.00
#>  - Infected     0.00  0.71  0.29
#>  - Recovered    0.00  0.00  1.00

Setup LFMCMC

# Extract the observed data from the model
obs_data <- get_today_total(model_sir)

# Define the LFMCMC simulation function
simfun <- function(params) {

  set_param(model_sir, "Recovery rate", params[1])
  set_param(model_sir, "Transmission rate", params[2])

  run(
    model_sir,
    ndays = 50
  )

  get_today_total(model_sir)

}

# Define the LFMCMC summary function
sumfun <- function(dat) {
  return(dat)
}

# Define the LFMCMC proposal function
# - Based on proposal_fun_normal from lfmcmc-meat.hpp
propfun <- function(params_prev) {
  res <- plogis(qlogis(params_prev) + rnorm(length(params_prev)))
  return(res)
}

# Define the LFMCMC kernel function
# - Based on kernel_fun_uniform from lfmcmc-meat.hpp
kernelfun <- function(stats_now, stats_obs, epsilon) {
  dnorm(sqrt(sum((stats_now - stats_obs)^2)))
}

# Create the LFMCMC model
lfmcmc_model <- LFMCMC(model_sir) |>
  set_simulation_fun(simfun) |>
  set_summary_fun(sumfun) |>
  set_proposal_fun(propfun) |>
  set_kernel_fun(kernelfun) |>
  set_observed_data(obs_data)

Run LFMCMC simulation

# Set initial parameters
par0 <- c(0.1, 0.5)
n_samp <- 2000
epsil <- 1.0

# Run the LFMCMC simulation
run_lfmcmc(
  lfmcmc = lfmcmc_model,
  params_init_ = par0,
  n_samples_ = n_samp,
  epsilon_ = epsil,
  seed = model_seed
)

# Print the results
set_stats_names(lfmcmc_model, get_states(model_sir))
set_par_names(lfmcmc_model, c("Immune recovery", "Infectiousness"))

print(lfmcmc_model)
#> ___________________________________________
#> 
#> LIKELIHOOD-FREE MARKOV CHAIN MONTE CARLO
#> 
#> N Samples : 2000
#> Elapsed t : 1.00s
#> 
#> Parameters:
#>   -Immune recovery :  0.31 [ 0.19,  0.40] (initial :  0.10)
#>   -Infectiousness  :  0.28 [ 0.19,  0.33] (initial :  0.50)
#> 
#> Statistics:
#>   -Susceptible :  284.71 [ 272.00,  307.00] (Observed:  285.00)
#>   -Infected    :    0.85 [   0.00,    0.00] (Observed:    0.00)
#>   -Recovered   :  713.94 [ 693.00,  728.00] (Observed:  715.00)
#> ___________________________________________