diff --git a/vignettes/approx-inference.Rmd b/vignettes/approx-inference.Rmd
new file mode 100644
index 000000000..6793965af
--- /dev/null
+++ b/vignettes/approx-inference.Rmd
@@ -0,0 +1,179 @@
+---
+title: "Approximate Bayesian inference in `epidist`"
+description: "..."
+output: 
+  bookdown::html_document2:
+    fig_caption: yes
+    code_folding: show
+    number_sections: true
+pkgdown:
+  as_is: true
+# csl: https://raw.githubusercontent.com/citation-style-language/styles/master/apa-numeric-superscript-brackets.csl
+link-citations: true
+vignette: >
+  %\VignetteIndexEntry{Getting started with epidist}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+bibliography: references.bib
+---
+
+```{r setup, include=FALSE}
+# exclude compile warnings from cmdstanr
+knitr::opts_chunk$set(
+  fig.path = "figures/epidist-",
+  cache = TRUE,
+  collapse = TRUE,
+  comment = "#>",
+  message = FALSE,
+  warning = FALSE,
+  error = FALSE
+)
+```
+
+```{r load-requirements}
+library(epidist)
+library(data.table)
+library(purrr)
+library(ggplot2)
+library(gt)
+library(dplyr)
+```
+
+```{r}
+meanlog <- 1.8
+sdlog <- 0.5
+obs_time <- 25
+sample_size <- 200
+
+obs_cens_trunc <- simulate_gillespie() |>
+  simulate_secondary(
+    meanlog = meanlog,
+    sdlog = sdlog
+  ) |>
+  observe_process() |>
+  filter_obs_by_obs_time(obs_time = obs_time)
+
+obs_cens_trunc_samp <- obs_cens_trunc[sample(seq_len(.N), sample_size, replace = FALSE)]
+
+data <- obs_cens_trunc_samp
+
+formula <- brms::bf(delay_central | vreal(obs_t, pwindow_upr, swindow_upr) ~ 1, sigma ~ 1)
+
+fn <- brms::brm
+
+family <- brms::custom_family(
+  "latent_lognormal",
+  dpars = c("mu", "sigma"),
+  links = c("identity", "log"),
+  lb = c(NA, 0),
+  ub = c(NA, NA),
+  type = "real",
+  vars = c("pwindow", "swindow", "vreal1"),
+  loop = FALSE
+)
+
+scode_functions <- "
+   \n real latent_lognormal_lpdf(vector y, vector mu, vector sigma,
+   \n vector pwindow, vector swindow,
+   \n array[] real obs_t) {
+   \n      int n = num_elements(y);
+   \n      vector[n] d = y - pwindow + swindow;
+   \n      vector[n] obs_time = to_vector(obs_t) - pwindow;
+   \n      return lognormal_lpdf(d | mu, sigma) -
+   \n        lognormal_lcdf(obs_time | mu, sigma);
+   \n      }
+   \n
+"
+
+scode_parameters <- "
+   \n    vector<lower = 0, upper = 1>[N] swindow_raw;
+   \n    vector<lower = 0, upper = 1>[N] pwindow_raw;
+   \n
+"
+
+scode_tparameters <- "
+   \n    vector<lower = 0>[N] pwindow;\n    vector<lower = 0>[N] swindow;
+   \n    swindow = to_vector(vreal3) .* swindow_raw;
+   \n    pwindow[noverlap] = to_vector(vreal2[noverlap]) .* pwindow_raw[noverlap];
+   \n    if (wN) {
+   \n      pwindow[woverlap] = swindow[woverlap] .* pwindow_raw[woverlap];
+   \n    }
+   \n
+"
+
+scode_priors <- "
+   \n    swindow_raw ~ uniform(0, 1);
+   \n    pwindow_raw ~ uniform(0, 1);
+   \n
+"
+
+data <- data.table::as.data.table(data)
+data[, id := 1:.N]
+data[, obs_t := obs_at - ptime_lwr]
+
+data[, pwindow_upr := ifelse(
+  stime_lwr < ptime_upr,
+  stime_upr - ptime_lwr,
+  ptime_upr - ptime_lwr
+)]
+
+data[, woverlap := as.numeric(stime_lwr < ptime_upr)]
+data[, swindow_upr := stime_upr - stime_lwr]
+data[, delay_central := stime_lwr - ptime_lwr]
+data[, row_id := 1:.N]
+
+if (nrow(data) > 1) {
+  data <- data[, id := as.factor(id)]
+}
+
+stanvars_functions <- brms::stanvar(
+  block = "functions", scode = scode_functions
+)
+
+stanvars_data <- brms::stanvar(
+  block = "data", scode = "int wN;",
+  x = nrow(data[woverlap > 0]),
+  name = "wN"
+) +
+
+brms::stanvar(
+  block = "data", scode = "array[N - wN] int noverlap;",
+  x = data[woverlap == 0][, row_id],
+  name = "noverlap"
+) +
+brms::stanvar(
+  block = "data", scode = "array[wN] int woverlap;",
+  x = data[woverlap > 0][, row_id],
+  name = "woverlap"
+)
+
+stanvars_parameters <- brms::stanvar(
+  block = "parameters", scode = scode_parameters
+)
+
+stanvars_tparameters <- brms::stanvar(
+  block = "tparameters", scode = scode_tparameters
+)
+
+stanvars_priors <- brms::stanvar(block = "model", scode = scode_priors)
+
+stanvars_all <- stanvars_functions + stanvars_data + stanvars_parameters + 
+  stanvars_tparameters + stanvars_priors
+
+fit_laplace <- fn(
+  formula = formula, family = family, stanvars = stanvars_all,
+  backend = "cmdstanr", data = data, algorithm = "laplace"
+)
+
+fit_pathfinder <- fn(
+  formula = formula, family = family, stanvars = stanvars_all,
+  backend = "cmdstanr", data = data, algorithm = "pathfinder"
+)
+
+fit_hmc <- fn(
+  formula = formula, family = family, stanvars = stanvars_all,
+  backend = "cmdstanr", data = data, algorithm = "sampling"
+)
+```
+
+## Bibliography {-}