Likelihood Free Markhov Chain Monte Carlo (LFMCMC)
Andrew Pulsipher
2024-11-25
Source:vignettes/likelihood-free-mcmc.Rmd
likelihood-free-mcmc.Rmd
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 : 699.00µs
#> Last run speed : 71.53 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.27 [ 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)
#> ___________________________________________