-
Notifications
You must be signed in to change notification settings - Fork 3
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
1 parent
18f61dd
commit 60982d7
Showing
1 changed file
with
32 additions
and
60 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,96 +1,68 @@ | ||
# Create Model to use in LFMCMC simulation | ||
model_seed <- 122 | ||
|
||
model_sir <- ModelSIR( | ||
name = "COVID-19", | ||
prevalence = .1, | ||
transmission_rate = .9, | ||
recovery_rate = .3 | ||
) | ||
# Test just this file: tinytest::run_test_file("inst/tinytest/test-lfmcmc.R") | ||
|
||
agents_smallworld( | ||
model_sir, | ||
n = 1000, | ||
k = 5, | ||
d = FALSE, | ||
p = 0.01 | ||
) | ||
# Set model parameters | ||
model_seed <- 122 | ||
|
||
# Create and run SIR Model for LFMCMC simulation ------------------------------- | ||
model_sir <- ModelSIR(name = "COVID-19", prevalence = .1, | ||
transmission_rate = .9, 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) | ||
|
||
run( | ||
model_sir, | ||
ndays = 50, | ||
seed = model_seed | ||
) | ||
# Create LFMCMC model ---------------------------------------------------------- | ||
lfmcmc_model <- LFMCMC(model_sir) | ||
|
||
# Check initialization | ||
expect_inherits(lfmcmc_model, "epiworld_lfmcmc") | ||
|
||
# Setup LFMCMC | ||
## Extract the observed data from the model | ||
# Extract observed data from the model | ||
obs_data <- unname(as.integer(get_today_total(model_sir))) | ||
|
||
## Define the LFMCMC simulation function | ||
simfun <- function(params) { | ||
expect_silent(set_observed_data(lfmcmc_model, obs_data)) | ||
|
||
# Define LFMCMC functions | ||
simfun <- function(params) { | ||
set_param(model_sir, "Recovery rate", params[1]) | ||
set_param(model_sir, "Transmission rate", params[2]) | ||
|
||
run( | ||
model_sir, | ||
ndays = 50 | ||
) | ||
|
||
run(model_sir, ndays = 50) | ||
res <- unname(as.integer(get_today_total(model_sir))) | ||
return(res) | ||
} | ||
|
||
## Define the LFMCMC summary function | ||
sumfun <- function(dat) { | ||
return(dat) | ||
} | ||
sumfun <- function(dat) { return(dat) } | ||
|
||
## Define the LFMCMC proposal function | ||
## - Based on proposal_fun_normal from lfmcmc-meat.hpp | ||
propfun <- function(params_prev) { | ||
res <- 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) { | ||
|
||
ans <- sum(mapply(function(v1, v2) (v1 - v2)^2, | ||
stats_obs, | ||
stats_now)) | ||
|
||
ans <- sum(mapply(function(v1, v2) (v1 - v2)^2, stats_obs, stats_now)) | ||
return(ifelse(sqrt(ans) < epsilon, 1.0, 0.0)) | ||
} | ||
|
||
## 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) | ||
# Check adding functions to LFMCMC | ||
expect_silent(set_simulation_fun(lfmcmc_model, simfun)) | ||
expect_silent(set_summary_fun(lfmcmc_model, sumfun)) | ||
expect_silent(set_proposal_fun(lfmcmc_model, propfun)) | ||
expect_silent(set_kernel_fun(lfmcmc_model, kernelfun)) | ||
|
||
# Run LFMCMC simulation | ||
## Set initial parameters | ||
# Create LFMCMC simulation ----------------------------------------------------- | ||
# Initial parameters | ||
par0 <- as.double(c(0.1, 0.5)) | ||
n_samp <- 2000 | ||
epsil <- as.double(1.0) | ||
|
||
## Run the LFMCMC simulation | ||
run_lfmcmc( | ||
expect_silent(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")) | ||
expect_silent(set_stats_names(lfmcmc_model, get_states(model_sir))) | ||
expect_silent(set_par_names(lfmcmc_model, c("Immune recovery", "Infectiousness"))) | ||
|
||
print(lfmcmc_model) | ||
expect_stdout(print(lfmcmc_model)) |