Lfmcmc rand engine (#51)

* Removing vscode cached

* Updating news file and correcting how seeds are set

* Adding another set of tests

* Modified vignette with better example

* Remove unneeded code in kernelfun in lfmcmc vignette

* Move test from test-lfmcmc-other-models.R to test-lfmcmc.R

* Add Andrew Pulsipher as co-author

* Fix LFMCMC unit tests

* Fix errors in examples in LFMCMC.R

* Document data type coercion in lfmcmc.R

* Add test of type coercion to test-lfmcmc.R

* Fix type coercion tests in test-lfmcmc.R to use correct values

---------

Co-authored-by: Andrew Pulsipher <[email protected]>

DESCRIPTION

@@ -7,6 +7,8 @@ Authors@R: c(
   email="[email protected]", comment = c(ORCID = "0000-0002-3171-0844")),
   person(given="Derek", family="Meyer", role=c("aut"),
   email="[email protected]", comment = c(ORCID = "0009-0005-1350-6988")),
+  person(given="Andrew", family="Pulsipher", role=c("aut"),
+  email="[email protected]", comment = c(ORCID = "0000-0002-0773-3210")),
   person(given="Susan", family="Holmes", role = "rev", comment =
          c(what = "JOSS reviewer", ORCID="0000-0002-2208-8168")),
   person(given="Abinash", family="Satapathy", role = "rev", comment =

NEWS.md

@@ -1,3 +1,21 @@
+# epiworldR 0.4-3 (development version)
+
+## New features
+
+* The package now includes the `LFMCMC` module that implements
+  the likelihood-free Markov Chain Monte Carlo algorithm. This
+  module is used to estimate the parameters of the models.
+
+* The new function `add_param()` allows the user to add parameters
+  to the model.
+
+* The new function `rm_globalevent()` allows the user to remove
+  global events from the model.
+
+* The function `today()` returns the current day (step) of the
+  simulation.
+
+
 # epiworldR 0.3-2
 
 * Starting version 0.3-0, `epiworldR` is versioned using the same version as the C++ library, `epiworld`.

R/LFMCMC.R

@@ -1,10 +1,12 @@
 #' Likelihood-Free Markhov Chain Monte Carlo (LFMCMC)
 #'
-#'
 #' @aliases epiworld_lfmcmc
+#' @param model A model of class [epiworld_model] or `NULL` (see details).
 #' @details
-#' Performs a Likelihood-Free Markhov Chain Monte Carlo simulation
-#' @param model A model of class [epiworld_model]
+#' Performs a Likelihood-Free Markhov Chain Monte Carlo simulation. When
+#' `model` is not `NULL`, the model uses the same random-number generator
+#' engine as the model. Otherwise, when `model` is `NULL`, a new random-number
+#' generator engine is created.
 #' @returns
 #' The `LFMCMC` function returns a model of class [epiworld_lfmcmc].
 #' @examples
@@ -24,14 +26,14 @@
 #'
 #' ## Setup LFMCMC
 #' # Extract the observed data from the model
-#' obs_data <- unname(as.integer(get_today_total(model_sir)))
+#' obs_data <- get_today_total(model_sir)
 #'
 #' # Define the 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)
-#'   res <- unname(as.integer(get_today_total(model_sir)))
+#'   res <- get_today_total(model_sir)
 #'   return(res)
 #' }
 #'
@@ -50,9 +52,9 @@
 #'
 #' ## Run LFMCMC simulation
 #' # Set initial parameters
-#' par0 <- as.double(c(0.1, 0.5))
+#' par0 <- c(0.1, 0.5)
 #' n_samp <- 2000
-#' epsil <- as.double(1.0)
+#' epsil <- 1.0
 #'
 #' # Run the LFMCMC simulation
 #' run_lfmcmc(
@@ -73,43 +75,74 @@
 #' get_params_mean(lfmcmc_model)
 #'
 #' @export
-LFMCMC <- function(model) {
-  if (!inherits(model, "epiworld_model"))
-    stop("model should be of class 'epiworld_model'. It is of class ", class(model))
+LFMCMC <- function(model = NULL) {
+
+  if ((length(model) > 0) && !inherits(model, "epiworld_model"))
+    stop(
+      "model should be of class 'epiworld_model'. It is of class ",
+      paste(class(model), collapse = "\", ")
+    )
 
   structure(
     LFMCMC_cpp(model),
     class = c("epiworld_lfmcmc")
   )
+
 }
 
 #' @rdname LFMCMC
 #' @param lfmcmc LFMCMC model
-#' @param params_init_ Initial model parameters
-#' @param n_samples_ Number of samples
-#' @param epsilon_ Epsilon parameter
+#' @param params_init_ Initial model parameters, treated as double
+#' @param n_samples_ Number of samples, treated as integer
+#' @param epsilon_ Epsilon parameter, treated as double
 #' @param seed Random engine seed
 #' @returns The simulated model of class [epiworld_lfmcmc].
 #' @export
-run_lfmcmc <- function(lfmcmc, params_init_, n_samples_, epsilon_, seed = NULL) UseMethod("run_lfmcmc")
+run_lfmcmc <- function(
+    lfmcmc, params_init_, n_samples_, epsilon_,
+    seed = NULL
+    ) {
+  UseMethod("run_lfmcmc")
+}
 
 #' @export
-run_lfmcmc.epiworld_lfmcmc <- function(lfmcmc, params_init_, n_samples_, epsilon_, seed = NULL) {
-  if (length(seed)) set.seed(seed)
-  run_lfmcmc_cpp(lfmcmc, params_init_, n_samples_, epsilon_, sample.int(1e4, 1))
+run_lfmcmc.epiworld_lfmcmc <- function(
+    lfmcmc, params_init_, n_samples_, epsilon_,
+    seed = NULL
+    ) {
+
+  if (length(seed))
+    set.seed(seed)
+
+  run_lfmcmc_cpp(
+    lfmcmc,
+    as.double(params_init_),
+    as.integer(n_samples_),
+    as.double(epsilon_),
+    sample.int(1e4, 1)
+  )
+
   invisible(lfmcmc)
+
 }
 
 #' @rdname LFMCMC
 #' @param lfmcmc LFMCMC model
-#' @param observed_data_ Observed data
+#' @param observed_data_ Observed data, treated as double
 #' @returns The lfmcmc model with the observed data added
 #' @export
-set_observed_data <- function(lfmcmc, observed_data_) UseMethod("set_observed_data")
+set_observed_data <- function(
+    lfmcmc, observed_data_
+    ) {
+  UseMethod("set_observed_data")
+}
 
 #' @export
 set_observed_data.epiworld_lfmcmc <- function(lfmcmc, observed_data_) {
-  set_observed_data_cpp(lfmcmc, observed_data_)
+  set_observed_data_cpp(
+    lfmcmc,
+    as.double(observed_data_)
+  )
   invisible(lfmcmc)
 }
 

inst/tinytest/test-lfmcmc.R

@@ -5,13 +5,15 @@ 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)
+                      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)
 
+# Check init of LFMCMC model without epiworld model ----------------------------
+expect_silent(lfmcmc_nomodel <- LFMCMC())
+
 # Check bad init of LFMCMC model -----------------------------------------------
-expect_error(lfmcmc_bad <- LFMCMC(), 'argument "model" is missing')
 expect_error(lfmcmc_bad <- LFMCMC(c("not_a_model")), "model should be of class 'epiworld_model'")
 
 # Create LFMCMC model ----------------------------------------------------------
@@ -22,7 +24,7 @@ expect_inherits(lfmcmc_model, "epiworld_lfmcmc")
 expect_length(class(lfmcmc_model), 1)
 
 # Extract observed data from the model
-obs_data <- unname(as.integer(get_today_total(model_sir)))
+obs_data <- get_today_total(model_sir)
 
 expect_silent(set_observed_data(lfmcmc_model, obs_data))
 
@@ -31,20 +33,19 @@ simfun <- function(params) {
   set_param(model_sir, "Recovery rate", params[1])
   set_param(model_sir, "Transmission rate", params[2])
   run(model_sir, ndays = 50)
-  res <- unname(as.integer(get_today_total(model_sir)))
+  res <- get_today_total(model_sir)
   return(res)
 }
 
 sumfun <- function(dat) { return(dat) }
 
 propfun <- function(params_prev) {
-  res <- params_prev + rnorm(length(params_prev), )
+  res <- plogis(qlogis(params_prev) + rnorm(length(params_prev)))
   return(res)
 }
 
 kernelfun <- function(stats_now, stats_obs, epsilon) {
-  ans <- sum(mapply(function(v1, v2) (v1 - v2)^2, stats_obs, stats_now))
-  return(ifelse(sqrt(ans) < epsilon, 1.0, 0.0))
+  dnorm(sqrt(sum((stats_now - stats_obs)^2)))
 }
 
 # Check adding functions to LFMCMC
@@ -55,9 +56,9 @@ expect_silent(set_kernel_fun(lfmcmc_model, kernelfun))
 
 # Run LFMCMC simulation --------------------------------------------------------
 # Initial parameters
-par0 <- as.double(c(0.1, 0.5))
+par0 <- c(0.1, 0.5)
 n_samp <- 2000
-epsil <- as.double(1.0)
+epsil <- 1.0
 
 expect_silent(run_lfmcmc(
   lfmcmc = lfmcmc_model,
@@ -72,8 +73,8 @@ expect_silent(set_par_names(lfmcmc_model, c("Immune recovery", "Infectiousness")
 
 expect_stdout(print(lfmcmc_model))
 
-expect_equal(get_stats_mean(lfmcmc_model), c(4.45, 2.6135, 992.4365))
-expect_equal(get_params_mean(lfmcmc_model), c(11.58421, 18.96851), tolerance = 0.00001)
+expect_equal(get_stats_mean(lfmcmc_model), c(284.7140, 0.8485, 713.9375))
+expect_equal(get_params_mean(lfmcmc_model), c(0.3132901, 0.2782186), tolerance = 0.00001)
 
 # Check LFMCMC using factory functions -----------------------------------------
 expect_silent(use_proposal_norm_reflective(lfmcmc_model))
@@ -87,6 +88,22 @@ expect_silent(run_lfmcmc(
   seed = model_seed
 ))
 
+# Check LFMCMC type coercion of parameters and observed data -------------------
+obs_data_int <- as.integer(obs_data)
+expect_silent(set_observed_data(lfmcmc_model, obs_data_int))
+
+par0_int <- as.integer(c(1, 5))
+n_samp_double <- as.double(2000.0)
+epsil_int <- as.integer(1)
+
+expect_silent(run_lfmcmc(
+  lfmcmc = lfmcmc_model,
+  params_init_ = par0_int,
+  n_samples_ = n_samp_double,
+  epsilon_ = epsil_int,
+  seed = model_seed
+))
+
 # Check running LFMCMC with missing parameters ---------------------------------
 expect_silent(run_lfmcmc(
   lfmcmc = lfmcmc_model,
@@ -120,3 +137,60 @@ expect_error(run_lfmcmc(
 ))
 
 expect_error(run_lfmcmc(lfmcmc = lfmcmc_model))
+
+# Check running LFMCMC without epiworld model ----------------------------------
+model_seed <- 4222
+set.seed(model_seed)
+Y <- rnorm(2000, mean = -5, sd = 2.5)
+
+# Define LFMCMC functions
+simfun <- function(par) {
+  rnorm(2000, mean = par[1], sd = par[2])
+}
+
+sumfun <- function(x) {
+  c(mean(x), sd(x))
+}
+
+propfun <- function(par) {
+
+  par_new <- par + rnorm(2, sd = 0.1)
+
+  # Reflecting par2
+  if (par_new[2] < 0) {
+    par_new[2] <- par[2] - (par_new[2] - par[2])
+  }
+
+  return(par_new)
+}
+
+kernelfun <- function(stats_now, stats_obs, epsilon) {
+
+  dnorm(sqrt(sum((stats_obs - stats_now)^2)))
+
+}
+
+# Setup LFMCMC
+lfmcmc_model <- LFMCMC() |>
+  set_simulation_fun(simfun) |>
+  set_summary_fun(sumfun) |>
+  set_proposal_fun(propfun) |>
+  set_kernel_fun(kernelfun) |>
+  set_observed_data(Y)
+
+# Run LFMCMC
+x <- run_lfmcmc(
+  lfmcmc = lfmcmc_model,
+  params_init_ = c(0, 1),
+  n_samples_ = 3000,
+  epsilon_ = 1.0,
+  seed = model_seed
+)
+
+x_means <- get_params_mean(x)
+
+expect_equivalent(
+  x_means,
+  c(-5, 2.5),
+  tolerance = 0.1
+)