diff --git a/.devcontainer/Dockerfile b/.devcontainer/Dockerfile
index 1c71422e..8f990d4c 100644
--- a/.devcontainer/Dockerfile
+++ b/.devcontainer/Dockerfile
@@ -10,4 +10,6 @@ RUN install2.r cpp11 roxygen2 tinytest data.table netplot \
 
 RUN install2.r languageserver
 
+RUN apt-get update && apt-get install --no-install-recommends -y valgrind gdb
+
 CMD ["bash"]
diff --git a/.gitignore b/.gitignore
index bc64c2f6..eb777ffc 100644
--- a/.gitignore
+++ b/.gitignore
@@ -20,3 +20,6 @@ src/Makevars
 images
 inst/doc
 docs
+
+config.log
+config.status
diff --git a/.vscode/settings.json b/.vscode/settings.json
index b4549178..296aa280 100644
--- a/.vscode/settings.json
+++ b/.vscode/settings.json
@@ -54,7 +54,8 @@
         "stdexcept": "cpp",
         "streambuf": "cpp",
         "typeinfo": "cpp",
-        "thread": "cpp"
+        "thread": "cpp",
+        "cinttypes": "cpp"
     },
     "editor.indentSize": "tabSize",
     "[r]": {
diff --git a/DESCRIPTION b/DESCRIPTION
index e98c661e..bd48843f 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: epiworldR
 Type: Package
 Title: Fast Agent-Based Epi Models
-Version: 0.3-2
+Version: 0.4-3
 Authors@R: c(
   person(given="George", family="Vega Yon", role=c("aut","cre"),
   email="g.vegayon@gmail.com", comment = c(ORCID = "0000-0002-3171-0844")),
diff --git a/Makefile b/Makefile
index 7e0f0bbd..4dfbef48 100644
--- a/Makefile
+++ b/Makefile
@@ -46,9 +46,14 @@ clean:
 	sed -i -E 's/^library\(epiworldRdev\)/library(epiworldR)/g' README.*
 
 docs:
-	Rscript --vanilla -e 'roxygen2::roxigenize()'
-
-.PHONY: build update check clean docs docker-debug
+	Rscript --vanilla -e 'roxygen2::roxygenize()'
 
 checkv: build
 	R CMD check --as-cran --use-valgrind epiworldR*.tar.gz
+
+# Builds and installs without vignettes
+dev: clean
+	R CMD build --no-build-vignettes .
+	R CMD INSTALL epiworldR_$(VERSION).tar.gz
+
+.PHONY: build update check clean docs docker-debug dev
diff --git a/NAMESPACE b/NAMESPACE
index 17f23b74..316e5938 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -62,6 +62,7 @@ S3method(print,epiworld_agents_tools)
 S3method(print,epiworld_entities)
 S3method(print,epiworld_entity)
 S3method(print,epiworld_globalevent)
+S3method(print,epiworld_lfmcmc)
 S3method(print,epiworld_model)
 S3method(print,epiworld_saver)
 S3method(print,epiworld_tool)
@@ -73,13 +74,22 @@ S3method(queuing_on,epiworld_model)
 S3method(queuing_on,epiworld_seirconn)
 S3method(queuing_on,epiworld_sirconn)
 S3method(run,epiworld_model)
+S3method(run_lfmcmc,epiworld_lfmcmc)
 S3method(run_multiple,epiworld_model)
+S3method(set_kernel_fun,epiworld_lfmcmc)
 S3method(set_name,epiworld_model)
+S3method(set_observed_data,epiworld_lfmcmc)
+S3method(set_par_names,epiworld_lfmcmc)
 S3method(set_param,epiworld_model)
+S3method(set_proposal_fun,epiworld_lfmcmc)
+S3method(set_simulation_fun,epiworld_lfmcmc)
+S3method(set_stats_names,epiworld_lfmcmc)
+S3method(set_summary_fun,epiworld_lfmcmc)
 S3method(size,epiworld_model)
 S3method(summary,epiworld_model)
 S3method(verbose_off,epiworld_model)
 S3method(verbose_on,epiworld_model)
+export(LFMCMC)
 export(ModelDiffNet)
 export(ModelSEIR)
 export(ModelSEIRCONN)
@@ -167,6 +177,7 @@ export(rm_entity)
 export(rm_tool)
 export(rm_virus)
 export(run)
+export(run_lfmcmc)
 export(run_multiple)
 export(run_multiple_get_results)
 export(set_agents_data)
@@ -179,9 +190,12 @@ export(set_distribution_virus)
 export(set_incubation)
 export(set_incubation_fun)
 export(set_incubation_ptr)
+export(set_kernel_fun)
 export(set_name)
 export(set_name_tool)
 export(set_name_virus)
+export(set_observed_data)
+export(set_par_names)
 export(set_param)
 export(set_prob_death)
 export(set_prob_death_fun)
@@ -192,9 +206,13 @@ export(set_prob_infecting_ptr)
 export(set_prob_recovery)
 export(set_prob_recovery_fun)
 export(set_prob_recovery_ptr)
+export(set_proposal_fun)
 export(set_recovery_enhancer)
 export(set_recovery_enhancer_fun)
 export(set_recovery_enhancer_ptr)
+export(set_simulation_fun)
+export(set_stats_names)
+export(set_summary_fun)
 export(set_susceptibility_reduction)
 export(set_susceptibility_reduction_fun)
 export(set_susceptibility_reduction_ptr)
@@ -204,6 +222,8 @@ export(set_transmission_reduction_ptr)
 export(size)
 export(tool)
 export(tool_fun_logit)
+export(use_kernel_fun_gaussian)
+export(use_proposal_norm_reflective)
 export(verbose_off)
 export(verbose_on)
 export(virus)
diff --git a/R/LFMCMC.R b/R/LFMCMC.R
new file mode 100644
index 00000000..efcc8a8c
--- /dev/null
+++ b/R/LFMCMC.R
@@ -0,0 +1,216 @@
+#' Likelihood-Free Markhov Chain Monte Carlo (LFMCMC)
+#'
+#'
+#' @aliases epiworld_lfmcmc
+#' @details
+#' Performs a Likelihood-Free Markhov Chain Monte Carlo simulation
+#' @param model A model of class [epiworld_model]
+#' @returns
+#' The `LFMCMC` function returns a model of class [epiworld_lfmcmc].
+#' @examples
+#' ## Setup an SIR model to use in the simulation
+#' model_seed <- 122
+#' 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)
+#'
+#' ## Setup LFMCMC
+#' # Extract the observed data from the model
+#' obs_data <- unname(as.integer(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)))
+#'   return(res)
+#' }
+#'
+#' # Define the summary function
+#' sumfun <- function(dat) {
+#'   return(dat)
+#' }
+#'
+#' # Create the LFMCMC model
+#' lfmcmc_model <- LFMCMC(model_sir) |>
+#'   set_simulation_fun(simfun) |>
+#'   set_summary_fun(sumfun) |>
+#'   use_proposal_norm_reflective() |>
+#'   use_kernel_fun_gaussian() |>
+#'   set_observed_data(obs_data)
+#'
+#' ## Run LFMCMC simulation
+#' # Set initial parameters
+#' par0 <- as.double(c(0.1, 0.5))
+#' n_samp <- 2000
+#' epsil <- as.double(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)
+#' @export
+LFMCMC <- function(model) {
+  if (!inherits(model, "epiworld_model"))
+    stop("model should be of class 'epiworld_model'. It is of class ", class(model))
+
+  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 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")
+
+#' @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))
+  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")
+
+#' @export
+set_observed_data.epiworld_lfmcmc <- function(lfmcmc, observed_data_) {
+  set_observed_data_cpp(lfmcmc, observed_data_)
+  invisible(lfmcmc)
+}
+
+#' @rdname LFMCMC
+#' @param lfmcmc LFMCMC model
+#' @param fun The LFMCMC proposal function
+#' @returns The lfmcmc model with the proposal function added
+#' @export
+set_proposal_fun <- function(lfmcmc, fun) UseMethod("set_proposal_fun")
+
+#' @export
+set_proposal_fun.epiworld_lfmcmc <- function(lfmcmc, fun) {
+  set_proposal_fun_cpp(lfmcmc, fun)
+  invisible(lfmcmc)
+}
+
+#' @rdname LFMCMC
+#' @param lfmcmc The LFMCMC model
+#' @returns The LFMCMC model with proposal function set to norm reflective
+#' @export
+use_proposal_norm_reflective <- function(lfmcmc) {
+  use_proposal_norm_reflective_cpp(lfmcmc)
+  invisible(lfmcmc)
+}
+
+#' @rdname LFMCMC
+#' @param lfmcmc LFMCMC model
+#' @param fun The LFMCMC simulation function
+#' @returns The lfmcmc model with the simulation function added
+#' @export
+set_simulation_fun <- function(lfmcmc, fun) UseMethod("set_simulation_fun")
+
+#' @export
+set_simulation_fun.epiworld_lfmcmc <- function(lfmcmc, fun) {
+  set_simulation_fun_cpp(lfmcmc, fun)
+  invisible(lfmcmc)
+}
+
+#' @rdname LFMCMC
+#' @param lfmcmc LFMCMC model
+#' @param fun The LFMCMC sumamry function
+#' @returns The lfmcmc model with the summary function added
+#' @export
+set_summary_fun <- function(lfmcmc, fun) UseMethod("set_summary_fun")
+
+#' @export
+set_summary_fun.epiworld_lfmcmc <- function(lfmcmc, fun) {
+  set_summary_fun_cpp(lfmcmc, fun)
+  invisible(lfmcmc)
+}
+
+#' @rdname LFMCMC
+#' @param lfmcmc LFMCMC model
+#' @param fun The LFMCMC kernel function
+#' @returns The lfmcmc model with the kernel function added
+#' @export
+set_kernel_fun <- function(lfmcmc, fun) UseMethod("set_kernel_fun")
+
+#' @export
+set_kernel_fun.epiworld_lfmcmc <- function(lfmcmc, fun) {
+  set_kernel_fun_cpp(lfmcmc, fun)
+  invisible(lfmcmc)
+}
+
+#' @rdname LFMCMC
+#' @param lfmcmc The LFMCMC model
+#' @returns The LFMCMC model with kernel function set to gaussian
+#' @export
+use_kernel_fun_gaussian <- function(lfmcmc) {
+  use_kernel_fun_gaussian_cpp(lfmcmc)
+  invisible(lfmcmc)
+}
+
+#' @rdname LFMCMC
+#' @param lfmcmc LFMCMC model
+#' @param names The model parameter names
+#' @returns The lfmcmc model with the parameter names added
+#' @export
+set_par_names <- function(lfmcmc, names) UseMethod("set_par_names")
+
+#' @export
+set_par_names.epiworld_lfmcmc <- function(lfmcmc, names) {
+  set_par_names_cpp(lfmcmc, names)
+  invisible(lfmcmc)
+}
+
+#' @rdname LFMCMC
+#' @param lfmcmc LFMCMC model
+#' @param names The model stats names
+#' @returns The lfmcmc model with the stats names added
+#' @export
+set_stats_names <- function(lfmcmc, names) UseMethod("set_stats_names")
+
+#' @export
+set_stats_names.epiworld_lfmcmc <- function(lfmcmc, names) {
+  set_stats_names_cpp(lfmcmc, names)
+  invisible(lfmcmc)
+}
+
+#' @rdname LFMCMC
+#' @param x LFMCMC model to print
+#' @param ... Ignored
+#' @returns The lfmcmc model
+#' @export
+print.epiworld_lfmcmc <- function(x, ...) {
+  print_lfmcmc_cpp(x)
+  invisible(x)
+}
diff --git a/R/cpp11.R b/R/cpp11.R
index b1276fdb..d699a510 100644
--- a/R/cpp11.R
+++ b/R/cpp11.R
@@ -220,6 +220,54 @@ ModelSEIRMixing_cpp <- function(name, n, prevalence, contact_rate, transmission_
   .Call(`_epiworldR_ModelSEIRMixing_cpp`, name, n, prevalence, contact_rate, transmission_rate, incubation_days, recovery_rate, contact_matrix)
 }
 
+LFMCMC_cpp <- function(model) {
+  .Call(`_epiworldR_LFMCMC_cpp`, model)
+}
+
+run_lfmcmc_cpp <- function(lfmcmc, params_init_, n_samples_, epsilon_, seed) {
+  .Call(`_epiworldR_run_lfmcmc_cpp`, lfmcmc, params_init_, n_samples_, epsilon_, seed)
+}
+
+set_observed_data_cpp <- function(lfmcmc, observed_data_) {
+  .Call(`_epiworldR_set_observed_data_cpp`, lfmcmc, observed_data_)
+}
+
+set_proposal_fun_cpp <- function(lfmcmc, fun) {
+  .Call(`_epiworldR_set_proposal_fun_cpp`, lfmcmc, fun)
+}
+
+use_proposal_norm_reflective_cpp <- function(lfmcmc) {
+  .Call(`_epiworldR_use_proposal_norm_reflective_cpp`, lfmcmc)
+}
+
+set_simulation_fun_cpp <- function(lfmcmc, fun) {
+  .Call(`_epiworldR_set_simulation_fun_cpp`, lfmcmc, fun)
+}
+
+set_summary_fun_cpp <- function(lfmcmc, fun) {
+  .Call(`_epiworldR_set_summary_fun_cpp`, lfmcmc, fun)
+}
+
+set_kernel_fun_cpp <- function(lfmcmc, fun) {
+  .Call(`_epiworldR_set_kernel_fun_cpp`, lfmcmc, fun)
+}
+
+use_kernel_fun_gaussian_cpp <- function(lfmcmc) {
+  .Call(`_epiworldR_use_kernel_fun_gaussian_cpp`, lfmcmc)
+}
+
+set_par_names_cpp <- function(lfmcmc, names) {
+  .Call(`_epiworldR_set_par_names_cpp`, lfmcmc, names)
+}
+
+set_stats_names_cpp <- function(lfmcmc, names) {
+  .Call(`_epiworldR_set_stats_names_cpp`, lfmcmc, names)
+}
+
+print_lfmcmc_cpp <- function(lfmcmc) {
+  .Call(`_epiworldR_print_lfmcmc_cpp`, lfmcmc)
+}
+
 print_cpp <- function(m, lite) {
   .Call(`_epiworldR_print_cpp`, m, lite)
 }
diff --git a/epiworldR.Rproj b/epiworldR.Rproj
index eaa6b818..d9f3f4e4 100644
--- a/epiworldR.Rproj
+++ b/epiworldR.Rproj
@@ -14,5 +14,6 @@ LaTeX: pdfLaTeX
 
 BuildType: Package
 PackageUseDevtools: Yes
+PackageCleanBeforeInstall: No
 PackageInstallArgs: --no-multiarch --with-keep.source
 PackageRoxygenize: rd,collate,namespace
diff --git a/inst/include/epiworld/database-bones.hpp b/inst/include/epiworld/database-bones.hpp
index 3d64ea00..152ee1a4 100644
--- a/inst/include/epiworld/database-bones.hpp
+++ b/inst/include/epiworld/database-bones.hpp
@@ -178,6 +178,10 @@ class DataBase {
         std::vector< int > & counts
     ) const;
 
+    void get_today_transition_matrix(
+        std::vector< int > & counts
+    ) const;
+
     void get_hist_total(
         std::vector< int > * date,
         std::vector< std::string > * state,
@@ -246,10 +250,13 @@ class DataBase {
         std::string fn_generation_time
         ) const;
     
+    /***
+     * @brief Record a transmission event
+     */
     void record_transmission(int i, int j, int virus, int i_expo_date);
 
-    size_t get_n_viruses() const;
-    size_t get_n_tools() const;
+    size_t get_n_viruses() const; ///< Get the number of viruses
+    size_t get_n_tools() const; ///< Get the number of tools
     
     void set_user_data(std::vector< std::string > names);
     void add_user_data(std::vector< epiworld_double > x);
@@ -288,7 +295,11 @@ class DataBase {
 
     /**
      * Calculates the generating time
-     * @param agent_id,virus_id,time,gentime vectors where to save the values agent_id
+     * @param agent_id,virus_id,time,gentime vectors where to save the values 
+     * 
+     * @details
+     * The generation time is the time between the infection of the source and 
+     * the infection of the target.
     */
    ///@{
     void generation_time(
@@ -296,11 +307,11 @@ class DataBase {
         std::vector< int > & virus_id,
         std::vector< int > & time,
         std::vector< int > & gentime
-    ) const;
+    ) const; ///< Get the generation time
 
     void generation_time(
         std::string fn
-    ) const;
+    ) const; ///< Write the generation time to a file
     ///@}
 
 };
diff --git a/inst/include/epiworld/database-meat.hpp b/inst/include/epiworld/database-meat.hpp
index 95e641bb..f6d358db 100644
--- a/inst/include/epiworld/database-meat.hpp
+++ b/inst/include/epiworld/database-meat.hpp
@@ -639,6 +639,18 @@ inline void DataBase<TSeq>::get_hist_tool(
 
 }
 
+template<typename TSeq>
+inline void DataBase<TSeq>::get_today_transition_matrix(
+    std::vector< int > & counts
+) const
+{
+
+    counts = transition_matrix;
+
+    return;
+
+}
+
 template<typename TSeq>
 inline void DataBase<TSeq>::get_hist_transition_matrix(
     std::vector< std::string > & state_from,
@@ -995,9 +1007,9 @@ inline void DataBase<TSeq>::write_data(
                         #ifdef EPI_DEBUG
                         EPI_GET_THREAD_ID() << " " <<
                         #endif
-                        i << " " <<
-                        model->states_labels[from] << " " <<
-                        model->states_labels[to] << " " <<
+                        i << " \"" <<
+                        model->states_labels[from] << "\" \"" <<
+                        model->states_labels[to] << "\" " <<
                         hist_transition_matrix[i * (ns * ns) + to * ns + from] << "\n";
                 
         }
diff --git a/inst/include/epiworld/epiworld.hpp b/inst/include/epiworld/epiworld.hpp
index 7e8f6624..6a060f7e 100644
--- a/inst/include/epiworld/epiworld.hpp
+++ b/inst/include/epiworld/epiworld.hpp
@@ -18,8 +18,8 @@
 
 /* Versioning */
 #define EPIWORLD_VERSION_MAJOR 0
-#define EPIWORLD_VERSION_MINOR 3
-#define EPIWORLD_VERSION_PATCH 2
+#define EPIWORLD_VERSION_MINOR 4
+#define EPIWORLD_VERSION_PATCH 3
 
 static const int epiworld_version_major = EPIWORLD_VERSION_MAJOR;
 static const int epiworld_version_minor = EPIWORLD_VERSION_MINOR;
@@ -33,7 +33,7 @@ namespace epiworld {
     #include "misc.hpp"
     #include "progress.hpp"
 
-    // #include "math/summary-stats.hpp"
+    #include "math/distributions.hpp"
 
     #include "math/lfmcmc.hpp"
 
diff --git a/inst/include/epiworld/math/distributions.hpp b/inst/include/epiworld/math/distributions.hpp
new file mode 100644
index 00000000..b32a012e
--- /dev/null
+++ b/inst/include/epiworld/math/distributions.hpp
@@ -0,0 +1,168 @@
+#ifndef EPIWORLD_MATH_DISTRIBUTIONS_HPP
+#define EPIWORLD_MATH_DISTRIBUTIONS_HPP
+
+// Implementing the factorial function
+/**
+ * @brief Compute the log of the factorial
+ * 
+ * @param n Number
+ * 
+ * @return The log of the factorial
+ */
+inline double log_factorial(int n)
+{
+    if (n == 0)
+        return 0.0;
+    return std::log(static_cast<double>(n)) + log_factorial(n-1);
+}
+
+/**
+ * @brief Compute the Poisson probability
+ * 
+ * @param k Number of events
+ * @param lambda Rate
+ * @param max_n Maximum number of events
+ * @param as_log Return the log of the probability
+ * 
+ * @return The Poisson probability
+ */
+inline double dpois(
+    int k,
+    double lambda,
+    int max_n = 100,
+    bool as_log = false
+    )
+{
+
+    if (max_n < k)
+        throw std::runtime_error("max_n must be greater than k");
+
+    double res = k * std::log(lambda) - lambda - log_factorial(
+        std::min(k, max_n)
+        );
+    
+    return as_log ? res : std::exp(res);
+}
+
+/**
+ * @brief Compute the probability of the generation interval
+ * 
+ * @details
+ * If `p_0_approx` is negative, it will be computed using the Poisson
+ * distribution. If `normalizing` is negative, it will be computed on the fly
+ * 
+ * @param g Generation interval
+ * @param S Population size
+ * @param p_c Probability of contact
+ * @param p_i Probability of infection
+ * @param p_r Probability of recovery
+ * @param p_0_approx Approximation of the probability of not being infected
+ * @param normalizing Normalizing constant
+ * @param max_contacts Maximum number of contacts
+ * @param max_days Maximum number of days
+ * 
+ * @return The probability of the generation interval
+ * 
+ */
+inline double dgenint(
+    int g,
+    double S,
+    double p_c,
+    double p_i,
+    double p_r,
+    double & p_0_approx,
+    double & normalizing,
+    int max_contacts = 200,
+    int max_days = 200
+    ) {
+
+    if ((g < 1) || (g > max_days))
+        return 0.0;
+
+    if (p_0_approx < 0.0)
+    {
+
+        p_0_approx = 0.0;
+        for (int i = 0; i < max_contacts; ++i)
+        {
+
+            p_0_approx += std::exp(
+                dpois(i, S * p_c, max_contacts, true) +
+                std::log(1.0 - p_i) * static_cast<double>(i)
+                );
+
+        }
+    }
+
+    double g_dbl = static_cast<double>(g);
+
+    if (normalizing < 0.0)
+    {
+
+        normalizing = 1.0;
+        double log1_p_r = std::log(1.0 - p_r);
+        double log_p_r = std::log(p_r);
+        double log_p_0_approx = std::log(p_0_approx);
+        for (int i = 1; i <= max_days; ++i)
+        {
+
+            double i_dbl = static_cast<double>(i);
+
+            normalizing -= std::exp(
+                log1_p_r * (i_dbl - 1.0) +
+                log_p_r +
+                log_p_0_approx * (i_dbl - 1.0)
+                );
+        }
+
+    }
+
+
+    return std::exp(
+        std::log(1 - p_r) * (g_dbl)+
+        std::log(p_0_approx) * (g_dbl - 1.0) +
+        std::log(1.0 - p_0_approx) -
+        std::log(normalizing)
+        );
+
+}
+
+// Mean of the generation interval
+/**
+ * @brief Compute the mean of the generation interval
+ * @param S Population size.
+ * @param p_c Probability of contact.
+ * @param p_i Probability of infection.
+ * @param p_r Probability of recovery.
+ * @param max_days Maximum number of days.
+ * @param max_n Maximum number of contacts.
+ * 
+ * @return The mean of the generation interval
+ */
+inline double gen_int_mean(
+    double S,
+    double p_c,
+    double p_i,
+    double p_r,
+    int max_days = 200,
+    int max_n = 200
+    ) {
+
+    double mean = 0.0;
+    double p_0_approx = -1.0;
+    double normalizing = -1.0;
+    for (int i = 1; i < max_days; ++i)
+    {
+        mean += 
+            static_cast<double>(i) *
+            dgenint(
+                i, S, p_c, p_i, p_r, p_0_approx, normalizing, max_n, max_days
+                );
+
+    }
+
+    return mean;
+
+}
+
+#endif
\ No newline at end of file
diff --git a/inst/include/epiworld/math/lfmcmc/lfmcmc-bones.hpp b/inst/include/epiworld/math/lfmcmc/lfmcmc-bones.hpp
index 6c494b95..7002aa24 100755
--- a/inst/include/epiworld/math/lfmcmc/lfmcmc-bones.hpp
+++ b/inst/include/epiworld/math/lfmcmc/lfmcmc-bones.hpp
@@ -116,7 +116,7 @@ class LFMCMC {
 private:
 
     // Random number sampling
-    std::mt19937 * engine = nullptr;
+    std::shared_ptr< std::mt19937 > engine = nullptr;
     
     std::shared_ptr< std::uniform_real_distribution<> > runifd =
         std::make_shared< std::uniform_real_distribution<> >(0.0, 1.0);
@@ -128,7 +128,7 @@ class LFMCMC {
         std::make_shared< std::gamma_distribution<> >();
 
     // Process data
-    TData * observed_data;
+    TData observed_data;
     
     // Information about the size of the problem
     size_t n_samples;
@@ -187,14 +187,15 @@ class LFMCMC {
     void run(
         std::vector< epiworld_double > param_init,
         size_t n_samples_,
-        epiworld_double epsilon_
+        epiworld_double epsilon_,
+        int seed = -1
         );
 
     LFMCMC() {};
-    LFMCMC(TData & observed_data_) : observed_data(&observed_data_) {};
+    LFMCMC(const TData & observed_data_) : observed_data(observed_data_) {};
     ~LFMCMC() {};
 
-    void set_observed_data(TData & observed_data_) {observed_data = &observed_data_;};
+    void set_observed_data(const TData & observed_data_) {observed_data = observed_data_;};
     void set_proposal_fun(LFMCMCProposalFun<TData> fun);
     void set_simulation_fun(LFMCMCSimFun<TData> fun);
     void set_summary_fun(LFMCMCSummaryFun<TData> fun);
@@ -206,8 +207,8 @@ class LFMCMC {
      * @param eng 
      */
     ///@{
-    void set_rand_engine(std::mt19937 & eng);
-    std::mt19937 & get_rand_endgine();
+    void set_rand_engine(std::shared_ptr< std::mt19937 > & eng);
+    std::shared_ptr< std::mt19937 > & get_rand_endgine();
     void seed(epiworld_fast_uint s);
     void set_rand_gamma(epiworld_double alpha, epiworld_double beta);
     epiworld_double runif();
diff --git a/inst/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp b/inst/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp
index 38cbf56b..363fb766 100755
--- a/inst/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp
+++ b/inst/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp
@@ -52,7 +52,7 @@ inline void LFMCMC<TData>::print()
     printf_epiworld("___________________________________________\n\n");
     printf_epiworld("LIKELIHOOD-FREE MARKOV CHAIN MONTE CARLO\n\n");
 
-    printf_epiworld("N Samples : %ld\n", n_samples);
+    printf_epiworld("N Samples : %zu\n", n_samples);
 
     std::string abbr;
     epiworld_double elapsed;
diff --git a/inst/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp b/inst/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp
index 688f4fd0..c7565121 100755
--- a/inst/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp
+++ b/inst/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp
@@ -205,7 +205,8 @@ template<typename TData>
 inline void LFMCMC<TData>::run(
     std::vector< epiworld_double > params_init_,
     size_t n_samples_,
-    epiworld_double epsilon_
+    epiworld_double epsilon_,
+    int seed
     )
 {
 
@@ -218,6 +219,9 @@ inline void LFMCMC<TData>::run(
     params_init  = params_init_;
     n_parameters = params_init_.size();
 
+    if (seed >= 0)
+        this->seed(seed);
+
     params_now.resize(n_parameters);
     params_prev.resize(n_parameters);
 
@@ -228,7 +232,7 @@ inline void LFMCMC<TData>::run(
     params_now  = params_init;
 
     // Computing the baseline sufficient statistics
-    summary_fun(observed_stats, *observed_data, this);
+    summary_fun(observed_stats, observed_data, this);
     n_statistics = observed_stats.size();
 
     // Reserving size
@@ -378,9 +382,9 @@ inline void LFMCMC<TData>::seed(epiworld_fast_uint s) {
 }
 
 template<typename TData>
-inline void LFMCMC<TData>::set_rand_engine(std::mt19937 & eng)
+inline void LFMCMC<TData>::set_rand_engine(std::shared_ptr< std::mt19937 > & eng)
 {
-    engine = &eng;
+    engine = eng;
 }
 
 template<typename TData>
@@ -390,9 +394,9 @@ inline void LFMCMC<TData>::set_rand_gamma(epiworld_double alpha, epiworld_double
 }
 
 template<typename TData>
-inline std::mt19937 & LFMCMC<TData>::get_rand_endgine()
+inline std::shared_ptr< std::mt19937 > & LFMCMC<TData>::get_rand_endgine()
 {
-    return *engine;
+    return engine;
 }
 
 // Step 1: Simulate data
diff --git a/inst/include/epiworld/model-bones.hpp b/inst/include/epiworld/model-bones.hpp
index 4029f12a..d4055094 100644
--- a/inst/include/epiworld/model-bones.hpp
+++ b/inst/include/epiworld/model-bones.hpp
@@ -140,7 +140,7 @@ class Model {
     std::vector< Entity<TSeq> > entities = {}; 
     std::vector< Entity<TSeq> > entities_backup = {};
 
-    std::mt19937 engine;
+    std::shared_ptr< std::mt19937 > engine = std::make_shared< std::mt19937 >();
     
     std::uniform_real_distribution<> runifd      =
         std::uniform_real_distribution<> (0.0, 1.0);
@@ -261,17 +261,6 @@ class Model {
 
     virtual ~Model() {};
 
-    void clone_population(
-        std::vector< Agent<TSeq> > & other_population,
-        std::vector< Entity<TSeq> > & other_entities,
-        Model<TSeq> * other_model,
-        bool & other_directed
-    ) const ;
-
-    void clone_population(
-        const Model<TSeq> & other_model
-    );
-
     /**
      * @name Set the backup object
      * @details `backup` can be used to restore the entire object
@@ -285,6 +274,7 @@ class Model {
     ///@}
 
     DataBase<TSeq> & get_db();
+    const DataBase<TSeq> & get_db() const;
     epiworld_double & operator()(std::string pname);
 
     size_t size() const;
@@ -296,8 +286,8 @@ class Model {
      * @param s Seed
      */
     ///@{
-    void set_rand_engine(std::mt19937 & eng);
-    std::mt19937 & get_rand_endgine();
+    void set_rand_engine(std::shared_ptr< std::mt19937 > & eng);
+    std::shared_ptr< std::mt19937 > & get_rand_endgine();
     void seed(size_t s);
     void set_rand_norm(epiworld_double mean, epiworld_double sd);
     void set_rand_unif(epiworld_double a, epiworld_double b);
@@ -609,7 +599,7 @@ class Model {
     // void set_param(size_t k, epiworld_double val);
     void set_param(std::string pname, epiworld_double val);
     // epiworld_double par(epiworld_fast_uint k);
-    epiworld_double par(std::string pname);
+    epiworld_double par(std::string pname) const;
     ///@}
 
     void get_elapsed(
diff --git a/inst/include/epiworld/model-meat.hpp b/inst/include/epiworld/model-meat.hpp
index d8d4e978..5df3befb 100644
--- a/inst/include/epiworld/model-meat.hpp
+++ b/inst/include/epiworld/model-meat.hpp
@@ -495,13 +495,6 @@ inline Model<TSeq> & Model<TSeq>::operator=(const Model<TSeq> & m)
         for (auto & p : population_backup)
             p.model = this;
 
-    for (auto & e : entities)
-        e.model = this;
-
-    if (entities_backup.size() != 0)
-        for (auto & e : entities_backup)
-            e.model = this;
-
     db = m.db;
     db.model = this;
     db.user_data.model = this;
@@ -566,6 +559,13 @@ inline DataBase<TSeq> & Model<TSeq>::get_db()
     return db;
 }
 
+template<typename TSeq>
+inline const DataBase<TSeq> & Model<TSeq>::get_db() const
+{
+    return db;
+}
+
+
 template<typename TSeq>
 inline std::vector<Agent<TSeq>> & Model<TSeq>::get_agents()
 {
@@ -674,12 +674,6 @@ inline void Model<TSeq>::agents_empty_graph(
 
 }
 
-// template<typename TSeq>
-// inline void Model<TSeq>::set_rand_engine(std::mt19937 & eng)
-// {
-//     engine = std::make_shared< std::mt19937 >(eng);
-// }
-
 template<typename TSeq>
 inline void Model<TSeq>::set_rand_gamma(epiworld_double alpha, epiworld_double beta)
 {
@@ -820,7 +814,7 @@ inline void Model<TSeq>::set_backup()
 // }
 
 template<typename TSeq>
-inline std::mt19937 & Model<TSeq>::get_rand_endgine()
+inline std::shared_ptr< std::mt19937 > & Model<TSeq>::get_rand_endgine()
 {
     return engine;
 }
@@ -828,86 +822,86 @@ inline std::mt19937 & Model<TSeq>::get_rand_endgine()
 template<typename TSeq>
 inline epiworld_double Model<TSeq>::runif() {
     // CHECK_INIT()
-    return runifd(engine);
+    return runifd(*engine);
 }
 
 template<typename TSeq>
 inline epiworld_double Model<TSeq>::runif(epiworld_double a, epiworld_double b) {
     // CHECK_INIT()
-    return runifd(engine) * (b - a) + a;
+    return runifd(*engine) * (b - a) + a;
 }
 
 template<typename TSeq>
 inline epiworld_double Model<TSeq>::rnorm() {
     // CHECK_INIT()
-    return rnormd(engine);
+    return rnormd(*engine);
 }
 
 template<typename TSeq>
 inline epiworld_double Model<TSeq>::rnorm(epiworld_double mean, epiworld_double sd) {
     // CHECK_INIT()
-    return rnormd(engine) * sd + mean;
+    return rnormd(*engine) * sd + mean;
 }
 
 template<typename TSeq>
 inline epiworld_double Model<TSeq>::rgamma() {
-    return rgammad(engine);
+    return rgammad(*engine);
 }
 
 template<typename TSeq>
 inline epiworld_double Model<TSeq>::rgamma(epiworld_double alpha, epiworld_double beta) {
     auto old_param = rgammad.param();
     rgammad.param(std::gamma_distribution<>::param_type(alpha, beta));
-    epiworld_double ans = rgammad(engine);
+    epiworld_double ans = rgammad(*engine);
     rgammad.param(old_param);
     return ans;
 }
 
 template<typename TSeq>
 inline epiworld_double Model<TSeq>::rexp() {
-    return rexpd(engine);
+    return rexpd(*engine);
 }
 
 template<typename TSeq>
 inline epiworld_double Model<TSeq>::rexp(epiworld_double lambda) {
     auto old_param = rexpd.param();
     rexpd.param(std::exponential_distribution<>::param_type(lambda));
-    epiworld_double ans = rexpd(engine);
+    epiworld_double ans = rexpd(*engine);
     rexpd.param(old_param);
     return ans;
 }
 
 template<typename TSeq>
 inline epiworld_double Model<TSeq>::rlognormal() {
-    return rlognormald(engine);
+    return rlognormald(*engine);
 }
 
 template<typename TSeq>
 inline epiworld_double Model<TSeq>::rlognormal(epiworld_double mean, epiworld_double shape) {
     auto old_param = rlognormald.param();
     rlognormald.param(std::lognormal_distribution<>::param_type(mean, shape));
-    epiworld_double ans = rlognormald(engine);
+    epiworld_double ans = rlognormald(*engine);
     rlognormald.param(old_param);
     return ans;
 }
 
 template<typename TSeq>
 inline int Model<TSeq>::rbinom() {
-    return rbinomd(engine);
+    return rbinomd(*engine);
 }
 
 template<typename TSeq>
 inline int Model<TSeq>::rbinom(int n, epiworld_double p) {
     auto old_param = rbinomd.param();
     rbinomd.param(std::binomial_distribution<>::param_type(n, p));
-    epiworld_double ans = rbinomd(engine);
+    epiworld_double ans = rbinomd(*engine);
     rbinomd.param(old_param);
     return ans;
 }
 
 template<typename TSeq>
 inline void Model<TSeq>::seed(size_t s) {
-    this->engine.seed(s);
+    this->engine->seed(s);
 }
 
 template<typename TSeq>
@@ -1297,7 +1291,7 @@ inline Model<TSeq> & Model<TSeq>::run(
     this->ndays = ndays;
 
     if (seed >= 0)
-        engine.seed(seed);
+        engine->seed(seed);
 
     array_double_tmp.resize(std::max(
         size(),
@@ -1410,15 +1404,6 @@ inline void Model<TSeq>::run_multiple(
 
     // Seeds will be reproducible by default
     std::vector< int > seeds_n(nexperiments);
-    // #ifdef EPI_DEBUG
-    // std::fill(
-    //     seeds_n.begin(),
-    //     seeds_n.end(),
-    //     std::floor(
-    //         runif() * static_cast<double>(std::numeric_limits<int>::max())
-    //     )
-    //     );
-    // #else
     for (auto & s : seeds_n)
     {
         s = static_cast<int>(
@@ -2083,9 +2068,12 @@ inline void Model<TSeq>::set_param(std::string pname, epiworld_double value)
 // }
 
 template<typename TSeq>
-inline epiworld_double Model<TSeq>::par(std::string pname)
+inline epiworld_double Model<TSeq>::par(std::string pname) const
 {
-    return parameters[pname];
+    const auto iter = parameters.find(pname);
+    if (iter == parameters.end())
+        throw std::logic_error("The parameter " + pname + " does not exists.");
+    return iter->second;
 }
 
 #define DURCAST(tunit,txtunit) {\
diff --git a/inst/include/epiworld/models/seirconnected.hpp b/inst/include/epiworld/models/seirconnected.hpp
index cd87e6c3..296f2f3c 100644
--- a/inst/include/epiworld/models/seirconnected.hpp
+++ b/inst/include/epiworld/models/seirconnected.hpp
@@ -60,6 +60,16 @@ class ModelSEIRCONN : public epiworld::Model<TSeq>
 
     size_t get_n_infected() const { return infected.size(); }
 
+    /***
+     * @brief Compute expected generation time
+     * @param max_days Maximum number of days.
+     * @param max_contacts Maximum number of contacts.
+     */
+    std::vector< double > generation_time_expected(
+        int max_days = 200,
+        int max_contacts = 200
+    ) const;
+
 };
 
 template<typename TSeq>
@@ -388,4 +398,61 @@ inline ModelSEIRCONN<TSeq> & ModelSEIRCONN<TSeq>::initial_states(
 
 }
 
+template<typename TSeq>
+inline std::vector< double > ModelSEIRCONN<TSeq>::generation_time_expected(
+    int max_days,
+    int max_contacts
+) const
+{
+
+    // Retrieving total counts
+    std::vector< int > h_date;
+    std::vector< std::string > h_state;
+    std::vector< int > h_counts;
+    const auto this_const = dynamic_cast<const ModelSEIRCONN<TSeq> *>(this);
+    this_const->get_db().get_hist_total(
+        &h_date,
+        &h_state,
+        &h_counts
+    );
+
+    // Retrieving information on susceptibles
+    std::vector< double > S(this_const->get_ndays(), 0.0);
+    for (size_t i = 0; i < h_date.size(); ++i)
+    {
+        if (h_state[i] == "Susceptible")
+            S[h_date[i]] += h_counts[i];
+    }
+
+    // Computing the expected number of days in exposed
+    double days_exposed = this_const->par("Avg. Incubation days");
+
+    // The generation time in the SEIR model starts from 2, as agents 
+    // spend at least one day in the exposed state, and 1 day in the 
+    // infectious state before starting transmitting.
+    std::vector< double > gen_times(
+        this_const->get_ndays(), 1.0 + days_exposed
+        );
+        
+    double p_c = this_const->par("Contact rate")/this_const->size();
+    double p_i = this_const->par("Prob. Transmission");
+    double p_r = this_const->par("Prob. Recovery");
+
+    for (size_t i = 0u; i < this_const->get_ndays(); ++i)
+    {
+        gen_times[i] += gen_int_mean(
+            S[i],
+            p_c,
+            p_i,
+            p_r,
+            max_days,
+            max_contacts
+        );
+
+    }
+
+    return gen_times;
+
+}
+
 #endif
diff --git a/inst/include/epiworld/models/sirconnected.hpp b/inst/include/epiworld/models/sirconnected.hpp
index c66984a1..df39bce4 100644
--- a/inst/include/epiworld/models/sirconnected.hpp
+++ b/inst/include/epiworld/models/sirconnected.hpp
@@ -66,6 +66,15 @@ class ModelSIRCONN : public epiworld::Model<TSeq>
         return infected.size();
     }
 
+    /***
+     * @brief Compute expected generation time
+     * @param max_days Maximum number of days.
+     * @param max_contacts Maximum number of contacts.
+     */
+    std::vector< double > generation_time_expected(
+        int max_days = 200,
+        int max_contacts = 200
+    ) const;
 
 };
 
@@ -361,5 +370,54 @@ inline ModelSIRCONN<TSeq> & ModelSIRCONN<TSeq>::initial_states(
 
 }
 
+template<typename TSeq>
+inline std::vector< double > ModelSIRCONN<TSeq>::generation_time_expected(
+    int max_days,
+    int max_contacts
+) const
+{
+
+    // Retrieving total counts
+    std::vector< int > h_date;
+    std::vector< std::string > h_state;
+    std::vector< int > h_counts;
+    const auto this_const = dynamic_cast<const ModelSIRCONN<TSeq> *>(this);
+    this_const->get_db().get_hist_total(
+        &h_date,
+        &h_state,
+        &h_counts
+    );
+
+    // Retrieving information on susceptibles
+    std::vector< double > S(this_const->get_ndays(), 0.0);
+    for (size_t i = 0; i < h_date.size(); ++i)
+    {
+        if (h_state[i] == "Susceptible")
+            S[h_date[i]] += h_counts[i];
+    }
+
+    // The generation time in the SIR model starts from 1, as agents 
+    // spend at least one day in the infected state before starting
+    // transmitting.
+    std::vector< double > gen_times(this_const->get_ndays(), 1.0);
+    double p_c = this_const->par("Contact rate")/this_const->size();
+    double p_i = this_const->par("Transmission rate");
+    double p_r = this_const->par("Recovery rate");
+    for (size_t i = 0u; i < this_const->get_ndays(); ++i)
+    {
+        gen_times[i] = gen_int_mean(
+            S[i],
+            p_c,
+            p_i,
+            p_r,
+            max_days,
+            max_contacts
+        );
+
+    }
+
+    return gen_times;
+
+}
 
 #endif
diff --git a/inst/include/epiworld/models/surveillance.hpp b/inst/include/epiworld/models/surveillance.hpp
index 6dca7e4e..8b1a2889 100644
--- a/inst/include/epiworld/models/surveillance.hpp
+++ b/inst/include/epiworld/models/surveillance.hpp
@@ -227,7 +227,7 @@ inline ModelSURV<TSeq>::ModelSURV(
 
         // How many will we find
         std::binomial_distribution<> bdist(m->size(), m->par("Surveilance prob."));
-        int nsampled = bdist(m->get_rand_endgine());
+        int nsampled = bdist(*m->get_rand_endgine());
 
         int to_go = nsampled + 1;
 
diff --git a/inst/include/epiworld/progress.hpp b/inst/include/epiworld/progress.hpp
index a0cd001c..0cad6137 100644
--- a/inst/include/epiworld/progress.hpp
+++ b/inst/include/epiworld/progress.hpp
@@ -31,9 +31,16 @@ class Progress {
 inline Progress::Progress(int n_, int width_) {
 
 
+    if (n_ < 0)
+        throw std::invalid_argument("n must be greater or equal than 0.");
+
+    if (width_ <= 0)
+        throw std::invalid_argument("width must be greater than 0");
+
     width     = std::max(7, width_ - 7);
     n         = n_;
-    step_size = static_cast<epiworld_double>(width)/static_cast<epiworld_double>(n);
+    step_size = n == 0? width : static_cast<epiworld_double>(width)/
+        static_cast<epiworld_double>(n);
     last_loc  = 0;
     i         = 0;
 
@@ -58,10 +65,9 @@ inline void Progress::next() {
 
     cur_loc = std::floor((++i) * step_size);
 
-
     #ifndef EPI_DEBUG
     for (int j = 0; j < (cur_loc - last_loc); ++j)
-    {
+    { 
         printf_epiworld("|");
     }
     #endif
diff --git a/inst/include/epiworld/random_graph.hpp b/inst/include/epiworld/random_graph.hpp
index dfe6883d..ee13b5f3 100644
--- a/inst/include/epiworld/random_graph.hpp
+++ b/inst/include/epiworld/random_graph.hpp
@@ -5,7 +5,7 @@ class RandGraph {
 
 private:
     std::shared_ptr< std::mt19937 > engine;
-    std::shared_ptr< std::uniform_real_distribution<> > unifd;
+    std::uniform_real_distribution<> unifd;
     int N = 0;
     bool initialized = false;
 
@@ -14,7 +14,7 @@ class RandGraph {
     RandGraph(int N_) : N(N_) {};
 
     void init(int s);
-    void set_rand_engine(std::mt19937 & e);
+    void set_rand_engine(std::shared_ptr< std::mt19937 > & e);
     epiworld_double runif();
 
 };
@@ -27,18 +27,15 @@ inline void RandGraph::init(int s) {
 
     engine->seed(s);
 
-    if (!unifd)
-        unifd = std::make_shared< std::uniform_real_distribution<> >(0, 1);
-
     initialized = true;
 
 
 }
 
-inline void RandGraph::set_rand_engine(std::mt19937 & e)
+inline void RandGraph::set_rand_engine(std::shared_ptr< std::mt19937 > & e)
 {
 
-    engine = std::make_shared< std::mt19937 >(e);
+    engine = e;
 
 }
 
@@ -47,7 +44,7 @@ inline epiworld_double RandGraph::runif() {
     if (!initialized)
         throw std::logic_error("The object has not been initialized");
 
-    return (*unifd)(engine);
+    return unifd(*engine);
 
 }
 
diff --git a/inst/tinytest/test-lfmcmc.R b/inst/tinytest/test-lfmcmc.R
new file mode 100644
index 00000000..470c4bbf
--- /dev/null
+++ b/inst/tinytest/test-lfmcmc.R
@@ -0,0 +1,96 @@
+# Create Model to use in LFMCMC simulation
+model_seed <- 122
+
+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
+)
+
+# Setup LFMCMC
+## Extract the observed data from the model
+obs_data <- unname(as.integer(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
+  )
+
+  res <- unname(as.integer(get_today_total(model_sir)))
+  return(res)
+}
+
+## 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 <- 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))
+
+  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)
+
+# Run LFMCMC simulation
+## Set initial parameters
+par0 <- as.double(c(0.1, 0.5))
+n_samp <- 2000
+epsil <- as.double(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)
\ No newline at end of file
diff --git a/man/LFMCMC.Rd b/man/LFMCMC.Rd
new file mode 100644
index 00000000..2bf04342
--- /dev/null
+++ b/man/LFMCMC.Rd
@@ -0,0 +1,158 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/LFMCMC.R
+\name{LFMCMC}
+\alias{LFMCMC}
+\alias{epiworld_lfmcmc}
+\alias{run_lfmcmc}
+\alias{set_observed_data}
+\alias{set_proposal_fun}
+\alias{use_proposal_norm_reflective}
+\alias{set_simulation_fun}
+\alias{set_summary_fun}
+\alias{set_kernel_fun}
+\alias{use_kernel_fun_gaussian}
+\alias{set_par_names}
+\alias{set_stats_names}
+\alias{print.epiworld_lfmcmc}
+\title{Likelihood-Free Markhov Chain Monte Carlo (LFMCMC)}
+\usage{
+LFMCMC(model)
+
+run_lfmcmc(lfmcmc, params_init_, n_samples_, epsilon_, seed = NULL)
+
+set_observed_data(lfmcmc, observed_data_)
+
+set_proposal_fun(lfmcmc, fun)
+
+use_proposal_norm_reflective(lfmcmc)
+
+set_simulation_fun(lfmcmc, fun)
+
+set_summary_fun(lfmcmc, fun)
+
+set_kernel_fun(lfmcmc, fun)
+
+use_kernel_fun_gaussian(lfmcmc)
+
+set_par_names(lfmcmc, names)
+
+set_stats_names(lfmcmc, names)
+
+\method{print}{epiworld_lfmcmc}(x, ...)
+}
+\arguments{
+\item{model}{A model of class \link{epiworld_model}}
+
+\item{lfmcmc}{LFMCMC model}
+
+\item{params_init_}{Initial model parameters}
+
+\item{n_samples_}{Number of samples}
+
+\item{epsilon_}{Epsilon parameter}
+
+\item{seed}{Random engine seed}
+
+\item{observed_data_}{Observed data}
+
+\item{fun}{The LFMCMC kernel function}
+
+\item{names}{The model stats names}
+
+\item{x}{LFMCMC model to print}
+
+\item{...}{Ignored}
+}
+\value{
+The \code{LFMCMC} function returns a model of class \link{epiworld_lfmcmc}.
+
+The simulated model of class \link{epiworld_lfmcmc}.
+
+The lfmcmc model with the observed data added
+
+The lfmcmc model with the proposal function added
+
+The LFMCMC model with proposal function set to norm reflective
+
+The lfmcmc model with the simulation function added
+
+The lfmcmc model with the summary function added
+
+The lfmcmc model with the kernel function added
+
+The LFMCMC model with kernel function set to gaussian
+
+The lfmcmc model with the parameter names added
+
+The lfmcmc model with the stats names added
+
+The lfmcmc model
+}
+\description{
+Likelihood-Free Markhov Chain Monte Carlo (LFMCMC)
+}
+\details{
+Performs a Likelihood-Free Markhov Chain Monte Carlo simulation
+}
+\examples{
+## Setup an SIR model to use in the simulation
+model_seed <- 122
+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)
+
+## Setup LFMCMC
+# Extract the observed data from the model
+obs_data <- unname(as.integer(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)))
+  return(res)
+}
+
+# Define the summary function
+sumfun <- function(dat) {
+  return(dat)
+}
+
+# Create the LFMCMC model
+lfmcmc_model <- LFMCMC(model_sir) |>
+  set_simulation_fun(simfun) |>
+  set_summary_fun(sumfun) |>
+  use_proposal_norm_reflective() |>
+  use_kernel_fun_gaussian() |>
+  set_observed_data(obs_data)
+
+## Run LFMCMC simulation
+# Set initial parameters
+par0 <- as.double(c(0.1, 0.5))
+n_samp <- 2000
+epsil <- as.double(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)
+}
diff --git a/src/Makevars.in b/src/Makevars.in
index fc3bb29e..845e17c4 100644
--- a/src/Makevars.in
+++ b/src/Makevars.in
@@ -1,3 +1,3 @@
 PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) @OPENMP_FLAG@
 PKG_CXXFLAGS=@OPENMP_FLAG@ -I../inst/include/ $(EPI_CONFIG) \
-    -Dprintf_epiworld=Rprintf
+    -Dprintf_epiworld=Rprintf -Depiworld_double=double
diff --git a/src/Makevars.win b/src/Makevars.win
index 15c6b642..31bc1854 100644
--- a/src/Makevars.win
+++ b/src/Makevars.win
@@ -2,7 +2,7 @@ PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
 
 # This is necesary since ARMADILLO now supports OpenMP
 PKG_CXXFLAGS=$(SHLIB_OPENMP_CXXFLAGS) -I../inst/include/ \
-    -Dprintf_epiworld=Rprintf
+    -Dprintf_epiworld=Rprintf -Depiworld_double=double
 
 # For testing
 #PKG_CXXFLAGS=-Wall
diff --git a/src/cpp11.cpp b/src/cpp11.cpp
index b8ad0851..a50bdc74 100644
--- a/src/cpp11.cpp
+++ b/src/cpp11.cpp
@@ -390,6 +390,90 @@ extern "C" SEXP _epiworldR_ModelSEIRMixing_cpp(SEXP name, SEXP n, SEXP prevalenc
     return cpp11::as_sexp(ModelSEIRMixing_cpp(cpp11::as_cpp<cpp11::decay_t<std::string>>(name), cpp11::as_cpp<cpp11::decay_t<unsigned int>>(n), cpp11::as_cpp<cpp11::decay_t<double>>(prevalence), cpp11::as_cpp<cpp11::decay_t<double>>(contact_rate), cpp11::as_cpp<cpp11::decay_t<double>>(transmission_rate), cpp11::as_cpp<cpp11::decay_t<double>>(incubation_days), cpp11::as_cpp<cpp11::decay_t<double>>(recovery_rate), cpp11::as_cpp<cpp11::decay_t<std::vector< double >>>(contact_matrix)));
   END_CPP11
 }
+// lfmcmc.cpp
+SEXP LFMCMC_cpp(SEXP model);
+extern "C" SEXP _epiworldR_LFMCMC_cpp(SEXP model) {
+  BEGIN_CPP11
+    return cpp11::as_sexp(LFMCMC_cpp(cpp11::as_cpp<cpp11::decay_t<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);
+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)));
+  END_CPP11
+}
+// lfmcmc.cpp
+SEXP set_observed_data_cpp(SEXP lfmcmc, std::vector< int > 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_)));
+  END_CPP11
+}
+// lfmcmc.cpp
+SEXP set_proposal_fun_cpp(SEXP lfmcmc, cpp11::function fun);
+extern "C" SEXP _epiworldR_set_proposal_fun_cpp(SEXP lfmcmc, SEXP fun) {
+  BEGIN_CPP11
+    return cpp11::as_sexp(set_proposal_fun_cpp(cpp11::as_cpp<cpp11::decay_t<SEXP>>(lfmcmc), cpp11::as_cpp<cpp11::decay_t<cpp11::function>>(fun)));
+  END_CPP11
+}
+// lfmcmc.cpp
+SEXP use_proposal_norm_reflective_cpp(SEXP lfmcmc);
+extern "C" SEXP _epiworldR_use_proposal_norm_reflective_cpp(SEXP lfmcmc) {
+  BEGIN_CPP11
+    return cpp11::as_sexp(use_proposal_norm_reflective_cpp(cpp11::as_cpp<cpp11::decay_t<SEXP>>(lfmcmc)));
+  END_CPP11
+}
+// lfmcmc.cpp
+SEXP set_simulation_fun_cpp(SEXP lfmcmc, cpp11::function fun);
+extern "C" SEXP _epiworldR_set_simulation_fun_cpp(SEXP lfmcmc, SEXP fun) {
+  BEGIN_CPP11
+    return cpp11::as_sexp(set_simulation_fun_cpp(cpp11::as_cpp<cpp11::decay_t<SEXP>>(lfmcmc), cpp11::as_cpp<cpp11::decay_t<cpp11::function>>(fun)));
+  END_CPP11
+}
+// lfmcmc.cpp
+SEXP set_summary_fun_cpp(SEXP lfmcmc, cpp11::function fun);
+extern "C" SEXP _epiworldR_set_summary_fun_cpp(SEXP lfmcmc, SEXP fun) {
+  BEGIN_CPP11
+    return cpp11::as_sexp(set_summary_fun_cpp(cpp11::as_cpp<cpp11::decay_t<SEXP>>(lfmcmc), cpp11::as_cpp<cpp11::decay_t<cpp11::function>>(fun)));
+  END_CPP11
+}
+// lfmcmc.cpp
+SEXP set_kernel_fun_cpp(SEXP lfmcmc, cpp11::function fun);
+extern "C" SEXP _epiworldR_set_kernel_fun_cpp(SEXP lfmcmc, SEXP fun) {
+  BEGIN_CPP11
+    return cpp11::as_sexp(set_kernel_fun_cpp(cpp11::as_cpp<cpp11::decay_t<SEXP>>(lfmcmc), cpp11::as_cpp<cpp11::decay_t<cpp11::function>>(fun)));
+  END_CPP11
+}
+// lfmcmc.cpp
+SEXP use_kernel_fun_gaussian_cpp(SEXP lfmcmc);
+extern "C" SEXP _epiworldR_use_kernel_fun_gaussian_cpp(SEXP lfmcmc) {
+  BEGIN_CPP11
+    return cpp11::as_sexp(use_kernel_fun_gaussian_cpp(cpp11::as_cpp<cpp11::decay_t<SEXP>>(lfmcmc)));
+  END_CPP11
+}
+// lfmcmc.cpp
+SEXP set_par_names_cpp(SEXP lfmcmc, std::vector< std::string > names);
+extern "C" SEXP _epiworldR_set_par_names_cpp(SEXP lfmcmc, SEXP names) {
+  BEGIN_CPP11
+    return cpp11::as_sexp(set_par_names_cpp(cpp11::as_cpp<cpp11::decay_t<SEXP>>(lfmcmc), cpp11::as_cpp<cpp11::decay_t<std::vector< std::string >>>(names)));
+  END_CPP11
+}
+// lfmcmc.cpp
+SEXP set_stats_names_cpp(SEXP lfmcmc, std::vector< std::string > names);
+extern "C" SEXP _epiworldR_set_stats_names_cpp(SEXP lfmcmc, SEXP names) {
+  BEGIN_CPP11
+    return cpp11::as_sexp(set_stats_names_cpp(cpp11::as_cpp<cpp11::decay_t<SEXP>>(lfmcmc), cpp11::as_cpp<cpp11::decay_t<std::vector< std::string >>>(names)));
+  END_CPP11
+}
+// lfmcmc.cpp
+SEXP print_lfmcmc_cpp(SEXP lfmcmc);
+extern "C" SEXP _epiworldR_print_lfmcmc_cpp(SEXP lfmcmc) {
+  BEGIN_CPP11
+    return cpp11::as_sexp(print_lfmcmc_cpp(cpp11::as_cpp<cpp11::decay_t<SEXP>>(lfmcmc)));
+  END_CPP11
+}
 // model.cpp
 SEXP print_cpp(SEXP m, bool lite);
 extern "C" SEXP _epiworldR_print_cpp(SEXP m, SEXP lite) {
@@ -911,6 +995,7 @@ extern "C" SEXP _epiworldR_distribute_virus_to_set_cpp(SEXP agents_ids) {
 
 extern "C" {
 static const R_CallMethodDef CallEntries[] = {
+    {"_epiworldR_LFMCMC_cpp",                           (DL_FUNC) &_epiworldR_LFMCMC_cpp,                            1},
     {"_epiworldR_ModelDiffNet_cpp",                     (DL_FUNC) &_epiworldR_ModelDiffNet_cpp,                      8},
     {"_epiworldR_ModelSEIRCONN_cpp",                    (DL_FUNC) &_epiworldR_ModelSEIRCONN_cpp,                     7},
     {"_epiworldR_ModelSEIRDCONN_cpp",                   (DL_FUNC) &_epiworldR_ModelSEIRDCONN_cpp,                    8},
@@ -990,6 +1075,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_epiworldR_print_cpp",                            (DL_FUNC) &_epiworldR_print_cpp,                             2},
     {"_epiworldR_print_entity_cpp",                     (DL_FUNC) &_epiworldR_print_entity_cpp,                      1},
     {"_epiworldR_print_global_action_cpp",              (DL_FUNC) &_epiworldR_print_global_action_cpp,               1},
+    {"_epiworldR_print_lfmcmc_cpp",                     (DL_FUNC) &_epiworldR_print_lfmcmc_cpp,                      1},
     {"_epiworldR_print_tool_cpp",                       (DL_FUNC) &_epiworldR_print_tool_cpp,                        1},
     {"_epiworldR_print_virus_cpp",                      (DL_FUNC) &_epiworldR_print_virus_cpp,                       1},
     {"_epiworldR_queuing_off_cpp",                      (DL_FUNC) &_epiworldR_queuing_off_cpp,                       1},
@@ -999,6 +1085,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_epiworldR_rm_tool_cpp",                          (DL_FUNC) &_epiworldR_rm_tool_cpp,                           2},
     {"_epiworldR_rm_virus_cpp",                         (DL_FUNC) &_epiworldR_rm_virus_cpp,                          2},
     {"_epiworldR_run_cpp",                              (DL_FUNC) &_epiworldR_run_cpp,                               3},
+    {"_epiworldR_run_lfmcmc_cpp",                       (DL_FUNC) &_epiworldR_run_lfmcmc_cpp,                        5},
     {"_epiworldR_run_multiple_cpp",                     (DL_FUNC) &_epiworldR_run_multiple_cpp,                      8},
     {"_epiworldR_set_agents_data_cpp",                  (DL_FUNC) &_epiworldR_set_agents_data_cpp,                   3},
     {"_epiworldR_set_death_reduction_cpp",              (DL_FUNC) &_epiworldR_set_death_reduction_cpp,               2},
@@ -1010,9 +1097,12 @@ static const R_CallMethodDef CallEntries[] = {
     {"_epiworldR_set_incubation_cpp",                   (DL_FUNC) &_epiworldR_set_incubation_cpp,                    2},
     {"_epiworldR_set_incubation_fun_cpp",               (DL_FUNC) &_epiworldR_set_incubation_fun_cpp,                3},
     {"_epiworldR_set_incubation_ptr_cpp",               (DL_FUNC) &_epiworldR_set_incubation_ptr_cpp,                3},
+    {"_epiworldR_set_kernel_fun_cpp",                   (DL_FUNC) &_epiworldR_set_kernel_fun_cpp,                    2},
     {"_epiworldR_set_name_cpp",                         (DL_FUNC) &_epiworldR_set_name_cpp,                          2},
     {"_epiworldR_set_name_tool_cpp",                    (DL_FUNC) &_epiworldR_set_name_tool_cpp,                     2},
     {"_epiworldR_set_name_virus_cpp",                   (DL_FUNC) &_epiworldR_set_name_virus_cpp,                    2},
+    {"_epiworldR_set_observed_data_cpp",                (DL_FUNC) &_epiworldR_set_observed_data_cpp,                 2},
+    {"_epiworldR_set_par_names_cpp",                    (DL_FUNC) &_epiworldR_set_par_names_cpp,                     2},
     {"_epiworldR_set_param_cpp",                        (DL_FUNC) &_epiworldR_set_param_cpp,                         3},
     {"_epiworldR_set_prob_death_cpp",                   (DL_FUNC) &_epiworldR_set_prob_death_cpp,                    2},
     {"_epiworldR_set_prob_death_fun_cpp",               (DL_FUNC) &_epiworldR_set_prob_death_fun_cpp,                3},
@@ -1023,9 +1113,13 @@ static const R_CallMethodDef CallEntries[] = {
     {"_epiworldR_set_prob_recovery_cpp",                (DL_FUNC) &_epiworldR_set_prob_recovery_cpp,                 2},
     {"_epiworldR_set_prob_recovery_fun_cpp",            (DL_FUNC) &_epiworldR_set_prob_recovery_fun_cpp,             3},
     {"_epiworldR_set_prob_recovery_ptr_cpp",            (DL_FUNC) &_epiworldR_set_prob_recovery_ptr_cpp,             3},
+    {"_epiworldR_set_proposal_fun_cpp",                 (DL_FUNC) &_epiworldR_set_proposal_fun_cpp,                  2},
     {"_epiworldR_set_recovery_enhancer_cpp",            (DL_FUNC) &_epiworldR_set_recovery_enhancer_cpp,             2},
     {"_epiworldR_set_recovery_enhancer_fun_cpp",        (DL_FUNC) &_epiworldR_set_recovery_enhancer_fun_cpp,         3},
     {"_epiworldR_set_recovery_enhancer_ptr_cpp",        (DL_FUNC) &_epiworldR_set_recovery_enhancer_ptr_cpp,         3},
+    {"_epiworldR_set_simulation_fun_cpp",               (DL_FUNC) &_epiworldR_set_simulation_fun_cpp,                2},
+    {"_epiworldR_set_stats_names_cpp",                  (DL_FUNC) &_epiworldR_set_stats_names_cpp,                   2},
+    {"_epiworldR_set_summary_fun_cpp",                  (DL_FUNC) &_epiworldR_set_summary_fun_cpp,                   2},
     {"_epiworldR_set_susceptibility_reduction_cpp",     (DL_FUNC) &_epiworldR_set_susceptibility_reduction_cpp,      2},
     {"_epiworldR_set_susceptibility_reduction_fun_cpp", (DL_FUNC) &_epiworldR_set_susceptibility_reduction_fun_cpp,  3},
     {"_epiworldR_set_susceptibility_reduction_ptr_cpp", (DL_FUNC) &_epiworldR_set_susceptibility_reduction_ptr_cpp,  3},
@@ -1035,6 +1129,8 @@ static const R_CallMethodDef CallEntries[] = {
     {"_epiworldR_size_cpp",                             (DL_FUNC) &_epiworldR_size_cpp,                              1},
     {"_epiworldR_tool_cpp",                             (DL_FUNC) &_epiworldR_tool_cpp,                              7},
     {"_epiworldR_tool_fun_logit_cpp",                   (DL_FUNC) &_epiworldR_tool_fun_logit_cpp,                    3},
+    {"_epiworldR_use_kernel_fun_gaussian_cpp",          (DL_FUNC) &_epiworldR_use_kernel_fun_gaussian_cpp,           1},
+    {"_epiworldR_use_proposal_norm_reflective_cpp",     (DL_FUNC) &_epiworldR_use_proposal_norm_reflective_cpp,      1},
     {"_epiworldR_verbose_off_cpp",                      (DL_FUNC) &_epiworldR_verbose_off_cpp,                       1},
     {"_epiworldR_verbose_on_cpp",                       (DL_FUNC) &_epiworldR_verbose_on_cpp,                        1},
     {"_epiworldR_virus_cpp",                            (DL_FUNC) &_epiworldR_virus_cpp,                             8},
diff --git a/src/lfmcmc.cpp b/src/lfmcmc.cpp
new file mode 100644
index 00000000..075d4017
--- /dev/null
+++ b/src/lfmcmc.cpp
@@ -0,0 +1,237 @@
+#include "cpp11.hpp"
+#include "cpp11/external_pointer.hpp"
+#include "cpp11/r_vector.hpp"
+#include "cpp11/sexp.hpp"
+#include <iostream>
+
+#include "epiworld-common.h"
+
+using namespace epiworld;
+
+#define TData_default std::vector< int >
+
+#define WrapLFMCMC(a) \
+    cpp11::external_pointer<LFMCMC<TData_default>> (a)
+
+// LFMCMC definitions:
+// https://github.com/UofUEpiBio/epiworld/tree/master/include/epiworld/math/lfmcmc
+
+// *************************************
+// LFMCMC Function
+// *************************************
+
+[[cpp11::register]]
+SEXP LFMCMC_cpp(
+    SEXP model
+) {
+    WrapLFMCMC(lfmcmc_ptr)(
+        new LFMCMC<TData_default>()
+    );
+
+    cpp11::external_pointer<Model<int>> modelptr(model);
+    lfmcmc_ptr->set_rand_engine(modelptr->get_rand_endgine());
+
+    return lfmcmc_ptr;
+}
+
+// *************************************
+// LFMCMC Run Function
+// *************************************
+
+[[cpp11::register]]
+SEXP run_lfmcmc_cpp(
+    SEXP lfmcmc,
+    std::vector<epiworld_double> params_init_,
+    size_t n_samples_,
+    epiworld_double epsilon_,
+    int seed
+) {
+    WrapLFMCMC(lfmcmc_ptr)(lfmcmc);
+    lfmcmc_ptr->run(params_init_, n_samples_, epsilon_, seed);
+    return lfmcmc;
+}
+
+// *************************************
+// LFMCMC Setup Functions
+// *************************************
+
+// observed_data_ should be of type TData
+[[cpp11::register]]
+SEXP set_observed_data_cpp(
+    SEXP lfmcmc,
+    std::vector< int > observed_data_
+) {
+    WrapLFMCMC(lfmcmc_ptr)(lfmcmc);
+    lfmcmc_ptr->set_observed_data(observed_data_);
+    return lfmcmc;
+}
+
+[[cpp11::register]]
+SEXP set_proposal_fun_cpp(
+    SEXP lfmcmc,
+    cpp11::function fun
+) {
+
+    LFMCMCProposalFun<TData_default> fun_call = [fun](
+        std::vector< epiworld_double >& params_now,
+        const std::vector< epiworld_double >& params_prev,
+        LFMCMC<TData_default>*
+        ) -> void {
+
+        auto params_doubles = cpp11::doubles(params_prev);
+
+        auto res_tmp = cpp11::doubles(fun(params_doubles));
+
+        std::copy(
+            res_tmp.begin(),
+            res_tmp.end(),
+            params_now.begin()
+            );
+
+        return;
+    };
+
+    WrapLFMCMC(lfmcmc_ptr)(lfmcmc);
+
+    lfmcmc_ptr->set_proposal_fun(fun_call);
+
+    return lfmcmc;
+}
+
+// Use proposal function defined in epiworld
+[[cpp11::register]]
+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));
+    return lfmcmc;
+}
+
+[[cpp11::register]]
+SEXP set_simulation_fun_cpp(
+    SEXP lfmcmc,
+    cpp11::function fun
+) {
+
+    LFMCMCSimFun<TData_default> fun_call = [fun](
+        const std::vector<epiworld_double>& params,
+        LFMCMC<TData_default>*
+        ) -> TData_default {
+
+        auto params_doubles = cpp11::doubles(params);
+
+        return cpp11::as_cpp<TData_default>(
+            cpp11::integers(fun(params_doubles))
+            );
+    };
+
+    WrapLFMCMC(lfmcmc_ptr)(lfmcmc);
+
+    lfmcmc_ptr->set_simulation_fun(fun_call);
+
+    return lfmcmc;
+}
+
+[[cpp11::register]]
+SEXP set_summary_fun_cpp(
+    SEXP lfmcmc,
+    cpp11::function fun
+) {
+
+    LFMCMCSummaryFun<TData_default> fun_call = [fun](
+        std::vector< epiworld_double >& res,
+        const TData_default& dat,
+        LFMCMC<TData_default>*
+        ) -> void {
+
+        if (res.size() == 0u)
+            res.resize(dat.size());
+
+        auto dat_int = cpp11::integers(dat);
+        auto res_tmp = cpp11::integers(fun(dat_int));
+
+        std::copy(res_tmp.begin(), res_tmp.end(), res.begin());
+
+        return;
+    };
+
+    WrapLFMCMC(lfmcmc_ptr)(lfmcmc);
+
+    lfmcmc_ptr->set_summary_fun(fun_call);
+
+    return lfmcmc;
+}
+
+[[cpp11::register]]
+SEXP set_kernel_fun_cpp(
+    SEXP lfmcmc,
+    cpp11::function fun
+) {
+
+    LFMCMCKernelFun<TData_default> fun_call = [fun](
+        const std::vector< epiworld_double >& stats_now,
+        const std::vector< epiworld_double >& stats_obs,
+        epiworld_double epsilon,
+        LFMCMC<TData_default>*
+        ) -> epiworld_double {
+
+        auto stats_now_doubles = cpp11::doubles(stats_now);
+        auto stats_obs_doubles = cpp11::doubles(stats_obs);
+
+        return cpp11::as_cpp<epiworld_double>(
+            fun(stats_now_doubles, stats_obs_doubles, epsilon)
+            );
+    };
+
+    WrapLFMCMC(lfmcmc_ptr)(lfmcmc);
+
+    lfmcmc_ptr->set_kernel_fun(fun_call);
+
+    return lfmcmc;
+}
+
+// Use kernel function defined in epiworld
+[[cpp11::register]]
+SEXP use_kernel_fun_gaussian_cpp(
+    SEXP lfmcmc
+) {
+    WrapLFMCMC(lfmcmc_ptr)(lfmcmc);
+    lfmcmc_ptr->set_kernel_fun(kernel_fun_gaussian<std::vector<int>>);
+    return lfmcmc;
+}
+
+// *************************************
+// LFMCMC Printing Functions
+// *************************************
+
+[[cpp11::register]]
+SEXP set_par_names_cpp(
+    SEXP lfmcmc,
+    std::vector< std::string > names
+) {
+    WrapLFMCMC(lfmcmc_ptr)(lfmcmc);
+    lfmcmc_ptr->set_par_names(names);
+    return lfmcmc;
+}
+
+[[cpp11::register]]
+SEXP set_stats_names_cpp(
+    SEXP lfmcmc,
+    std::vector< std::string > names
+) {
+    WrapLFMCMC(lfmcmc_ptr)(lfmcmc);
+    lfmcmc_ptr->set_stats_names(names);
+    return lfmcmc;
+}
+
+[[cpp11::register]]
+SEXP print_lfmcmc_cpp(
+    SEXP lfmcmc
+) {
+    WrapLFMCMC(lfmcmc_ptr)(lfmcmc);
+    lfmcmc_ptr->print();
+    return lfmcmc;
+}
+
+#undef WrapLFMCMC
diff --git a/vignettes/likelihood-free-mcmc.Rmd b/vignettes/likelihood-free-mcmc.Rmd
new file mode 100644
index 00000000..c5e87fe4
--- /dev/null
+++ b/vignettes/likelihood-free-mcmc.Rmd
@@ -0,0 +1,131 @@
+---
+title: "Likelihood Free Markhov Chain Monte Carlo (LFMCMC)"
+author:
+  - Andrew Pulsipher
+date: "`r Sys.Date()`"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{LFMCMC}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+```{r setup, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>", out.width = "80%", fig.width = 7, fig.height = 5,
+  fig.align = "center"
+)
+```
+# 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.
+```{r sir-setup}
+library(epiworldR)
+
+model_seed <- 122
+
+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
+)
+
+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)))
+
+# 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
+  )
+
+  res <- unname(as.integer(get_today_total(model_sir)))
+  return(res)
+}
+
+# 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 <- 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))
+
+  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)
+```
+
+## Run LFMCMC simulation
+```{r lfmcmc-run}
+# Set initial parameters
+par0 <- as.double(c(0.1, 0.5))
+n_samp <- 2000
+epsil <- as.double(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)
+```