Lfmcmc rand engine

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
@@ -73,14 +75,19 @@
 #' 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
@@ -91,25 +98,51 @@ LFMCMC <- function(model) {
 #' @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
 #' @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-other-models.R

@@ -0,0 +1,53 @@
+
+set.seed(4222)
+Y <- rnorm(2000, mean = -5, sd = 2.5)
+
+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)))
+
+}
+
+
+lfmcmc_model <- LFMCMC() |>
+  set_simulation_fun(simfun) |>
+  set_summary_fun(sumfun) |>
+  set_proposal_fun(propfun) |>
+  set_kernel_fun(kernelfun) |>
+  set_observed_data(Y)
+
+x <- run_lfmcmc(
+  lfmcmc = lfmcmc_model,
+  params_init_ = c(0, 1),
+  n_samples_ = 3000,
+  epsilon_ = 1.0,
+  seed = 4222
+)
+
+x_means <- get_params_mean(x)
+
+expect_equivalent(
+  x_means,
+  c(-5, 2.5),
+  tolerance = 0.1
+)

src/cpp11.cpp

@@ -398,17 +398,17 @@ extern "C" SEXP _epiworldR_LFMCMC_cpp(SEXP model) {
   END_CPP11
 }
 // lfmcmc.cpp
-SEXP run_lfmcmc_cpp(SEXP lfmcmc, std::vector<epiworld_double> params_init_, size_t n_samples_, epiworld_double epsilon_, int seed);
+SEXP run_lfmcmc_cpp(SEXP lfmcmc, std::vector<double> params_init_, size_t n_samples_, double epsilon_, int seed);
 extern "C" SEXP _epiworldR_run_lfmcmc_cpp(SEXP lfmcmc, SEXP params_init_, SEXP n_samples_, SEXP epsilon_, SEXP seed) {
   BEGIN_CPP11
-    return cpp11::as_sexp(run_lfmcmc_cpp(cpp11::as_cpp<cpp11::decay_t<SEXP>>(lfmcmc), cpp11::as_cpp<cpp11::decay_t<std::vector<epiworld_double>>>(params_init_), cpp11::as_cpp<cpp11::decay_t<size_t>>(n_samples_), cpp11::as_cpp<cpp11::decay_t<epiworld_double>>(epsilon_), cpp11::as_cpp<cpp11::decay_t<int>>(seed)));
+    return cpp11::as_sexp(run_lfmcmc_cpp(cpp11::as_cpp<cpp11::decay_t<SEXP>>(lfmcmc), cpp11::as_cpp<cpp11::decay_t<std::vector<double>>>(params_init_), cpp11::as_cpp<cpp11::decay_t<size_t>>(n_samples_), cpp11::as_cpp<cpp11::decay_t<double>>(epsilon_), cpp11::as_cpp<cpp11::decay_t<int>>(seed)));
   END_CPP11
 }
 // lfmcmc.cpp
-SEXP set_observed_data_cpp(SEXP lfmcmc, std::vector< int > observed_data_);
+SEXP set_observed_data_cpp(SEXP lfmcmc, std::vector< double > observed_data_);
 extern "C" SEXP _epiworldR_set_observed_data_cpp(SEXP lfmcmc, SEXP observed_data_) {
   BEGIN_CPP11
-    return cpp11::as_sexp(set_observed_data_cpp(cpp11::as_cpp<cpp11::decay_t<SEXP>>(lfmcmc), cpp11::as_cpp<cpp11::decay_t<std::vector< int >>>(observed_data_)));
+    return cpp11::as_sexp(set_observed_data_cpp(cpp11::as_cpp<cpp11::decay_t<SEXP>>(lfmcmc), cpp11::as_cpp<cpp11::decay_t<std::vector< double >>>(observed_data_)));
   END_CPP11
 }
 // lfmcmc.cpp

src/lfmcmc.cpp

@@ -9,7 +9,7 @@
 
 using namespace epiworld;
 
-#define TData_default std::vector< int >
+#define TData_default std::vector< double >
 
 #define WrapLFMCMC(a) \
     cpp11::external_pointer<LFMCMC<TData_default>> (a)
@@ -29,8 +29,13 @@ SEXP LFMCMC_cpp(
         new LFMCMC<TData_default>()
     );
 
-    cpp11::external_pointer<Model<int>> modelptr(model);
-    lfmcmc_ptr->set_rand_engine(modelptr->get_rand_endgine());
+    if (Rf_inherits(model, "epiworld_model")) {
+        cpp11::external_pointer<Model<int>> modelptr(model);
+        lfmcmc_ptr->set_rand_engine(modelptr->get_rand_endgine());
+    } else {
+        auto new_ptr = std::make_shared<std::mt19937>(std::mt19937());
+        lfmcmc_ptr->set_rand_engine(new_ptr);
+    }
 
     return lfmcmc_ptr;
 }
@@ -42,9 +47,9 @@ SEXP LFMCMC_cpp(
 [[cpp11::register]]
 SEXP run_lfmcmc_cpp(
     SEXP lfmcmc,
-    std::vector<epiworld_double> params_init_,
+    std::vector<double> params_init_,
     size_t n_samples_,
-    epiworld_double epsilon_,
+    double epsilon_,
     int seed
 ) {
     WrapLFMCMC(lfmcmc_ptr)(lfmcmc);
@@ -60,7 +65,7 @@ SEXP run_lfmcmc_cpp(
 [[cpp11::register]]
 SEXP set_observed_data_cpp(
     SEXP lfmcmc,
-    std::vector< int > observed_data_
+    std::vector< double > observed_data_
 ) {
     WrapLFMCMC(lfmcmc_ptr)(lfmcmc);
     lfmcmc_ptr->set_observed_data(observed_data_);
@@ -105,7 +110,7 @@ SEXP use_proposal_norm_reflective_cpp(
     SEXP lfmcmc
 ) {
     WrapLFMCMC(lfmcmc_ptr)(lfmcmc);
-    lfmcmc_ptr->set_proposal_fun(make_proposal_norm_reflective<std::vector<int>>(.5, 0, 1));
+    lfmcmc_ptr->set_proposal_fun(make_proposal_norm_reflective<std::vector<double>>(.5, 0, 1));
     return lfmcmc;
 }
 
@@ -123,7 +128,7 @@ SEXP set_simulation_fun_cpp(
         auto params_doubles = cpp11::doubles(params);
 
         return cpp11::as_cpp<TData_default>(
-            cpp11::integers(fun(params_doubles))
+            cpp11::doubles(fun(params_doubles))
             );
     };
 
@@ -146,8 +151,8 @@ SEXP set_summary_fun_cpp(
         LFMCMC<TData_default>*
         ) -> void {
 
-        auto dat_int = cpp11::integers(dat);
-        auto res_tmp = cpp11::integers(fun(dat_int));
+        auto dat_int = cpp11::doubles(dat);
+        auto res_tmp = cpp11::doubles(fun(dat_int));
 
         if (res.size() == 0u)
             res.resize(res_tmp.size());
@@ -198,7 +203,7 @@ SEXP use_kernel_fun_gaussian_cpp(
     SEXP lfmcmc
 ) {
     WrapLFMCMC(lfmcmc_ptr)(lfmcmc);
-    lfmcmc_ptr->set_kernel_fun(kernel_fun_gaussian<std::vector<int>>);
+    lfmcmc_ptr->set_kernel_fun(kernel_fun_gaussian<std::vector<double>>);
     return lfmcmc;
 }
 

vignettes/likelihood-free-mcmc.Rmd

@@ -32,7 +32,7 @@ model_seed <- 122
 model_sir <- ModelSIR(
   name = "COVID-19",
   prevalence = .1,
-  transmission_rate = .9,
+  transmission_rate = .3,
   recovery_rate = .3
 )
 
@@ -58,7 +58,7 @@ summary(model_sir)
 ## Setup LFMCMC
 ```{r lfmcmc-setup}
 # 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 LFMCMC simulation function
 simfun <- function(params) {
@@ -71,8 +71,8 @@ simfun <- function(params) {
     ndays = 50
   )
 
-  res <- unname(as.integer(get_today_total(model_sir)))
-  return(res)
+  get_today_total(model_sir)
+
 }
 
 # Define the LFMCMC summary function
@@ -83,7 +83,7 @@ sumfun <- function(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), )
+  res <- plogis(qlogis(params_prev) + rnorm(length(params_prev)))
   return(res)
 }
 
@@ -95,7 +95,7 @@ kernelfun <- function(stats_now, stats_obs, epsilon) {
     stats_obs,
     stats_now))
 
-  return(ifelse(sqrt(ans) < epsilon, 1.0, 0.0))
+  dnorm(sqrt(sum((stats_now - stats_obs)^2)))
 }
 
 # Create the LFMCMC model