diff --git a/.gitignore b/.gitignore
index 65d4fd1cc..b9c70bc28 100644
--- a/.gitignore
+++ b/.gitignore
@@ -25,3 +25,7 @@ data/models/*reparam
 docs
 /doc/
 /Meta/
+vignettes/**_cache/
+*.pdf
+.vscode/
+vignettes/figures/
diff --git a/DESCRIPTION b/DESCRIPTION
index d0b4f6e05..c7162d733 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -34,7 +34,8 @@ Imports:
     rstan (>= 2.26.0),
     dplyr,
     tibble,
-    lubridate
+    lubridate,
+    primarycensored
 Suggests:
     bookdown,
     testthat (>= 3.0.0),
diff --git a/NAMESPACE b/NAMESPACE
index 8fb43ebb5..5fa0b68b8 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -38,6 +38,8 @@ export(epidist_family_prior)
 export(epidist_family_reparam)
 export(epidist_formula)
 export(epidist_formula_model)
+export(epidist_gen_posterior_epred)
+export(epidist_gen_posterior_predict)
 export(epidist_model_prior)
 export(epidist_prior)
 export(epidist_stancode)
diff --git a/R/gen.R b/R/gen.R
new file mode 100644
index 000000000..ebef86b3b
--- /dev/null
+++ b/R/gen.R
@@ -0,0 +1,73 @@
+#' Create a function to draw from the posterior predictive distribution for a
+#' latent model
+#'
+#' This function creates a function that draws from the posterior predictive
+#' distribution for a latent model using [primarycensored::rpcens()] to handle
+#' censoring and truncation. The returned function takes a `prep` argument from
+#' `brms` and returns posterior predictions. This is used internally by
+#' [brms::posterior_predict()] to generate predictions for latent models.
+#'
+#' @inheritParams epidist_family
+#'
+#' @return A function that takes a `prep` argument from brms and returns a
+#' matrix of posterior predictions, with one row per posterior draw and one
+#' column per observation. The `prep` object must have the following variables:
+#' * `vreal1`: relative observation time
+#' * `vreal2`: primary event window
+#' * `vreal3`: secondary event window
+#'
+#' @seealso [brms::posterior_predict()] for details on how this is used within
+#' `brms`, [primarycensored::rpcens()] for details on the censoring approach
+#' @autoglobal
+#' @family gen
+#' @export
+epidist_gen_posterior_predict <- function(family) {
+  dist_fn <- .get_brms_fn("posterior_predict", family)
+
+  rdist <- function(n, i, prep, ...) {
+    prep$ndraws <- n
+    do.call(dist_fn, list(i = i, prep = prep))
+  }
+
+  .predict <- function(i, prep, ...) {
+    relative_obs_time <- prep$data$vreal1[i]
+    pwindow <- prep$data$vreal2[i]
+    swindow <- prep$data$vreal3[i]
+
+    as.matrix(primarycensored::rpcens(
+      n = prep$ndraws,
+      rdist = rdist,
+      rprimary = stats::runif,
+      pwindow = prep$data$vreal2[i],
+      swindow = prep$data$vreal3[i],
+      D = prep$data$vreal1[i],
+      i = i,
+      prep = prep
+    ))
+  }
+  return(.predict)
+}
+
+#' Create a function to draw from the expected value of the posterior predictive
+#' distribution for a latent model
+#'
+#' This function creates a function that calculates the expected value of the
+#' posterior predictive distribution for a latent model. The returned function
+#' takes a `prep` argument (from brms) and returns posterior expected values.
+#' This is used internally by [brms::posterior_epred()] to calculate expected
+#' values for latent models.
+#'
+#' @inheritParams epidist_family
+#'
+#' @return A function that takes a prep argument from brms and returns a matrix
+#' of posterior expected values, with one row per posterior draw and one column
+#' per observation.
+#'
+#' @seealso [brms::posterior_epred()] for details on how this is used within
+#' `brms`.
+#' @autoglobal
+#' @family gen
+#' @export
+epidist_gen_posterior_epred <- function(family) {
+  .get_brms_fn("posterior_epred", family)
+}
diff --git a/R/globals.R b/R/globals.R
index 41593e849..bfa8d919d 100644
--- a/R/globals.R
+++ b/R/globals.R
@@ -1,16 +1,9 @@
 # Generated by roxyglobals: do not edit by hand
 
 utils::globalVariables(c(
-  ".data", # <epidist_diagnostics>
   "samples", # <epidist_diagnostics>
-  ".data", # <as_epidist_latent_model.epidist_linelist_data>
   "woverlap", # <epidist_stancode.epidist_latent_model>
-  ".data", # <as_epidist_naive_model.epidist_linelist_data>
-  ".data", # <add_mean_sd.lognormal_samples>
-  ".data", # <add_mean_sd.gamma_samples>
   "rlnorm", # <simulate_secondary>
-  ".data", # <simulate_secondary>
-  ".data", # <.replace_prior>
   "prior_new", # <.replace_prior>
   "source_new", # <.replace_prior>
   NULL
diff --git a/R/latent_gamma.R b/R/latent_gamma.R
deleted file mode 100644
index 46b7c2d87..000000000
--- a/R/latent_gamma.R
+++ /dev/null
@@ -1,84 +0,0 @@
-#' Draws from the posterior predictive distribution of the `latent_gamma` family
-#'
-#' See [brms::posterior_predict()].
-#'
-#' @param i The index of the observation to predict
-#' @param prep The result of a call to [brms::posterior_predict()]
-#' @param ... Additional arguments
-#' @autoglobal
-#' @keywords internal
-posterior_predict_latent_gamma <- function(i, prep, ...) { # nolint: object_length_linter
-  mu <- brms::get_dpar(prep, "mu", i = i)
-  shape <- brms::get_dpar(prep, "shape", i = i)
-
-  relative_obs_time <- prep$data$vreal1[i]
-  pwindow <- prep$data$vreal2[i]
-  swindow <- prep$data$vreal3[i]
-
-  .predict <- function(s) {
-    d_censored <- relative_obs_time + 1
-    # while loop to impose the truncation
-    while (d_censored > relative_obs_time) {
-      p_latent <- stats::runif(1, 0, 1) * pwindow
-      d_latent <- stats::rgamma(1, shape = shape[s], scale = mu[s] / shape[s])
-      s_latent <- p_latent + d_latent
-      p_censored <- .floor_mult(p_latent, pwindow)
-      s_censored <- .floor_mult(s_latent, swindow)
-      d_censored <- s_censored - p_censored
-    }
-    return(d_censored)
-  }
-
-  # Within brms this is a helper function called rblapply
-  do.call(rbind, lapply(seq_len(prep$ndraws), .predict))
-}
-
-#' Draws from the expected value of the posterior predictive distribution of the
-#' `latent_gamma` family
-#'
-#' See [brms::posterior_epred()].
-#'
-#' @param prep The result of a call to [`brms::prepare_predictions`]
-#' @autoglobal
-#' @keywords internal
-posterior_epred_latent_gamma <- function(prep) { # nolint: object_length_linter
-  brms::get_dpar(prep, "mu")
-}
-
-#' Calculate the pointwise log likelihood of the `latent_gamma` family
-#'
-#' See [brms::log_lik()].
-#'
-#' @param i The index of the observation to calculate the log likelihood of
-#' @param prep The result of a call to [brms::prepare_predictions()]
-#' @autoglobal
-#' @keywords internal
-log_lik_latent_gamma <- function(i, prep) {
-  mu <- brms::get_dpar(prep, "mu", i = i)
-  shape <- brms::get_dpar(prep, "shape", i = i)
-  y <- prep$data$Y[i]
-  relative_obs_time <- prep$data$vreal1[i]
-  pwindow <- prep$data$vreal2[i]
-  swindow <- prep$data$vreal3[i]
-
-  swindow_raw <- stats::runif(prep$ndraws)
-  pwindow_raw <- stats::runif(prep$ndraws)
-
-  swindow <- swindow_raw * swindow
-
-  # For no overlap calculate as usual, for overlap ensure pwindow < swindow
-  if (i %in% prep$data$noverlap) {
-    pwindow <- pwindow_raw * pwindow
-  } else {
-    pwindow <- pwindow_raw * swindow
-  }
-
-  d <- y - pwindow + swindow
-  obs_time <- relative_obs_time - pwindow
-  lpdf <- stats::dgamma(d, shape = shape, scale = mu / shape, log = TRUE)
-  lcdf <- stats::pgamma(
-    obs_time,
-    shape = shape, scale = mu / shape, log.p = TRUE
-  )
-  return(lpdf - lcdf)
-}
diff --git a/R/latent_lognormal.R b/R/latent_lognormal.R
deleted file mode 100644
index a9de72c26..000000000
--- a/R/latent_lognormal.R
+++ /dev/null
@@ -1,87 +0,0 @@
-#' Draws from the posterior predictive distribution of the `latent_lognormal`
-#' family
-#'
-#' See [brms::posterior_predict()].
-#'
-#' @param i The index of the observation to predict
-#' @param prep The result of a call to [brms::posterior_predict()]
-#' @param ... Additional arguments
-#' @autoglobal
-#' @keywords internal
-posterior_predict_latent_lognormal <- function(i, prep, ...) { # nolint: object_length_linter
-  mu <- brms::get_dpar(prep, "mu", i = i)
-  sigma <- brms::get_dpar(prep, "sigma", i = i)
-
-  relative_obs_time <- prep$data$vreal1[i]
-  pwindow <- prep$data$vreal2[i]
-  swindow <- prep$data$vreal3[i]
-
-  .predict <- function(s) {
-    d_censored <- relative_obs_time + 1
-    # while loop to impose the truncation
-    while (d_censored > relative_obs_time) {
-      p_latent <- stats::runif(1, 0, 1) * pwindow
-      d_latent <- stats::rlnorm(1, meanlog = mu[s], sdlog = sigma[s])
-      s_latent <- p_latent + d_latent
-      p_censored <- .floor_mult(p_latent, pwindow)
-      s_censored <- .floor_mult(s_latent, swindow)
-      d_censored <- s_censored - p_censored
-    }
-    return(d_censored)
-  }
-
-  # Within brms this is a helper function called rblapply
-  do.call(rbind, lapply(seq_len(prep$ndraws), .predict))
-}
-
-#' Draws from the expected value of the posterior predictive distribution of the
-#' `latent_gamma` family
-#'
-#' See [brms::posterior_epred()].
-#'
-#' @param prep The result of a call to [`brms::prepare_predictions`]
-#' @autoglobal
-#' @keywords internal
-posterior_epred_latent_lognormal <- function(prep) { # nolint: object_length_linter
-  mu <- brms::get_dpar(prep, "mu")
-  sigma <- brms::get_dpar(prep, "sigma")
-  exp(mu + sigma^2 / 2)
-}
-
-#' Calculate the pointwise log likelihood of the `latent_gamma` family
-#'
-#' See [brms::log_lik()].
-#'
-#' @param i The index of the observation to calculate the log likelihood of
-#' @param prep The result of a call to [brms::prepare_predictions()]
-#' @autoglobal
-#' @keywords internal
-log_lik_latent_lognormal <- function(i, prep) {
-  mu <- brms::get_dpar(prep, "mu", i = i)
-  sigma <- brms::get_dpar(prep, "sigma", i = i)
-  y <- prep$data$Y[i]
-  relative_obs_time <- prep$data$vreal1[i]
-  pwindow <- prep$data$vreal2[i]
-  swindow <- prep$data$vreal3[i]
-
-  # Generates values of the swindow_raw and pwindow_raw, but really these should
-  # be extracted from prep or the fitted raws somehow. See:
-  # https://github.com/epinowcast/epidist/issues/267
-  swindow_raw <- stats::runif(prep$ndraws)
-  pwindow_raw <- stats::runif(prep$ndraws)
-
-  swindow <- swindow_raw * swindow
-
-  # For no overlap calculate as usual, for overlap ensure pwindow < swindow
-  if (i %in% prep$data$noverlap) {
-    pwindow <- pwindow_raw * pwindow
-  } else {
-    pwindow <- pwindow_raw * swindow
-  }
-
-  d <- y - pwindow + swindow
-  obs_time <- relative_obs_time - pwindow
-  lpdf <- stats::dlnorm(d, meanlog = mu, sdlog = sigma, log = TRUE)
-  lcdf <- stats::plnorm(obs_time, meanlog = mu, sdlog = sigma, log.p = TRUE)
-  return(lpdf - lcdf)
-}
diff --git a/R/latent_model.R b/R/latent_model.R
index 99ccdd5c9..96215b971 100644
--- a/R/latent_model.R
+++ b/R/latent_model.R
@@ -89,12 +89,76 @@ epidist_family_model.epidist_latent_model <- function(
     ub = c(NA, as.numeric(lapply(family$other_bounds, "[[", "ub"))),
     type = family$type,
     vars = c("pwindow", "swindow", "vreal1"),
-    loop = FALSE
+    loop = FALSE,
+    log_lik = epidist_gen_log_lik_latent(family),
+    posterior_predict = epidist_gen_posterior_predict(family),
+    posterior_epred = epidist_gen_posterior_epred(family)
   )
   custom_family$reparm <- family$reparm
   return(custom_family)
 }
 
+#' Create a function to calculate the pointwise log likelihood of the latent
+#' model
+#'
+#' This function creates a log likelihood function that accounts for the latent
+#' variables in the model, including primary and secondary event windows and
+#' their overlap. The returned function calculates the log likelihood for a
+#' single observation by augmenting the data with the latent variables and
+#' using the underlying brms log likelihood function.
+#'
+#' @seealso [brms::log_lik()] for details on the brms log likelihood interface.
+#'
+#' @inheritParams epidist_family
+#'
+#' @return A function that calculates the log likelihood for a single
+#' observation. The prep object must have the following variables:
+#' * `vreal1`: relative observation time
+#' * `vreal2`: primary event window
+#' * `vreal3`: secondary event window
+#'
+#' @family latent_model
+#' @autoglobal
+epidist_gen_log_lik_latent <- function(family) {
+  # Get internal brms log_lik function
+  log_lik_brms <- .get_brms_fn("log_lik", family)
+
+  .log_lik <- function(i, prep) {
+    y <- prep$data$Y[i]
+    relative_obs_time <- prep$data$vreal1[i]
+    pwindow <- prep$data$vreal2[i]
+    swindow <- prep$data$vreal3[i]
+
+    # Generates values of the swindow_raw and pwindow_raw, but really these
+    # should be extracted from prep or the fitted raws somehow. See:
+    # https://github.com/epinowcast/epidist/issues/267
+    swindow_raw <- stats::runif(prep$ndraws)
+    pwindow_raw <- stats::runif(prep$ndraws)
+
+    swindow <- swindow_raw * swindow
+
+    # For no overlap calculate as usual, for overlap ensure pwindow < swindow
+    if (i %in% prep$data$noverlap) {
+      pwindow <- pwindow_raw * pwindow
+    } else {
+      pwindow <- pwindow_raw * swindow
+    }
+
+    d <- y - pwindow + swindow
+    obs_time <- relative_obs_time - pwindow
+    # Create brms truncation upper bound
+    prep$data$ub <- rep(obs_time, length(prep$data$Y))
+    # Update augmented data
+    prep$data$Y <- rep(d, length(prep$data$Y))
+
+    # Call internal brms log_lik function with augmented data
+    lpdf <- log_lik_brms(i, prep)
+    return(lpdf)
+  }
+
+  return(.log_lik)
+}
+
 #' Define the model-specific component of an `epidist` custom formula
 #'
 #' @inheritParams epidist_formula_model
diff --git a/R/utils.R b/R/utils.R
index 00629dab5..9fd6a419e 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -157,3 +157,22 @@
 
   return(df)
 }
+
+#' Get a brms function by prefix and family
+#'
+#' Helper function to get internal brms functions by constructing their name
+#' from a prefix and family. Used to get functions like `log_lik_*`,
+#' `posterior_predict_*` etc.
+#'
+#' @param prefix Character string prefix of the brms function to get (e.g.
+#' "log_lik")
+#'
+#' @inheritParams epidist_family
+#' @return The requested brms function
+#' @keywords internal
+.get_brms_fn <- function(prefix, family) {
+  get(
+    paste0(prefix, "_", tolower(family$family)),
+    asNamespace("brms")
+  )
+}
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 0f393849a..51ef0149f 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -72,6 +72,10 @@ reference:
   desc: Functions for specifying custom Stan code to put into `brms`
   contents:
   - has_concept("stan")
+- title: Generator model functions
+  desc: Generator functions for creating family specific model utilities.
+  contents:
+  - has_concept("gen")
 - title: Simulation
   desc: Tools for simulating datasets
   contents:
diff --git a/man/as_epidist_latent_model.Rd b/man/as_epidist_latent_model.Rd
index fb2329eda..2aa6e600b 100644
--- a/man/as_epidist_latent_model.Rd
+++ b/man/as_epidist_latent_model.Rd
@@ -17,6 +17,7 @@ Other latent_model:
 \code{\link{as_epidist_latent_model.epidist_linelist_data}()},
 \code{\link{epidist_family_model.epidist_latent_model}()},
 \code{\link{epidist_formula_model.epidist_latent_model}()},
+\code{\link{epidist_gen_log_lik_latent}()},
 \code{\link{is_epidist_latent_model}()},
 \code{\link{new_epidist_latent_model}()}
 }
diff --git a/man/as_epidist_latent_model.epidist_linelist_data.Rd b/man/as_epidist_latent_model.epidist_linelist_data.Rd
index ca3c6acf0..03e839d66 100644
--- a/man/as_epidist_latent_model.epidist_linelist_data.Rd
+++ b/man/as_epidist_latent_model.epidist_linelist_data.Rd
@@ -17,6 +17,7 @@ Other latent_model:
 \code{\link{as_epidist_latent_model}()},
 \code{\link{epidist_family_model.epidist_latent_model}()},
 \code{\link{epidist_formula_model.epidist_latent_model}()},
+\code{\link{epidist_gen_log_lik_latent}()},
 \code{\link{is_epidist_latent_model}()},
 \code{\link{new_epidist_latent_model}()}
 }
diff --git a/man/dot-get_brms_fn.Rd b/man/dot-get_brms_fn.Rd
new file mode 100644
index 000000000..414c7a6e0
--- /dev/null
+++ b/man/dot-get_brms_fn.Rd
@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{.get_brms_fn}
+\alias{.get_brms_fn}
+\title{Get a brms function by prefix and family}
+\usage{
+.get_brms_fn(prefix, family)
+}
+\arguments{
+\item{prefix}{Character string prefix of the brms function to get (e.g.
+"log_lik")}
+
+\item{family}{A description of the response distribution and link function to
+be used in the model. Every family function has a link argument allowing
+users to specify the link function to be applied on the response variable.
+If not specified, default links are used. For details of all supported
+families see \code{\link[=brmsfamily]{brmsfamily()}}. Commonly used, such as \code{\link[=lognormal]{lognormal()}}, are also
+reexported as part of \code{epidist}.}
+}
+\value{
+The requested brms function
+}
+\description{
+Helper function to get internal brms functions by constructing their name
+from a prefix and family. Used to get functions like \verb{log_lik_*},
+\verb{posterior_predict_*} etc.
+}
+\keyword{internal}
diff --git a/man/epidist_family_model.epidist_latent_model.Rd b/man/epidist_family_model.epidist_latent_model.Rd
index ae3dae26b..6794976a2 100644
--- a/man/epidist_family_model.epidist_latent_model.Rd
+++ b/man/epidist_family_model.epidist_latent_model.Rd
@@ -22,6 +22,7 @@ Other latent_model:
 \code{\link{as_epidist_latent_model}()},
 \code{\link{as_epidist_latent_model.epidist_linelist_data}()},
 \code{\link{epidist_formula_model.epidist_latent_model}()},
+\code{\link{epidist_gen_log_lik_latent}()},
 \code{\link{is_epidist_latent_model}()},
 \code{\link{new_epidist_latent_model}()}
 }
diff --git a/man/epidist_formula_model.epidist_latent_model.Rd b/man/epidist_formula_model.epidist_latent_model.Rd
index 5811233df..6a5cdfc46 100644
--- a/man/epidist_formula_model.epidist_latent_model.Rd
+++ b/man/epidist_formula_model.epidist_latent_model.Rd
@@ -25,6 +25,7 @@ Other latent_model:
 \code{\link{as_epidist_latent_model}()},
 \code{\link{as_epidist_latent_model.epidist_linelist_data}()},
 \code{\link{epidist_family_model.epidist_latent_model}()},
+\code{\link{epidist_gen_log_lik_latent}()},
 \code{\link{is_epidist_latent_model}()},
 \code{\link{new_epidist_latent_model}()}
 }
diff --git a/man/epidist_gen_log_lik_latent.Rd b/man/epidist_gen_log_lik_latent.Rd
new file mode 100644
index 000000000..5fc1e1dee
--- /dev/null
+++ b/man/epidist_gen_log_lik_latent.Rd
@@ -0,0 +1,45 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/latent_model.R
+\name{epidist_gen_log_lik_latent}
+\alias{epidist_gen_log_lik_latent}
+\title{Create a function to calculate the pointwise log likelihood of the latent
+model}
+\usage{
+epidist_gen_log_lik_latent(family)
+}
+\arguments{
+\item{family}{A description of the response distribution and link function to
+be used in the model. Every family function has a link argument allowing
+users to specify the link function to be applied on the response variable.
+If not specified, default links are used. For details of all supported
+families see \code{\link[=brmsfamily]{brmsfamily()}}. Commonly used, such as \code{\link[=lognormal]{lognormal()}}, are also
+reexported as part of \code{epidist}.}
+}
+\value{
+A function that calculates the log likelihood for a single
+observation. The prep object must have the following variables:
+\itemize{
+\item \code{vreal1}: relative observation time
+\item \code{vreal2}: primary event window
+\item \code{vreal3}: secondary event window
+}
+}
+\description{
+This function creates a log likelihood function that accounts for the latent
+variables in the model, including primary and secondary event windows and
+their overlap. The returned function calculates the log likelihood for a
+single observation by augmenting the data with the latent variables and
+using the underlying brms log likelihood function.
+}
+\seealso{
+\code{\link[brms:log_lik.brmsfit]{brms::log_lik()}} for details on the brms log likelihood interface.
+
+Other latent_model: 
+\code{\link{as_epidist_latent_model}()},
+\code{\link{as_epidist_latent_model.epidist_linelist_data}()},
+\code{\link{epidist_family_model.epidist_latent_model}()},
+\code{\link{epidist_formula_model.epidist_latent_model}()},
+\code{\link{is_epidist_latent_model}()},
+\code{\link{new_epidist_latent_model}()}
+}
+\concept{latent_model}
diff --git a/man/epidist_gen_posterior_epred.Rd b/man/epidist_gen_posterior_epred.Rd
new file mode 100644
index 000000000..05d072457
--- /dev/null
+++ b/man/epidist_gen_posterior_epred.Rd
@@ -0,0 +1,37 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/gen.R
+\name{epidist_gen_posterior_epred}
+\alias{epidist_gen_posterior_epred}
+\title{Create a function to draw from the expected value of the posterior predictive
+distribution for a latent model}
+\usage{
+epidist_gen_posterior_epred(family)
+}
+\arguments{
+\item{family}{A description of the response distribution and link function to
+be used in the model. Every family function has a link argument allowing
+users to specify the link function to be applied on the response variable.
+If not specified, default links are used. For details of all supported
+families see \code{\link[=brmsfamily]{brmsfamily()}}. Commonly used, such as \code{\link[=lognormal]{lognormal()}}, are also
+reexported as part of \code{epidist}.}
+}
+\value{
+A function that takes a prep argument from brms and returns a matrix
+of posterior expected values, with one row per posterior draw and one column
+per observation.
+}
+\description{
+This function creates a function that calculates the expected value of the
+posterior predictive distribution for a latent model. The returned function
+takes a \code{prep} argument (from brms) and returns posterior expected values.
+This is used internally by \code{\link[brms:posterior_epred.brmsfit]{brms::posterior_epred()}} to calculate expected
+values for latent models.
+}
+\seealso{
+\code{\link[brms:posterior_epred.brmsfit]{brms::posterior_epred()}} for details on how this is used within
+\code{brms}.
+
+Other gen: 
+\code{\link{epidist_gen_posterior_predict}()}
+}
+\concept{gen}
diff --git a/man/epidist_gen_posterior_predict.Rd b/man/epidist_gen_posterior_predict.Rd
new file mode 100644
index 000000000..b572d4c6e
--- /dev/null
+++ b/man/epidist_gen_posterior_predict.Rd
@@ -0,0 +1,42 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/gen.R
+\name{epidist_gen_posterior_predict}
+\alias{epidist_gen_posterior_predict}
+\title{Create a function to draw from the posterior predictive distribution for a
+latent model}
+\usage{
+epidist_gen_posterior_predict(family)
+}
+\arguments{
+\item{family}{A description of the response distribution and link function to
+be used in the model. Every family function has a link argument allowing
+users to specify the link function to be applied on the response variable.
+If not specified, default links are used. For details of all supported
+families see \code{\link[=brmsfamily]{brmsfamily()}}. Commonly used, such as \code{\link[=lognormal]{lognormal()}}, are also
+reexported as part of \code{epidist}.}
+}
+\value{
+A function that takes a \code{prep} argument from brms and returns a
+matrix of posterior predictions, with one row per posterior draw and one
+column per observation. The \code{prep} object must have the following variables:
+\itemize{
+\item \code{vreal1}: relative observation time
+\item \code{vreal2}: primary event window
+\item \code{vreal3}: secondary event window
+}
+}
+\description{
+This function creates a function that draws from the posterior predictive
+distribution for a latent model using \code{\link[primarycensored:rprimarycensored]{primarycensored::rpcens()}} to handle
+censoring and truncation. The returned function takes a \code{prep} argument from
+\code{brms} and returns posterior predictions. This is used internally by
+\code{\link[brms:posterior_predict.brmsfit]{brms::posterior_predict()}} to generate predictions for latent models.
+}
+\seealso{
+\code{\link[brms:posterior_predict.brmsfit]{brms::posterior_predict()}} for details on how this is used within
+\code{brms}, \code{\link[primarycensored:rprimarycensored]{primarycensored::rpcens()}} for details on the censoring approach
+
+Other gen: 
+\code{\link{epidist_gen_posterior_epred}()}
+}
+\concept{gen}
diff --git a/man/is_epidist_latent_model.Rd b/man/is_epidist_latent_model.Rd
index b9c4115ae..647725836 100644
--- a/man/is_epidist_latent_model.Rd
+++ b/man/is_epidist_latent_model.Rd
@@ -18,6 +18,7 @@ Other latent_model:
 \code{\link{as_epidist_latent_model.epidist_linelist_data}()},
 \code{\link{epidist_family_model.epidist_latent_model}()},
 \code{\link{epidist_formula_model.epidist_latent_model}()},
+\code{\link{epidist_gen_log_lik_latent}()},
 \code{\link{new_epidist_latent_model}()}
 }
 \concept{latent_model}
diff --git a/man/log_lik_latent_gamma.Rd b/man/log_lik_latent_gamma.Rd
deleted file mode 100644
index 5b460412b..000000000
--- a/man/log_lik_latent_gamma.Rd
+++ /dev/null
@@ -1,17 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/latent_gamma.R
-\name{log_lik_latent_gamma}
-\alias{log_lik_latent_gamma}
-\title{Calculate the pointwise log likelihood of the \code{latent_gamma} family}
-\usage{
-log_lik_latent_gamma(i, prep)
-}
-\arguments{
-\item{i}{The index of the observation to calculate the log likelihood of}
-
-\item{prep}{The result of a call to \code{\link[brms:prepare_predictions]{brms::prepare_predictions()}}}
-}
-\description{
-See \code{\link[brms:log_lik.brmsfit]{brms::log_lik()}}.
-}
-\keyword{internal}
diff --git a/man/log_lik_latent_lognormal.Rd b/man/log_lik_latent_lognormal.Rd
deleted file mode 100644
index 4cac6907c..000000000
--- a/man/log_lik_latent_lognormal.Rd
+++ /dev/null
@@ -1,17 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/latent_lognormal.R
-\name{log_lik_latent_lognormal}
-\alias{log_lik_latent_lognormal}
-\title{Calculate the pointwise log likelihood of the \code{latent_gamma} family}
-\usage{
-log_lik_latent_lognormal(i, prep)
-}
-\arguments{
-\item{i}{The index of the observation to calculate the log likelihood of}
-
-\item{prep}{The result of a call to \code{\link[brms:prepare_predictions]{brms::prepare_predictions()}}}
-}
-\description{
-See \code{\link[brms:log_lik.brmsfit]{brms::log_lik()}}.
-}
-\keyword{internal}
diff --git a/man/new_epidist_latent_model.Rd b/man/new_epidist_latent_model.Rd
index 02841df23..fe5e672a7 100644
--- a/man/new_epidist_latent_model.Rd
+++ b/man/new_epidist_latent_model.Rd
@@ -21,6 +21,7 @@ Other latent_model:
 \code{\link{as_epidist_latent_model.epidist_linelist_data}()},
 \code{\link{epidist_family_model.epidist_latent_model}()},
 \code{\link{epidist_formula_model.epidist_latent_model}()},
+\code{\link{epidist_gen_log_lik_latent}()},
 \code{\link{is_epidist_latent_model}()}
 }
 \concept{latent_model}
diff --git a/man/posterior_epred_latent_gamma.Rd b/man/posterior_epred_latent_gamma.Rd
deleted file mode 100644
index 1bc3a4e14..000000000
--- a/man/posterior_epred_latent_gamma.Rd
+++ /dev/null
@@ -1,16 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/latent_gamma.R
-\name{posterior_epred_latent_gamma}
-\alias{posterior_epred_latent_gamma}
-\title{Draws from the expected value of the posterior predictive distribution of the
-\code{latent_gamma} family}
-\usage{
-posterior_epred_latent_gamma(prep)
-}
-\arguments{
-\item{prep}{The result of a call to \code{\link[brms:prepare_predictions]{brms::prepare_predictions}}}
-}
-\description{
-See \code{\link[brms:posterior_epred.brmsfit]{brms::posterior_epred()}}.
-}
-\keyword{internal}
diff --git a/man/posterior_epred_latent_lognormal.Rd b/man/posterior_epred_latent_lognormal.Rd
deleted file mode 100644
index c8fb0d6ac..000000000
--- a/man/posterior_epred_latent_lognormal.Rd
+++ /dev/null
@@ -1,16 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/latent_lognormal.R
-\name{posterior_epred_latent_lognormal}
-\alias{posterior_epred_latent_lognormal}
-\title{Draws from the expected value of the posterior predictive distribution of the
-\code{latent_gamma} family}
-\usage{
-posterior_epred_latent_lognormal(prep)
-}
-\arguments{
-\item{prep}{The result of a call to \code{\link[brms:prepare_predictions]{brms::prepare_predictions}}}
-}
-\description{
-See \code{\link[brms:posterior_epred.brmsfit]{brms::posterior_epred()}}.
-}
-\keyword{internal}
diff --git a/man/posterior_predict_latent_gamma.Rd b/man/posterior_predict_latent_gamma.Rd
deleted file mode 100644
index cc0a1d3dd..000000000
--- a/man/posterior_predict_latent_gamma.Rd
+++ /dev/null
@@ -1,19 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/latent_gamma.R
-\name{posterior_predict_latent_gamma}
-\alias{posterior_predict_latent_gamma}
-\title{Draws from the posterior predictive distribution of the \code{latent_gamma} family}
-\usage{
-posterior_predict_latent_gamma(i, prep, ...)
-}
-\arguments{
-\item{i}{The index of the observation to predict}
-
-\item{prep}{The result of a call to \code{\link[brms:posterior_predict.brmsfit]{brms::posterior_predict()}}}
-
-\item{...}{Additional arguments}
-}
-\description{
-See \code{\link[brms:posterior_predict.brmsfit]{brms::posterior_predict()}}.
-}
-\keyword{internal}
diff --git a/man/posterior_predict_latent_lognormal.Rd b/man/posterior_predict_latent_lognormal.Rd
deleted file mode 100644
index bb5612305..000000000
--- a/man/posterior_predict_latent_lognormal.Rd
+++ /dev/null
@@ -1,20 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/latent_lognormal.R
-\name{posterior_predict_latent_lognormal}
-\alias{posterior_predict_latent_lognormal}
-\title{Draws from the posterior predictive distribution of the \code{latent_lognormal}
-family}
-\usage{
-posterior_predict_latent_lognormal(i, prep, ...)
-}
-\arguments{
-\item{i}{The index of the observation to predict}
-
-\item{prep}{The result of a call to \code{\link[brms:posterior_predict.brmsfit]{brms::posterior_predict()}}}
-
-\item{...}{Additional arguments}
-}
-\description{
-See \code{\link[brms:posterior_predict.brmsfit]{brms::posterior_predict()}}.
-}
-\keyword{internal}
diff --git a/tests/testthat/test-gen.R b/tests/testthat/test-gen.R
new file mode 100644
index 000000000..5dc4c5686
--- /dev/null
+++ b/tests/testthat/test-gen.R
@@ -0,0 +1,117 @@
+test_that("epidist_gen_posterior_predict returns a function that outputs positive integers with length equal to draws", { # nolint: line_length_linter.
+  skip_on_cran()
+  # Test lognormal
+  prep <- brms::prepare_predictions(fit)
+  i <- 1
+  predict_fn <- epidist_gen_posterior_predict(lognormal())
+  pred_i <- predict_fn(i = i, prep)
+  expect_identical(floor(pred_i), pred_i)
+  expect_length(pred_i, prep$ndraws)
+  expect_gte(min(pred_i), 0)
+
+  # Test gamma
+  prep_gamma <- brms::prepare_predictions(fit_gamma)
+  predict_fn_gamma <- epidist_gen_posterior_predict(Gamma())
+  pred_i_gamma <- predict_fn_gamma(i = i, prep_gamma)
+  expect_identical(floor(pred_i_gamma), pred_i_gamma)
+  expect_length(pred_i_gamma, prep_gamma$ndraws)
+  expect_gte(min(pred_i_gamma), 0)
+})
+
+test_that("epidist_gen_posterior_predict returns a function that errors for i out of bounds", { # nolint: line_length_linter.
+  skip_on_cran()
+  # Test lognormal
+  prep <- brms::prepare_predictions(fit)
+  i_out_of_bounds <- length(prep$data$Y) + 1
+  predict_fn <- epidist_gen_posterior_predict(lognormal())
+  expect_warning(
+    expect_error(
+      predict_fn(i = i_out_of_bounds, prep)
+    )
+  )
+
+  # Test gamma
+  prep_gamma <- brms::prepare_predictions(fit_gamma)
+  i_out_of_bounds_gamma <- length(prep_gamma$data$Y) + 1
+  predict_fn_gamma <- epidist_gen_posterior_predict(Gamma())
+  expect_warning(
+    expect_error(predict_fn_gamma(i = i_out_of_bounds_gamma, prep_gamma))
+  )
+})
+
+test_that("epidist_gen_posterior_predict returns a function that can generate predictions with no censoring", { # nolint: line_length_linter.
+  skip_on_cran()
+  # Test lognormal
+  predict_fn <- epidist_gen_posterior_predict(lognormal())
+  draws <- data.frame(relative_obs_time = 1000, pwindow = 0, swindow = 0) |>
+    tidybayes::add_predicted_draws(fit, ndraws = 100)
+  expect_identical(draws$.draw, 1:100)
+  pred <- draws$.prediction
+  expect_gte(min(pred), 0)
+  expect_true(all(abs(pred - round(pred)) > .Machine$double.eps^0.5))
+
+  # Test gamma
+  predict_fn_gamma <- epidist_gen_posterior_predict(Gamma())
+  draws_gamma <- data.frame(
+    relative_obs_time = 1000, pwindow = 0, swindow = 0
+  ) |>
+    tidybayes::add_predicted_draws(fit_gamma, ndraws = 100)
+  expect_identical(draws_gamma$.draw, 1:100)
+  pred_gamma <- draws_gamma$.prediction
+  expect_gte(min(pred_gamma), 0)
+  expect_true(
+    all(abs(pred_gamma - round(pred_gamma)) > .Machine$double.eps^0.5)
+  )
+})
+
+test_that("epidist_gen_posterior_predict returns a function that predicts delays in the 95% credible interval", { # nolint: line_length_linter.
+  skip_on_cran()
+  # Test lognormal
+  prep <- brms::prepare_predictions(fit)
+  prep$ndraws <- 1000 # Down from the 4000 for time saving
+  predict_fn <- epidist_gen_posterior_predict(lognormal())
+  q <- purrr::map_vec(seq_along(prep$data$Y), function(i) {
+    y <- predict_fn(i, prep)
+    ecdf <- ecdf(y)
+    q <- ecdf(prep$data$Y[i])
+    return(q)
+  })
+  expect_lt(quantile(q, 0.1), 0.3)
+  expect_gt(quantile(q, 0.9), 0.7)
+  expect_lt(min(q), 0.1)
+  expect_gt(max(q), 0.9)
+  expect_lt(mean(q), 0.65)
+  expect_gt(mean(q), 0.35)
+
+  # Test gamma
+  prep_gamma <- brms::prepare_predictions(fit_gamma)
+  prep_gamma$ndraws <- 1000
+  predict_fn_gamma <- epidist_gen_posterior_predict(Gamma())
+  q_gamma <- purrr::map_vec(seq_along(prep_gamma$data$Y), function(i) {
+    y <- predict_fn_gamma(i, prep_gamma)
+    ecdf <- ecdf(y)
+    q <- ecdf(prep_gamma$data$Y[i])
+    return(q)
+  })
+  expect_lt(quantile(q_gamma, 0.1), 0.3)
+  expect_gt(quantile(q_gamma, 0.9), 0.7)
+  expect_lt(min(q_gamma), 0.1)
+  expect_gt(max(q_gamma), 0.9)
+  expect_lt(mean(q_gamma), 0.65)
+  expect_gt(mean(q_gamma), 0.35)
+})
+
+test_that("epidist_gen_posterior_epred returns a function that creates arrays with correct dimensions", { # nolint: line_length_linter.
+  skip_on_cran()
+  # Test lognormal
+  epred <- prep_obs |>
+    tidybayes::add_epred_draws(fit)
+  expect_equal(mean(epred$.epred), 5.97, tolerance = 0.1)
+  expect_gte(min(epred$.epred), 0)
+
+  # Test gamma
+  epred_gamma <- prep_obs |>
+    tidybayes::add_epred_draws(fit_gamma)
+  expect_equal(mean(epred_gamma$.epred), 6.56, tolerance = 0.1)
+  expect_gte(min(epred_gamma$.epred), 0)
+})
diff --git a/tests/testthat/test-latent_gamma.R b/tests/testthat/test-latent_gamma.R
deleted file mode 100644
index e91bbf4a5..000000000
--- a/tests/testthat/test-latent_gamma.R
+++ /dev/null
@@ -1,64 +0,0 @@
-test_that("posterior_predict_latent_gamma outputs positive integers with length equal to draws", { # nolint: line_length_linter.
-  skip_on_cran()
-  prep <- brms::prepare_predictions(fit_gamma)
-  i <- 1
-  pred_i <- posterior_predict_latent_gamma(i = i, prep)
-  expect_identical(floor(pred_i), pred_i)
-  expect_length(pred_i, prep$ndraws)
-  expect_gte(min(pred_i), 0)
-})
-
-test_that("posterior_predict_latent_gamma errors for i out of bounds", { # nolint: line_length_linter.
-  skip_on_cran()
-  prep <- brms::prepare_predictions(fit_gamma)
-  i_out_of_bounds <- length(prep$data$Y) + 1
-  expect_error(posterior_predict_latent_gamma(i = i_out_of_bounds, prep))
-})
-
-test_that("posterior_predict_latent_gamma can generate predictions with no censoring", { # nolint: line_length_linter.
-  skip_on_cran()
-  draws <- data.frame(relative_obs_time = 1000, pwindow = 0, swindow = 0) |>
-    tidybayes::add_predicted_draws(fit_gamma, ndraws = 100)
-  expect_identical(draws$.draw, 1:100)
-  pred <- draws$.prediction
-  expect_gte(min(pred), 0)
-  expect_true(all(abs(pred - round(pred)) > .Machine$double.eps^0.5))
-})
-
-test_that("posterior_predict_latent_gamma predicts delays for which the data is in the 95% credible interval", { # nolint: line_length_linter.
-  skip_on_cran()
-  prep <- brms::prepare_predictions(fit_gamma)
-  prep$ndraws <- 1000 # Down from the 4000 for time saving
-  q <- purrr::map_vec(seq_along(prep$data$Y), function(i) {
-    y <- posterior_predict_latent_gamma(i, prep)
-    ecdf <- ecdf(y)
-    q <- ecdf(prep$data$Y[i])
-    return(q)
-  })
-  expect_lt(quantile(q, 0.1), 0.3)
-  expect_gt(quantile(q, 0.9), 0.7)
-  expect_lt(min(q), 0.1)
-  expect_gt(max(q), 0.9)
-  expect_lt(mean(q), 0.65)
-  expect_gt(mean(q), 0.35)
-})
-
-test_that("posterior_epred_latent_gamma creates a array of non-negative numbers with the correct dimensions", { # nolint: line_length_linter.
-  skip_on_cran()
-  prep <- brms::prepare_predictions(fit_gamma)
-  epred <- posterior_epred_latent_gamma(prep)
-  expect_setequal(class(epred), c("matrix", "array"))
-  expect_identical(nrow(epred), prep$ndraws)
-  expect_identical(ncol(epred), length(prep$data$Y))
-  expect_gte(min(epred), 0)
-})
-
-test_that("log_lik_latent_gamma produces a vector with length ndraws of finite non-NA numbers", { # nolint: line_length_linter.
-  skip_on_cran()
-  prep <- brms::prepare_predictions(fit_gamma)
-  i <- 1
-  log_lik <- log_lik_latent_gamma(i, prep)
-  expect_length(log_lik, prep$ndraws)
-  expect_false(anyNA(log_lik))
-  expect_true(all(is.finite(log_lik)))
-})
diff --git a/tests/testthat/test-latent_lognormal.R b/tests/testthat/test-latent_lognormal.R
deleted file mode 100644
index 5ec8d7489..000000000
--- a/tests/testthat/test-latent_lognormal.R
+++ /dev/null
@@ -1,64 +0,0 @@
-test_that("posterior_predict_latent_lognormal outputs positive integers with length equal to draws", { # nolint: line_length_linter.
-  skip_on_cran()
-  prep <- brms::prepare_predictions(fit)
-  i <- 1
-  pred_i <- posterior_predict_latent_lognormal(i = i, prep)
-  expect_identical(floor(pred_i), pred_i)
-  expect_length(pred_i, prep$ndraws)
-  expect_gte(min(pred_i), 0)
-})
-
-test_that("posterior_predict_latent_lognormal errors for i out of bounds", { # nolint: line_length_linter.
-  skip_on_cran()
-  prep <- brms::prepare_predictions(fit)
-  i_out_of_bounds <- length(prep$data$Y) + 1
-  expect_error(posterior_predict_latent_lognormal(i = i_out_of_bounds, prep))
-})
-
-test_that("posterior_predict_latent_lognormal can generate predictions with no censoring", { # nolint: line_length_linter.
-  skip_on_cran()
-  draws <- data.frame(relative_obs_time = 1000, pwindow = 0, swindow = 0) |>
-    tidybayes::add_predicted_draws(fit, ndraws = 100)
-  expect_identical(draws$.draw, 1:100)
-  pred <- draws$.prediction
-  expect_gte(min(pred), 0)
-  expect_true(all(abs(pred - round(pred)) > .Machine$double.eps^0.5))
-})
-
-test_that("posterior_predict_latent_lognormal predicts delays for which the data is in the 95% credible interval", { # nolint: line_length_linter.
-  skip_on_cran()
-  prep <- brms::prepare_predictions(fit)
-  prep$ndraws <- 1000 # Down from the 4000 for time saving
-  q <- purrr::map_vec(seq_along(prep$data$Y), function(i) {
-    y <- posterior_predict_latent_lognormal(i, prep)
-    ecdf <- ecdf(y)
-    q <- ecdf(prep$data$Y[i])
-    return(q)
-  })
-  expect_lt(quantile(q, 0.1), 0.3)
-  expect_gt(quantile(q, 0.9), 0.7)
-  expect_lt(min(q), 0.1)
-  expect_gt(max(q), 0.9)
-  expect_lt(mean(q), 0.65)
-  expect_gt(mean(q), 0.35)
-})
-
-test_that("posterior_epred_latent_lognormal creates a array of non-negative numbers with the correct dimensions", { # nolint: line_length_linter.
-  skip_on_cran()
-  prep <- brms::prepare_predictions(fit)
-  epred <- posterior_epred_latent_lognormal(prep)
-  expect_setequal(class(epred), c("matrix", "array"))
-  expect_identical(nrow(epred), prep$ndraws)
-  expect_identical(ncol(epred), length(prep$data$Y))
-  expect_gte(min(epred), 0)
-})
-
-test_that("log_lik_latent_lognormal produces a vector with length ndraws of finite non-NA numbers", { # nolint: line_length_linter.
-  skip_on_cran()
-  prep <- brms::prepare_predictions(fit)
-  i <- 1
-  log_lik <- log_lik_latent_lognormal(i, prep)
-  expect_length(log_lik, prep$ndraws)
-  expect_false(anyNA(log_lik))
-  expect_true(all(is.finite(log_lik)))
-})
diff --git a/tests/testthat/test-latent_model.R b/tests/testthat/test-latent_model.R
index 579af69fc..c1a64560f 100644
--- a/tests/testthat/test-latent_model.R
+++ b/tests/testthat/test-latent_model.R
@@ -56,3 +56,23 @@ test_that("epidist_stancode.epidist_latent_model produces valid stanvars", { # n
   )
   expect_s3_class(stancode, "stanvars")
 })
+
+test_that("epidist_gen_log_lik_latent returns a function that produces valid log likelihoods", { # nolint: line_length_linter.
+  skip_on_cran()
+  # Test lognormal
+  prep <- brms::prepare_predictions(fit)
+  i <- 1
+  log_lik_fn <- epidist_gen_log_lik_latent(lognormal())
+  log_lik <- log_lik_fn(i = i, prep)
+  expect_length(log_lik, prep$ndraws)
+  expect_false(anyNA(log_lik))
+  expect_true(all(is.finite(log_lik)))
+
+  # Test gamma
+  prep_gamma <- brms::prepare_predictions(fit_gamma)
+  log_lik_fn_gamma <- epidist_gen_log_lik_latent(Gamma())
+  log_lik_gamma <- log_lik_fn_gamma(i = i, prep_gamma)
+  expect_length(log_lik_gamma, prep_gamma$ndraws)
+  expect_false(anyNA(log_lik_gamma))
+  expect_true(all(is.finite(log_lik_gamma)))
+})
diff --git a/vignettes/ebola.Rmd b/vignettes/ebola.Rmd
index 0e3592597..bd1c199fe 100644
--- a/vignettes/ebola.Rmd
+++ b/vignettes/ebola.Rmd
@@ -239,8 +239,9 @@ fit <- epidist(
   formula = mu ~ 1,
   family = lognormal(),
   algorithm = "sampling",
-  refresh = 0,
-  silent = 2,
+  chains = 2,
+  cores = 2,
+  refresh = as.integer(interactive()),
   seed = 1,
   backend = "cmdstanr"
 )
@@ -264,8 +265,9 @@ fit_sex <- epidist(
   formula = bf(mu ~ 1 + sex, sigma ~ 1 + sex),
   family = lognormal(),
   algorithm = "sampling",
-  refresh = 0,
-  silent = 2,
+  chains = 2,
+  cores = 2,
+  refresh = as.integer(interactive()),
   seed = 1,
   backend = "cmdstanr"
 )
@@ -294,8 +296,9 @@ fit_sex_district <- epidist(
   ),
   family = lognormal(),
   algorithm = "sampling",
-  refresh = 0,
-  silent = 2,
+  chains = 2,
+  cores = 2,
+  refresh = as.integer(interactive()),
   seed = 1,
   backend = "cmdstanr"
 )
@@ -312,14 +315,13 @@ ranef(fit_sex_district)
 
 To go further than summaries of the fitted model, we recommend using the `tidybayes` package.
 For example, to obtain the posterior expectation of the delay distribution, under no censoring or truncation, we may use the `modelr::data_grid()` function in combination with the `tidybayes::add_epred_draws()` function.
-The `tidybayes::add_epred_draws()` function uses the (internal) `posterior_epred_latent_lognormal()` function implemented in `epidist` for the `latent_lognormal` custom `brms` family.
+The `tidybayes::add_epred_draws()` function uses the `epidist_gen_posterior_predict()` function to generate a posterior prediction function for the `lognormal()` distribution.
 
 In Figure \@ref(fig:epred) we show the posterior expectation of the delay distribution for each of the three fitted models.
 Figure \@ref(fig:epred)B illustrates the higher mean of men as compared with women.
 
 ```{r}
 epred_draws <- obs_prep |>
-  as.data.frame() |>
   data_grid(NA) |>
   mutate(relative_obs_time = NA, pwindow = NA, swindow = NA) |>
   add_epred_draws(fit, dpar = TRUE)
@@ -331,7 +333,6 @@ epred_base_figure <- epred_draws |>
   theme_minimal()
 
 epred_draws_sex <- obs_prep |>
-  as.data.frame() |>
   data_grid(sex) |>
   mutate(relative_obs_time = NA, pwindow = NA, swindow = NA) |>
   add_epred_draws(fit_sex, dpar = TRUE)
@@ -343,7 +344,6 @@ epred_sex_figure <- epred_draws_sex |>
   theme_minimal()
 
 epred_draws_sex_district <- obs_prep |>
-  as.data.frame() |>
   data_grid(sex, district) |>
   mutate(relative_obs_time = NA, pwindow = NA, swindow = NA) |>
   add_epred_draws(fit_sex_district, dpar = TRUE)
diff --git a/vignettes/epidist.Rmd b/vignettes/epidist.Rmd
index 86fa9b254..399652a84 100644
--- a/vignettes/epidist.Rmd
+++ b/vignettes/epidist.Rmd
@@ -285,9 +285,14 @@ The parameters of the model are inferred using Bayesian inference.
 In particular, we use the the No-U-Turn Sampler (NUTS) Markov chain Monte Carlo (MCMC) algorithm via the [`brms`](https://paul-buerkner.github.io/brms/) R package [@brms].
 
 ```{r}
-fit <- epidist(data = data, chains = 2, cores = 2, refresh = 0)
+fit <- epidist(
+  data = data, chains = 2, cores = 2, refresh = as.integer(interactive())
+)
 ```
 
+**Note that here we use the default `rstan` backend but we generally recommend using the `cmdstanr` backend for faster sampling and additional features.**
+**This can be set using `backend = "cmdstanr"` after following the installing CmdStan instructions in the README.**
+
 The `fit` object is a `brmsfit` object containing MCMC samples from each of the parameters in the model, shown in the table below.
 Many of these parameters (e.g. `swindow` and `pwindow`) are the so called latent variables, and have lengths corresponding to the `sample_size`.
 
@@ -305,8 +310,13 @@ data.frame(
 Users familiar with Stan and `brms`, can work with `fit` directly.
 Any tool that supports `brms` fitted model objects will be compatible with `fit`.
 
-**Note that here we use the default `rstan` backend but we generally recommend using the `cmdstanr` backend for faster sampling and additional features.**
-**This can be set using `backend = "cmdstanr"` after following the installing CmdStan instructions in the README.**
+For example, we can use the built in `summary()` function
+
+```{r}
+summary(fit)
+```
+
+to see that the model has converged and recovered the delay distribution distribution paramters reasonably well.
 
 The `epidist` package also provides functions to make common post-processing tasks easy.
 For example, individual predictions of the lognormal delay parameters can be extracted using:
@@ -315,31 +325,42 @@ For example, individual predictions of the lognormal delay parameters can be ext
 pred <- predict_delay_parameters(fit)
 ```
 
-Figure \@ref(fig:fitted-lognormal) shows the lognormal delay distribution obtained using the average of the `mu` and `sigma` draws.
+Figure \@ref(fig:fitted-lognormal) shows the lognormal delay distribution obtained for 100 draws from the posterior distribution of the `mu` and `sigma` parameters.
 Whereas in Figure \@ref(fig:obs-est) the histogram of censored, truncated, sampled data was substantially different to the underlying delay distribution, using `epidist()` we have obtained a much closer match to the truth.
 
-(ref:fitted-lognormal) A fitted delay distribution (in pink) as compared with the true delay distribution (in black).
+(ref:fitted-lognormal) The fitted delay distribution (pink lines, 100 draws from the posterior) compared to the true underlying delay distribution (black line). Each pink line represents one possible delay distribution based on the fitted model parameters.
 
 ```{r fitted-lognormal, fig.cap="(ref:fitted-lognormal)"}
+# Sample 100 random draws from the posterior
+set.seed(123)
+sample_draws <- sample(nrow(pred), 100)
+
 ggplot() +
+  # Plot the true distribution
   geom_function(
     data = data.frame(x = c(0, 30)),
     aes(x = x),
     fun = dlnorm,
+    linewidth = 1.5,
     args = list(
       meanlog = secondary_dist[["mu"]],
       sdlog = secondary_dist[["sigma"]]
     )
   ) +
-  geom_function(
-    data = data.frame(x = c(0, 30)),
-    aes(x = x), fun = dlnorm,
-    args = list(
-      meanlog = mean(pred$mu),
-      sdlog = mean(pred$sigma)
-    ),
-    col = "#CC79A7"
-  ) +
+  # Plot 100 sampled posterior distributions
+  lapply(sample_draws, \(i) {
+    geom_function(
+      data = data.frame(x = c(0, 30)),
+      aes(x = x),
+      fun = dlnorm,
+      args = list(
+        meanlog = pred$mu[i],
+        sdlog = pred$sigma[i]
+      ),
+      alpha = 0.1,
+      color = "#CC79A7"
+    )
+  }) +
   labs(
     x = "Delay between primary and secondary event (days)",
     y = "Probability density"