diff --git a/DESCRIPTION b/DESCRIPTION
index 08c1f7dd3..41fd5e0d8 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -148,7 +148,7 @@ Encoding: UTF-8
 Language: en-GB
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 NeedsCompilation: yes
 SystemRequirements: GNU make
        C++17
diff --git a/NAMESPACE b/NAMESPACE
index 3d7391fe8..35d3d7585 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -227,6 +227,10 @@ importFrom(stats,pexp)
 importFrom(stats,pgamma)
 importFrom(stats,plnorm)
 importFrom(stats,pnorm)
+importFrom(stats,qexp)
+importFrom(stats,qgamma)
+importFrom(stats,qlnorm)
+importFrom(stats,qnorm)
 importFrom(stats,quantile)
 importFrom(stats,quasipoisson)
 importFrom(stats,rexp)
diff --git a/R/deprecated.R b/R/deprecated.R
index c7928e360..08abbdd70 100644
--- a/R/deprecated.R
+++ b/R/deprecated.R
@@ -447,3 +447,189 @@ rstan_opts <- function(object = NULL,
     "stan_opts()"
   )
 }
+
+#' Distribution Skeleton
+#'
+#' @description `r lifecycle::badge("deprecated")`
+#' This function acts as a skeleton for a truncated distribution defined by
+#' model type, maximum value and model parameters. It is designed to be used
+#' with the output from [get_dist()].
+#'
+#' @param n Numeric vector, number of samples to take (or days for the
+#' probability density).
+#'
+#' @param dist Logical, defaults to `FALSE`. Should the probability density be
+#' returned rather than a number of samples.
+#'
+#' @param cum Logical, defaults to `TRUE`. If `dist = TRUE` should the returned
+#' distribution be cumulative.
+#'
+#' @param model Character string, defining the model to be used. Supported
+#'   options are exponential ("exp"), gamma ("gamma"), and log normal
+#'   ("lognormal")
+#'
+#' @param discrete Logical,  defaults to `FALSE`. Should the probability
+#'   distribution be discretised. In this case each entry of the probability
+#'   mass function corresponds to the 2-length interval ending at the entry
+#'   except for the first interval that covers (0, 1).  That is, the probability
+#'   mass function is a vector where the first entry corresponds to the integral
+#'   over the (0,1] interval of the continuous distribution, the second entry
+#'   corresponds to the (0,2] interval, the third entry corresponds to the (1,
+#'   3] interval etc.
+#'
+#' @param params A list of parameters values (by name) required for each model.
+#' For the exponential model this is a rate parameter and for the gamma model
+#' this is alpha and beta.
+#'
+#' @param max_value Numeric, the maximum value to allow. Defaults to 120.
+#' Samples outside of this range are resampled.
+#'
+#' @return A vector of samples or a probability distribution.
+#' @export
+#' @examples
+#'
+#' ## Exponential model
+#' # sample
+#' dist_skel(10, model = "exp", params = list(rate = 1))
+#'
+#' # cumulative prob density
+#' dist_skel(1:10, model = "exp", dist = TRUE, params = list(rate = 1))
+#'
+#' # probability density
+#' dist_skel(1:10,
+#'   model = "exp", dist = TRUE,
+#'   cum = FALSE, params = list(rate = 1)
+#' )
+#'
+#' ## Gamma model
+#' # sample
+#' dist_skel(10, model = "gamma", params = list(shape = 1, rate = 0.5))
+#'
+#' # cumulative prob density
+#' dist_skel(0:10,
+#'   model = "gamma", dist = TRUE,
+#'   params = list(shape = 1, rate = 0.5)
+#' )
+#'
+#' # probability density
+#' dist_skel(0:10,
+#'   model = "gamma", dist = TRUE,
+#'   cum = FALSE, params = list(shape = 2, rate = 0.5)
+#' )
+#'
+#' ## Log normal model
+#' # sample
+#' dist_skel(10,
+#'   model = "lognormal", params = list(meanlog = log(5), sdlog = log(2))
+#' )
+#'
+#' # cumulative prob density
+#' dist_skel(0:10,
+#'   model = "lognormal", dist = TRUE,
+#'   params = list(meanlog = log(5), sdlog = log(2))
+#' )
+#'
+#' # probability density
+#' dist_skel(0:10,
+#'   model = "lognormal", dist = TRUE, cum = FALSE,
+#'   params = list(meanlog = log(5), sdlog = log(2))
+#' )
+dist_skel <- function(n, dist = FALSE, cum = TRUE, model,
+                      discrete = FALSE, params, max_value = 120) {
+  lifecycle::deprecate_warn(
+    "1.6.0", "dist_skel()"
+  )
+  ## define unnormalised support function
+  if (model == "exp") {
+    updist <- function(n) {
+      pexp(n, params[["rate"]])
+    }
+  } else if (model == "gamma") {
+    updist <- function(n) {
+      pgamma(n, params[["shape"]], params[["rate"]])
+    }
+  } else if (model == "lognormal") {
+    updist <- function(n) {
+      plnorm(n, params[["meanlog"]], params[["sdlog"]])
+    }
+  } else if (model == "normal") {
+    updist <- function(n) {
+      pnorm(n, params[["mean"]], params[["sd"]])
+    }
+  } else if (model == "fixed") {
+    updist <- function(n) {
+      as.integer(n > params[["value"]])
+    }
+  }
+
+  if (discrete) {
+    cmf <- c(0, updist(1),
+      updist(seq_len(max_value)) + updist(seq_len(max_value) + 1)
+    ) /
+      (updist(max_value) + updist(max_value + 1))
+    pmf <- diff(cmf)
+    rdist <- function(n) {
+      sample(
+        x = seq_len(max_value + 1) - 1, size = n, prob = pmf, replace = TRUE
+      )
+    }
+    pdist <- function(n) {
+      cmf[n + 1]
+    }
+    ddist <- function(n) {
+      pmf[n + 1]
+    }
+  } else {
+    pdist <- function(n) {
+      updist(n) / updist(max_value + 1)
+    }
+    ddist <- function(n) {
+      pdist(n + 1) - pdist(n)
+    }
+    if (model == "exp") {
+      rdist <- function(n) {
+        rexp(n, params[["rate"]])
+      }
+    } else if (model == "gamma") {
+      rdist <- function(n) {
+        rgamma(n, params[["shape"]], params[["rate"]])
+      }
+    } else if (model == "lognormal") {
+      rdist <- function(n) {
+        rlnorm(n, params[["meanlog"]], params[["sdlog"]])
+      }
+    }
+  }
+
+  # define internal sampling function
+  inner_skel <- function(n, dist = FALSE, cum = TRUE, max_value = NULL) {
+    if (dist) {
+      if (cum) {
+        ret <- pdist(n)
+      } else {
+        ret <- ddist(n)
+      }
+      ret[ret > 1] <- NA_real_
+      return(ret)
+    } else {
+      rdist(n)
+    }
+  }
+
+  # define truncation wrapper
+  truncated_skel <- function(n, dist, cum, max_value) {
+    n <- inner_skel(n, dist, cum, max_value)
+    if (!dist) {
+      while (any(!is.na(n) & n >= max_value)) {
+        n <- ifelse(n >= max_value, inner_skel(n), n)
+      }
+
+      n <- as.integer(n)
+    }
+    return(n)
+  }
+
+  # call function
+  sample <- truncated_skel(n, dist = dist, cum = cum, max_value = max_value)
+  return(sample)
+}
diff --git a/R/dist_spec.R b/R/dist_spec.R
index 48f717257..359214e5e 100644
--- a/R/dist_spec.R
+++ b/R/dist_spec.R
@@ -1,31 +1,32 @@
-#' Distribution Skeleton
+#' Discretised probability mass function
 #'
 #' @description `r lifecycle::badge("questioning")`
-#' This function acts as a skeleton for a truncated distribution defined by
-#' model type, maximum value and model parameters. It is designed to be used
-#' with the output from [get_dist()].
-#'
-#' @param n Numeric vector, number of samples to take (or days for the
-#' probability density).
-#'
-#' @param dist Logical, defaults to `FALSE`. Should the probability density be
-#' returned rather than a number of samples.
-#'
-#' @param cum Logical, defaults to `TRUE`. If `dist = TRUE` should the returned
-#' distribution be cumulative.
-#'
-#' @param model Character string, defining the model to be used. Supported
-#'   options are exponential ("exp"), gamma ("gamma"), and log normal
-#'   ("lognormal")
-#'
-#' @param discrete Logical,  defaults to `FALSE`. Should the probability
-#'   distribution be discretised. In this case each entry of the probability
-#'   mass function corresponds to the 2-length interval ending at the entry
-#'   except for the first interval that covers (0, 1).  That is, the probability
-#'   mass function is a vector where the first entry corresponds to the integral
-#'   over the (0,1] interval of the continuous distribution, the second entry
-#'   corresponds to the (0,2] interval, the third entry corresponds to the (1,
-#'   3] interval etc.
+#' This function returns the probability mass function of a discretised and
+#' truncated distribution defined by distribution type, maximum value and model
+#' parameters.
+#'
+#' # Methodological details
+#'
+#' The probability mass function of the discretised probability distribution is
+#'   a vector where the first entry corresponds to the integral over the (0,1]
+#'   interval of the corresponding continuous distribution (probability of
+#'   integer 0), the second entry corresponds to the (0,2] interval (probability
+#'   mass of integer 1), the third entry corresponds to the (1, 3] interval
+#'   (probability mass of integer 2), etc. This approximates the true
+#'   probability mass function of a double censored distribution which arises
+#'   from the difference of two censored events.
+#'
+#' @references
+#' Charniga, K., et al. “Best practices for estimating and reporting
+#'   epidemiological delay distributions of infectious diseases using public
+#'   health surveillance and healthcare data”, *arXiv e-prints*, 2024.
+#'   <https://doi.org/10.48550/arXiv.2405.08841>
+#' Park,  S. W.,  et al.,  "Estimating epidemiological delay distributions for
+#'   infectious diseases", *medRxiv*, 2024.
+#'   <https://doi.org/10.1101/2024.01.12.24301247>
+#'
+#' @param distribution A character string representing the distribution to be
+#'   used (one of "exp", "gamma", "lognormal", "normal" or "fixed")
 #'
 #' @param params A list of parameters values (by name) required for each model.
 #' For the exponential model this is a rate parameter and for the gamma model
@@ -34,151 +35,67 @@
 #' @param max_value Numeric, the maximum value to allow. Defaults to 120.
 #' Samples outside of this range are resampled.
 #'
-#' @return A vector of samples or a probability distribution.
-#' @export
-#' @examples
-#'
-#' ## Exponential model
-#' # sample
-#' dist_skel(10, model = "exp", params = list(rate = 1))
-#'
-#' # cumulative prob density
-#' dist_skel(1:10, model = "exp", dist = TRUE, params = list(rate = 1))
-#'
-#' # probability density
-#' dist_skel(1:10,
-#'   model = "exp", dist = TRUE,
-#'   cum = FALSE, params = list(rate = 1)
-#' )
-#'
-#' ## Gamma model
-#' # sample
-#' dist_skel(10, model = "gamma", params = list(shape = 1, rate = 0.5))
-#'
-#' # cumulative prob density
-#' dist_skel(0:10,
-#'   model = "gamma", dist = TRUE,
-#'   params = list(shape = 1, rate = 0.5)
-#' )
-#'
-#' # probability density
-#' dist_skel(0:10,
-#'   model = "gamma", dist = TRUE,
-#'   cum = FALSE, params = list(shape = 2, rate = 0.5)
-#' )
-#'
-#' ## Log normal model
-#' # sample
-#' dist_skel(10,
-#'   model = "lognormal", params = list(meanlog = log(5), sdlog = log(2))
-#' )
-#'
-#' # cumulative prob density
-#' dist_skel(0:10,
-#'   model = "lognormal", dist = TRUE,
-#'   params = list(meanlog = log(5), sdlog = log(2))
-#' )
-#'
-#' # probability density
-#' dist_skel(0:10,
-#'   model = "lognormal", dist = TRUE, cum = FALSE,
-#'   params = list(meanlog = log(5), sdlog = log(2))
-#' )
-dist_skel <- function(n, dist = FALSE, cum = TRUE, model,
-                      discrete = FALSE, params, max_value = 120) {
-  ## define unnormalised support function
-  if (model == "exp") {
+#' @param width Numeric, the width of each discrete bin.
+#
+#' @return A vector representing a probability distribution.
+#' @keywords internal
+#' @inheritParams apply_tolerance
+#' @importFrom stats pexp pgamma plnorm pnorm qexp qgamma qlnorm qnorm
+#' @importFrom rlang arg_match
+discrete_pmf <- function(distribution =
+                           c("exp", "gamma", "lognormal", "normal", "fixed"),
+                         params, max_value, tolerance, width) {
+  distribution <- arg_match(distribution)
+  ## define unnormalised support function and cumulative density function
+  if (distribution == "exp") {
     updist <- function(n) {
       pexp(n, params[["rate"]])
     }
-  } else if (model == "gamma") {
+    qdist <- qexp
+  } else if (distribution == "gamma") {
     updist <- function(n) {
       pgamma(n, params[["shape"]], params[["rate"]])
     }
-  } else if (model == "lognormal") {
+    qdist <- qgamma
+  } else if (distribution == "lognormal") {
     updist <- function(n) {
       plnorm(n, params[["meanlog"]], params[["sdlog"]])
     }
-  } else if (model == "normal") {
+    qdist <- qlnorm
+  } else if (distribution == "normal") {
     updist <- function(n) {
       pnorm(n, params[["mean"]], params[["sd"]])
     }
-  } else if (model == "fixed") {
+    qdist <- qnorm
+  } else if (distribution == "fixed") {
     updist <- function(n) {
       as.integer(n > params[["value"]])
     }
+    qdist <- function(p, value) return(value)
   }
 
-  if (discrete) {
-    cmf <- c(0, updist(1),
-      updist(seq_len(max_value)) + updist(seq_len(max_value) + 1)
-    ) /
-      (updist(max_value) + updist(max_value + 1))
-    pmf <- diff(cmf)
-    rdist <- function(n) {
-      sample(
-        x = seq_len(max_value + 1) - 1, size = n, prob = pmf, replace = TRUE
-      )
-    }
-    pdist <- function(n) {
-      cmf[n + 1]
-    }
-    ddist <- function(n) {
-      pmf[n + 1]
-    }
-  } else {
-    pdist <- function(n) {
-      updist(n) / updist(max_value + 1)
-    }
-    ddist <- function(n) {
-      pdist(n + 1) - pdist(n)
-    }
-    if (model == "exp") {
-      rdist <- function(n) {
-        rexp(n, params[["rate"]])
-      }
-    } else if (model == "gamma") {
-      rdist <- function(n) {
-        rgamma(n, params[["shape"]], params[["rate"]])
-      }
-    } else if (model == "lognormal") {
-      rdist <- function(n) {
-        rlnorm(n, params[["meanlog"]], params[["sdlog"]])
-      }
+  ## apply tolerance if given
+  if (!missing(tolerance)) {
+    ## tolerance_max
+    tol_max <- do.call(qdist, c(list(p = 1 - tolerance), params))
+    if (missing(max_value) || tol_max < max_value) {
+      max_value <- tol_max
     }
   }
 
-  # define internal sampling function
-  inner_skel <- function(n, dist = FALSE, cum = TRUE, max_value = NULL) {
-    if (dist) {
-      if (cum) {
-        ret <- pdist(n)
-      } else {
-        ret <- ddist(n)
-      }
-      ret[ret > 1] <- NA_real_
-      return(ret)
-    } else {
-      rdist(n)
-    }
+  ## determine pmf
+  max_value <- ceiling(max_value)
+  if (max_value < width) {
+    cmf <- c(0, 1)
+  } else {
+    x <- seq(width, max_value, by = width)
+    cmf <- c(0, updist(width), (updist(x) + updist(x + width))) /
+      (updist(max_value) + updist(max_value + width))
   }
 
-  # define truncation wrapper
-  truncated_skel <- function(n, dist, cum, max_value) {
-    n <- inner_skel(n, dist, cum, max_value)
-    if (!dist) {
-      while (any(!is.na(n) & n >= max_value)) {
-        n <- ifelse(n >= max_value, inner_skel(n), n)
-      }
-
-      n <- as.integer(n)
-    }
-    return(n)
-  }
+  pmf <- diff(cmf)
 
-  # call function
-  sample <- truncated_skel(n, dist = dist, cum = cum, max_value = max_value)
-  return(sample)
+  return(pmf)
 }
 
 #' Creates a delay distribution as the sum of two other delay distributions.
@@ -407,9 +324,8 @@ max.dist_spec <- function(x, ...) {
 #' Discretise a <dist_spec>
 #'
 #' @description `r lifecycle::badge("experimental")`
-#' By default it will discretise all the distributions it can discretise
-#' (i.e. those with finite support and constant parameters).
-#' @title Discretise a <dist_spec>
+#'
+#' @inherit discrete_pmf sections references
 #' @param x A `<dist_spec>`
 #' @param strict Logical; If `TRUE` (default) an error will be thrown if a
 #' distribution cannot be discretised (e.g., because no finite maximum has been
@@ -431,9 +347,12 @@ discretise <- function(x, strict = TRUE) {
   if (!is(x, "dist_spec")) {
     stop("Can only discretise a <dist_spec>.")
   }
-  ## check max
+  tolerance <- attr(x, "tolerance")
+  if (is.null(tolerance)) {
+    tolerance <- 0
+  }
   max_x <- max(x)
-  if (any(is.infinite(max_x)) && strict) {
+  if (any(is.infinite(max_x)) && !(tolerance > 0) && strict) {
     stop("Cannot discretise a distribution with infinite support.")
   }
   ## discretise
@@ -442,13 +361,12 @@ discretise <- function(x, strict = TRUE) {
     if (y$distribution == "nonparametric") {
       return(y)
     } else {
-      if (all(vapply(y$parameters, is.numeric, logical(1))) &&
-          is.finite(max_x[id])) {
-        z <- list(pmf = dist_skel(
-          n = seq_len(max_x[id] + 1) - 1, dist = TRUE, cum = FALSE,
-          model = y$distribution, params = y$parameters,
-          max_value = max_x[id], discrete = TRUE
-        ))
+      if (all(vapply(y$parameters, is.numeric, logical(1)))) {
+        z <- list(
+          pmf = discrete_pmf(
+            y$distribution, y$parameters, y$max, tolerance, width = 1
+          )
+        )
         z$distribution <- "nonparametric"
         return(z)
       } else if (strict) {
@@ -523,7 +441,9 @@ collapse <- function(x) {
 #' This removes any part of the tail of the nonparametric distributions in the
 #' <dist_spec> where the probability mass is below the threshold level.
 #' @param x A `<dist_spec>`
-#' @param tolerance Numeric; the desired tolerance level.
+#' @param tolerance Numeric; the desired tolerance level. Any part of the
+#' cumulative distribution function beyond 1 minus this tolerance level is
+#' removed.
 #' @return A `<dist_spec>` where probability masses below the threshold level
 #' have been removed
 #' @export
diff --git a/man/apply_tolerance.Rd b/man/apply_tolerance.Rd
index b7e2c81b6..16d365443 100644
--- a/man/apply_tolerance.Rd
+++ b/man/apply_tolerance.Rd
@@ -9,7 +9,9 @@ apply_tolerance(x, tolerance)
 \arguments{
 \item{x}{A \verb{<dist_spec>}}
 
-\item{tolerance}{Numeric; the desired tolerance level.}
+\item{tolerance}{Numeric; the desired tolerance level. Any part of the
+cumulative distribution function beyond 1 minus this tolerance level is
+removed.}
 }
 \value{
 A \verb{<dist_spec>} where probability masses below the threshold level
diff --git a/man/delay_opts.Rd b/man/delay_opts.Rd
index d6e2dc381..7b0628673 100644
--- a/man/delay_opts.Rd
+++ b/man/delay_opts.Rd
@@ -20,7 +20,9 @@ a fixed distribution with all mass at 0, i.e. no delay.}
 
 \item{fixed}{deprecated; use \code{dist} instead}
 
-\item{tolerance}{Numeric; the desired tolerance level.}
+\item{tolerance}{Numeric; the desired tolerance level. Any part of the
+cumulative distribution function beyond 1 minus this tolerance level is
+removed.}
 
 \item{weight_prior}{Logical; if TRUE (default), any priors given in \code{dist}
 will be weighted by the number of observation data points, in doing so
diff --git a/man/discrete_pmf.Rd b/man/discrete_pmf.Rd
new file mode 100644
index 000000000..aa1965e05
--- /dev/null
+++ b/man/discrete_pmf.Rd
@@ -0,0 +1,61 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/dist_spec.R
+\name{discrete_pmf}
+\alias{discrete_pmf}
+\title{Discretised probability mass function}
+\usage{
+discrete_pmf(
+  distribution = c("exp", "gamma", "lognormal", "normal", "fixed"),
+  params,
+  max_value,
+  tolerance,
+  width
+)
+}
+\arguments{
+\item{distribution}{A character string representing the distribution to be
+used (one of "exp", "gamma", "lognormal", "normal" or "fixed")}
+
+\item{params}{A list of parameters values (by name) required for each model.
+For the exponential model this is a rate parameter and for the gamma model
+this is alpha and beta.}
+
+\item{max_value}{Numeric, the maximum value to allow. Defaults to 120.
+Samples outside of this range are resampled.}
+
+\item{tolerance}{Numeric; the desired tolerance level. Any part of the
+cumulative distribution function beyond 1 minus this tolerance level is
+removed.}
+
+\item{width}{Numeric, the width of each discrete bin.}
+}
+\value{
+A vector representing a probability distribution.
+}
+\description{
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#questioning}{\figure{lifecycle-questioning.svg}{options: alt='[Questioning]'}}}{\strong{[Questioning]}}
+This function returns the probability mass function of a discretised and
+truncated distribution defined by distribution type, maximum value and model
+parameters.
+}
+\section{Methodological details}{
+The probability mass function of the discretised probability distribution is
+a vector where the first entry corresponds to the integral over the (0,1]
+interval of the corresponding continuous distribution (probability of
+integer 0), the second entry corresponds to the (0,2] interval (probability
+mass of integer 1), the third entry corresponds to the (1, 3] interval
+(probability mass of integer 2), etc. This approximates the true
+probability mass function of a double censored distribution which arises
+from the difference of two censored events.
+}
+
+\references{
+Charniga, K., et al. “Best practices for estimating and reporting
+epidemiological delay distributions of infectious diseases using public
+health surveillance and healthcare data”, \emph{arXiv e-prints}, 2024.
+\url{https://doi.org/10.48550/arXiv.2405.08841}
+Park,  S. W.,  et al.,  "Estimating epidemiological delay distributions for
+infectious diseases", \emph{medRxiv}, 2024.
+\url{https://doi.org/10.1101/2024.01.12.24301247}
+}
+\keyword{internal}
diff --git a/man/discretise.Rd b/man/discretise.Rd
index dff1ddc2a..141e96e01 100644
--- a/man/discretise.Rd
+++ b/man/discretise.Rd
@@ -23,12 +23,18 @@ nonparametric.
 }
 \description{
 \ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}}
-By default it will discretise all the distributions it can discretise
-(i.e. those with finite support and constant parameters).
 }
-\details{
-Discretise a <dist_spec>
+\section{Methodological details}{
+The probability mass function of the discretised probability distribution is
+a vector where the first entry corresponds to the integral over the (0,1]
+interval of the corresponding continuous distribution (probability of
+integer 0), the second entry corresponds to the (0,2] interval (probability
+mass of integer 1), the third entry corresponds to the (1, 3] interval
+(probability mass of integer 2), etc. This approximates the true
+probability mass function of a double censored distribution which arises
+from the difference of two censored events.
 }
+
 \examples{
 # A fixed gamma distribution with mean 5 and sd 1.
 dist1 <- Gamma(mean = 5, sd = 1, max = 20)
@@ -39,3 +45,12 @@ dist2 <- LogNormal(mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20)
 # The maxf the sum of two distributions
 discretise(dist1 + dist2, strict = FALSE)
 }
+\references{
+Charniga, K., et al. “Best practices for estimating and reporting
+epidemiological delay distributions of infectious diseases using public
+health surveillance and healthcare data”, \emph{arXiv e-prints}, 2024.
+\url{https://doi.org/10.48550/arXiv.2405.08841}
+Park,  S. W.,  et al.,  "Estimating epidemiological delay distributions for
+infectious diseases", \emph{medRxiv}, 2024.
+\url{https://doi.org/10.1101/2024.01.12.24301247}
+}
diff --git a/man/dist_skel.Rd b/man/dist_skel.Rd
index 2ce7f40f3..bb589b227 100644
--- a/man/dist_skel.Rd
+++ b/man/dist_skel.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/dist_spec.R
+% Please edit documentation in R/deprecated.R
 \name{dist_skel}
 \alias{dist_skel}
 \title{Distribution Skeleton}
@@ -48,7 +48,7 @@ Samples outside of this range are resampled.}
 A vector of samples or a probability distribution.
 }
 \description{
-\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#questioning}{\figure{lifecycle-questioning.svg}{options: alt='[Questioning]'}}}{\strong{[Questioning]}}
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[Deprecated]}}
 This function acts as a skeleton for a truncated distribution defined by
 model type, maximum value and model parameters. It is designed to be used
 with the output from \code{\link[=get_dist]{get_dist()}}.
diff --git a/man/generation_time_opts.Rd b/man/generation_time_opts.Rd
index 5cf19c2f1..89d555113 100644
--- a/man/generation_time_opts.Rd
+++ b/man/generation_time_opts.Rd
@@ -41,7 +41,9 @@ distribution is given a fixed generation time of 1 will be assumed.}
 
 \item{fixed}{deprecated; use \code{dist} instead}
 
-\item{tolerance}{Numeric; the desired tolerance level.}
+\item{tolerance}{Numeric; the desired tolerance level. Any part of the
+cumulative distribution function beyond 1 minus this tolerance level is
+removed.}
 
 \item{weight_prior}{Logical; if TRUE (default), any priors given in \code{dist}
 will be weighted by the number of observation data points, in doing so
diff --git a/man/trunc_opts.Rd b/man/trunc_opts.Rd
index 7205e3566..404f45c98 100644
--- a/man/trunc_opts.Rd
+++ b/man/trunc_opts.Rd
@@ -14,7 +14,9 @@ interface in \code{EpiNow2} (See \code{?EpiNow2::Distributions}) or estimated us
 for use here out-of-box. Default is a fixed distribution with maximum 0, i.e.
 no truncation.}
 
-\item{tolerance}{Numeric; the desired tolerance level.}
+\item{tolerance}{Numeric; the desired tolerance level. Any part of the
+cumulative distribution function beyond 1 minus this tolerance level is
+removed.}
 
 \item{weight_prior}{Logical; if TRUE, the truncation prior will be weighted
 by the number of observation data points, in doing so approximately placing