diff --git a/NEWS.md b/NEWS.md
index 93cdeb875..7e6a5030f 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -30,6 +30,7 @@
 
 - Brought the docs on `alpha_sd` up to date with the code change from prior PR #853. By @zsusswein in #862 and reviewed by @jamesmbaazam.
 - The `...` argument in `estimate_secondary()` has been removed because it was not used. By @jamesmbaazam in #894 and reviewed by @.
+- All examples now use the natural parameters of distributions rather than the mean and standard deviation when specifying uncertain distributions. This is to eliminate warnings and encourage best practice. By @jamesmbaazam in #893 and reviewed by @sbfnk.
 
 # EpiNow2 1.6.1
 
diff --git a/R/dist_spec.R b/R/dist_spec.R
index 5b9a155ba..cce218198 100644
--- a/R/dist_spec.R
+++ b/R/dist_spec.R
@@ -116,9 +116,12 @@ discrete_pmf <- function(distribution =
 #' )
 #' dist1 + dist1
 #'
-#' # An uncertain gamma distribution with mean 3 and sd 2
+#' # An uncertain gamma distribution with shape and rate normally distributed
+#' # as Normal(3, 0.5) and Normal(2, 0.5) respectively
 #' dist2 <- Gamma(
-#'   mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20
+#'   shape = Normal(3, 0.5),
+#'   rate = Normal(2, 0.5),
+#'   max = 20
 #' )
 #' dist1 + dist2
 `+.dist_spec` <- function(e1, e2) {
@@ -197,9 +200,12 @@ discrete_pmf <- function(distribution =
 #' )
 #' dist1 + dist1
 #'
-#' # An uncertain gamma distribution with mean 3 and sd 2
+#' # An uncertain gamma distribution with shape and rate normally distributed
+#' # as Normal(3, 0.5) and Normal(2, 0.5) respectively
 #' dist2 <- Gamma(
-#'   mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20
+#'   shape = Normal(3, 0.5),
+#'   rate = Normal(2, 0.5),
+#'   max = 20
 #' )
 #' c(dist1, dist2)
 c.dist_spec <- function(...) {
@@ -257,9 +263,12 @@ c.dist_spec <- function(...) {
 #' dist1 <- LogNormal(mean = 5, sd = 1, max = 20)
 #' mean(dist1)
 #'
-#' # An uncertain gamma distribution with mean 3 and sd 2
+#' # An uncertain gamma distribution with shape and rate normally distributed
+#' # as Normal(3, 0.5) and Normal(2, 0.5) respectively
 #' dist2 <- Gamma(
-#'  mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20
+#'   shape = Normal(3, 0.5),
+#'   rate = Normal(2, 0.5),
+#'   max = 20
 #' )
 #' mean(dist2)
 #'
@@ -393,8 +402,13 @@ sd.default <- function(x, ...) {
 #' dist1 <- Gamma(mean = 5, sd = 1, max = 20)
 #' max(dist1)
 #'
-#' # An uncertain lognormal distribution with mean 3 and sd 2
-#' dist2 <- LogNormal(mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20)
+#' # An uncertain lognormal distribution with meanlog and sdlog normally
+#' # distributed as Normal(3, 0.5) and Normal(2, 0.5) respectively
+#' dist2 <- LogNormal(
+#'   meanlog = Normal(3, 0.5),
+#'   sdlog = Normal(2, 0.5),
+#'   max = 20
+#' )
 #' max(dist2)
 #'
 #' # The max the sum of two distributions
@@ -447,8 +461,13 @@ discretise <- function(x, ...) {
 #' # A fixed gamma distribution with mean 5 and sd 1.
 #' dist1 <- Gamma(mean = 5, sd = 1, max = 20)
 #'
-#' # An uncertain lognormal distribution with mean 3 and sd 2
-#' dist2 <- LogNormal(mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20)
+#' # An uncertain lognormal distribution with meanlog and sdlog normally
+#' # distributed as Normal(3, 0.5) and Normal(2, 0.5) respectively
+#' dist2 <- LogNormal(
+#'   meanlog = Normal(3, 0.5),
+#'   sdlog = Normal(2, 0.5),
+#'   max = 20
+#' )
 #'
 #' # The maxf the sum of two distributions
 #' discretise(dist1 + dist2, strict = FALSE)
@@ -531,11 +550,16 @@ collapse <- function(x, ...) {
 #' # A fixed gamma distribution with mean 5 and sd 1.
 #' dist1 <- Gamma(mean = 5, sd = 1, max = 20)
 #'
-#' # An uncertain lognormal distribution with mean 3 and sd 2
-#' dist2 <- LogNormal(mean = 3, sd = 2, max = 20)
+#' # An uncertain lognormal distribution with meanlog and sdlog normally
+#' # distributed as Normal(3, 0.5) and Normal(2, 0.5) respectively
+#' dist2 <- LogNormal(
+#'   meanlog = Normal(3, 0.5),
+#'   sdlog = Normal(2, 0.5),
+#'   max = 20
+#' )
 #'
 #' # The maxf the sum of two distributions
-#' collapse(discretise(dist1 + dist2))
+#' collapse(discretise(dist1 + dist2, strict = FALSE))
 collapse.dist_spec <- function(x, ...) {
   return(x)
 }
@@ -588,9 +612,10 @@ collapse.multi_dist_spec <- function(x, ...) {
 #' dist1 <- LogNormal(mean = 1.5, sd = 0.5, max = 20)
 #' print(dist1)
 #'
-#' # An uncertain gamma distribution with mean 3 and sd 2
+#' # An uncertain gamma distribution with shape and rate normally distributed
+#' # as Normal(3, 0.5) and Normal(2, 0.5) respectively
 #' dist2 <- Gamma(
-#'   mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20
+#'   shape = Normal(3, 0.5), rate = Normal(2, 0.5), max = 20
 #' )
 #' print(dist2)
 print.dist_spec <- function(x, ...) {
@@ -680,9 +705,12 @@ print.dist_spec <- function(x, ...) {
 #' # Plot discretised distribution with 0.01 day discretisation window
 #' plot(dist1, res = 0.01, cumulative = FALSE)
 #'
-#' # An uncertain gamma distribution with mean 3 and sd 2
+#' # An uncertain gamma distribution with shape and rate normally distributed
+#' # as Normal(3, 0.5) and Normal(2, 0.5) respectively
 #' dist2 <- Gamma(
-#'   mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20
+#'   shape = Normal(3, 0.5),
+#'   rate = Normal(2, 0.5),
+#'   max = 20
 #' )
 #' plot(dist2)
 #'
@@ -781,9 +809,12 @@ plot.dist_spec <- function(x, samples = 50L, res = 1, cumulative = TRUE, ...) {
 #' @examples
 #' dist1 <- LogNormal(mean = 1.6, sd = 0.5, max = 20)
 #'
-#' # An uncertain gamma distribution with mean 3 and sd 2
+#' # An uncertain gamma distribution with shape and rate normally distributed
+#' # as Normal(3, 0.5) and Normal(2, 0.5) respectively
 #' dist2 <- Gamma(
-#'   mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20
+#'   shape = Normal(3, 0.5),
+#'   rate = Normal(2, 0.5),
+#'   max = 20
 #' )
 #'
 #' # Multiple distributions
@@ -828,9 +859,12 @@ fix_parameters <- function(x, ...) {
 #' @importFrom rlang arg_match
 #' @method fix_parameters dist_spec
 #' @examples
-#' # An uncertain gamma distribution with mean 3 and sd 2
-#' dist <- LogNormal(
-#'   meanlog = Normal(3, 0.5), sdlog = Normal(2, 0.5), max = 20
+#' # An uncertain gamma distribution with shape and rate normally distributed
+#' # as Normal(3, 0.5) and Normal(2, 0.5) respectively
+#' dist <- Gamma(
+#'   shape = Normal(3, 0.5),
+#'   rate = Normal(2, 0.5),
+#'   max = 20
 #' )
 #'
 #' fix_parameters(dist)
@@ -891,8 +925,13 @@ is_constrained <- function(x, ...) {
 #' # A fixed gamma distribution with mean 5 and sd 1.
 #' dist1 <- Gamma(mean = 5, sd = 1, max = 20)
 #'
-#' # An uncertain lognormal distribution with mean 3 and sd 2
-#' dist2 <- LogNormal(mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20)
+#' # An uncertain lognormal distribution with meanlog and sdlog normally
+#' # distributed as Normal(3, 0.5) and Normal(2, 0.5) respectively
+#' dist2 <- LogNormal(
+#'   meanlog = Normal(3, 0.5),
+#'   sdlog = Normal(2, 0.5),
+#'   max = 20
+#' )
 #'
 #' # both distributions are constrained and therefore so is the sum
 #' is_constrained(dist1 + dist2)
@@ -952,7 +991,8 @@ is_constrained.multi_dist_spec <- function(x, ...) {
 #' @examples
 #' LogNormal(mean = 4, sd = 1)
 #' LogNormal(mean = 4, sd = 1, max = 10)
-#' LogNormal(mean = Normal(4, 1), sd = 1, max = 10)
+#' # If specifying uncertain parameters, use the natural parameters
+#' LogNormal(meanlog = Normal(1.5, 0.5), sdlog = 0.25, max = 10)
 LogNormal <- function(meanlog, sdlog, mean, sd, ...) {
   params <- as.list(environment())
   return(new_dist_spec(params, "lognormal", ...))
diff --git a/R/opts.R b/R/opts.R
index a791ce882..b1aacd852 100644
--- a/R/opts.R
+++ b/R/opts.R
@@ -42,8 +42,8 @@
 #' # An uncertain gamma distributed generation time
 #' generation_time_opts(
 #'   Gamma(
-#'     mean = Normal(mean = 3, sd = 1),
-#'     sd = Normal(mean = 2, sd = 0.5),
+#'     shape = Normal(mean = 3, sd = 1),
+#'     rate = Normal(mean = 2, sd = 0.5),
 #'     max = 14
 #'   )
 #' )
@@ -186,7 +186,11 @@ secondary_opts <- function(type = c("incidence", "prevalence"), ...) {
 #' delay_opts()
 #'
 #' # A single delay that has uncertainty
-#' delay <- LogNormal(mean = Normal(1, 0.2), sd = Normal(0.5, 0.1), max = 14)
+#' delay <- LogNormal(
+#'   meanlog = Normal(1, 0.2),
+#'   sdlog = Normal(0.5, 0.1),
+#'   max = 14
+#' )
 #' delay_opts(delay)
 #'
 #' # A single delay without uncertainty
diff --git a/man/Distributions.Rd b/man/Distributions.Rd
index 583b0193d..add503e54 100644
--- a/man/Distributions.Rd
+++ b/man/Distributions.Rd
@@ -77,7 +77,8 @@ probability mass function is given directly as a numeric vector.
 \examples{
 LogNormal(mean = 4, sd = 1)
 LogNormal(mean = 4, sd = 1, max = 10)
-LogNormal(mean = Normal(4, 1), sd = 1, max = 10)
+# If specifying uncertain parameters, use the natural parameters
+LogNormal(meanlog = Normal(1.5, 0.5), sdlog = 0.25, max = 10)
 Gamma(mean = 4, sd = 1)
 Gamma(shape = 16, rate = 4)
 Gamma(shape = Normal(16, 2), rate = Normal(4, 1))
diff --git a/man/c.dist_spec.Rd b/man/c.dist_spec.Rd
index 45cffbf55..ae603e9e2 100644
--- a/man/c.dist_spec.Rd
+++ b/man/c.dist_spec.Rd
@@ -27,9 +27,12 @@ dist1 <- LogNormal(
 )
 dist1 + dist1
 
-# An uncertain gamma distribution with mean 3 and sd 2
+# An uncertain gamma distribution with shape and rate normally distributed
+# as Normal(3, 0.5) and Normal(2, 0.5) respectively
 dist2 <- Gamma(
-  mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20
+  shape = Normal(3, 0.5),
+  rate = Normal(2, 0.5),
+  max = 20
 )
 c(dist1, dist2)
 }
diff --git a/man/collapse.Rd b/man/collapse.Rd
index 97ed0d1c2..9715d9c33 100644
--- a/man/collapse.Rd
+++ b/man/collapse.Rd
@@ -25,9 +25,14 @@ in the <dist_spec>.
 # A fixed gamma distribution with mean 5 and sd 1.
 dist1 <- Gamma(mean = 5, sd = 1, max = 20)
 
-# An uncertain lognormal distribution with mean 3 and sd 2
-dist2 <- LogNormal(mean = 3, sd = 2, max = 20)
+# An uncertain lognormal distribution with meanlog and sdlog normally
+# distributed as Normal(3, 0.5) and Normal(2, 0.5) respectively
+dist2 <- LogNormal(
+  meanlog = Normal(3, 0.5),
+  sdlog = Normal(2, 0.5),
+  max = 20
+)
 
 # The maxf the sum of two distributions
-collapse(discretise(dist1 + dist2))
+collapse(discretise(dist1 + dist2, strict = FALSE))
 }
diff --git a/man/delay_opts.Rd b/man/delay_opts.Rd
index b28a47b72..dc383a541 100644
--- a/man/delay_opts.Rd
+++ b/man/delay_opts.Rd
@@ -44,7 +44,11 @@ functions.
 delay_opts()
 
 # A single delay that has uncertainty
-delay <- LogNormal(mean = Normal(1, 0.2), sd = Normal(0.5, 0.1), max = 14)
+delay <- LogNormal(
+  meanlog = Normal(1, 0.2),
+  sdlog = Normal(0.5, 0.1),
+  max = 14
+)
 delay_opts(delay)
 
 # A single delay without uncertainty
diff --git a/man/discretise.Rd b/man/discretise.Rd
index e6c2cda1a..a76646d4a 100644
--- a/man/discretise.Rd
+++ b/man/discretise.Rd
@@ -42,8 +42,13 @@ from the difference of two censored events.
 # A fixed gamma distribution with mean 5 and sd 1.
 dist1 <- Gamma(mean = 5, sd = 1, max = 20)
 
-# An uncertain lognormal distribution with mean 3 and sd 2
-dist2 <- LogNormal(mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20)
+# An uncertain lognormal distribution with meanlog and sdlog normally
+# distributed as Normal(3, 0.5) and Normal(2, 0.5) respectively
+dist2 <- LogNormal(
+  meanlog = Normal(3, 0.5),
+  sdlog = Normal(2, 0.5),
+  max = 20
+)
 
 # The maxf the sum of two distributions
 discretise(dist1 + dist2, strict = FALSE)
diff --git a/man/extract_single_dist.Rd b/man/extract_single_dist.Rd
index e3407cf6b..e6761e96c 100644
--- a/man/extract_single_dist.Rd
+++ b/man/extract_single_dist.Rd
@@ -20,9 +20,12 @@ A single \code{dist_spec} object
 \examples{
 dist1 <- LogNormal(mean = 1.6, sd = 0.5, max = 20)
 
-# An uncertain gamma distribution with mean 3 and sd 2
+# An uncertain gamma distribution with shape and rate normally distributed
+# as Normal(3, 0.5) and Normal(2, 0.5) respectively
 dist2 <- Gamma(
-  mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20
+  shape = Normal(3, 0.5),
+  rate = Normal(2, 0.5),
+  max = 20
 )
 
 # Multiple distributions
diff --git a/man/fix_parameters.Rd b/man/fix_parameters.Rd
index 04a8a4136..a43df5bbd 100644
--- a/man/fix_parameters.Rd
+++ b/man/fix_parameters.Rd
@@ -25,9 +25,12 @@ If the given \verb{<dist_spec>} has any uncertainty, it is removed and the
 corresponding distribution converted into a fixed one.
 }
 \examples{
-# An uncertain gamma distribution with mean 3 and sd 2
-dist <- LogNormal(
-  meanlog = Normal(3, 0.5), sdlog = Normal(2, 0.5), max = 20
+# An uncertain gamma distribution with shape and rate normally distributed
+# as Normal(3, 0.5) and Normal(2, 0.5) respectively
+dist <- Gamma(
+  shape = Normal(3, 0.5),
+  rate = Normal(2, 0.5),
+  max = 20
 )
 
 fix_parameters(dist)
diff --git a/man/generation_time_opts.Rd b/man/generation_time_opts.Rd
index 1200101c7..8306d7a61 100644
--- a/man/generation_time_opts.Rd
+++ b/man/generation_time_opts.Rd
@@ -79,8 +79,8 @@ generation_time_opts(Gamma(mean = 3, sd = 2, max = 14))
 # An uncertain gamma distributed generation time
 generation_time_opts(
   Gamma(
-    mean = Normal(mean = 3, sd = 1),
-    sd = Normal(mean = 2, sd = 0.5),
+    shape = Normal(mean = 3, sd = 1),
+    rate = Normal(mean = 2, sd = 0.5),
     max = 14
   )
 )
diff --git a/man/is_constrained.Rd b/man/is_constrained.Rd
index 7a097ab6d..39803ee62 100644
--- a/man/is_constrained.Rd
+++ b/man/is_constrained.Rd
@@ -23,8 +23,13 @@ Logical; TRUE if \code{x} is constrained
 # A fixed gamma distribution with mean 5 and sd 1.
 dist1 <- Gamma(mean = 5, sd = 1, max = 20)
 
-# An uncertain lognormal distribution with mean 3 and sd 2
-dist2 <- LogNormal(mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20)
+# An uncertain lognormal distribution with meanlog and sdlog normally
+# distributed as Normal(3, 0.5) and Normal(2, 0.5) respectively
+dist2 <- LogNormal(
+  meanlog = Normal(3, 0.5),
+  sdlog = Normal(2, 0.5),
+  max = 20
+)
 
 # both distributions are constrained and therefore so is the sum
 is_constrained(dist1 + dist2)
diff --git a/man/max.dist_spec.Rd b/man/max.dist_spec.Rd
index c81f811ef..516a0cc3a 100644
--- a/man/max.dist_spec.Rd
+++ b/man/max.dist_spec.Rd
@@ -25,8 +25,13 @@ in parameters)
 dist1 <- Gamma(mean = 5, sd = 1, max = 20)
 max(dist1)
 
-# An uncertain lognormal distribution with mean 3 and sd 2
-dist2 <- LogNormal(mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20)
+# An uncertain lognormal distribution with meanlog and sdlog normally
+# distributed as Normal(3, 0.5) and Normal(2, 0.5) respectively
+dist2 <- LogNormal(
+  meanlog = Normal(3, 0.5),
+  sdlog = Normal(2, 0.5),
+  max = 20
+)
 max(dist2)
 
 # The max the sum of two distributions
diff --git a/man/mean.dist_spec.Rd b/man/mean.dist_spec.Rd
index 5448ca1b8..8be98f2d0 100644
--- a/man/mean.dist_spec.Rd
+++ b/man/mean.dist_spec.Rd
@@ -25,9 +25,12 @@ distributions combined in the passed <dist_spec>.
 dist1 <- LogNormal(mean = 5, sd = 1, max = 20)
 mean(dist1)
 
-# An uncertain gamma distribution with mean 3 and sd 2
+# An uncertain gamma distribution with shape and rate normally distributed
+# as Normal(3, 0.5) and Normal(2, 0.5) respectively
 dist2 <- Gamma(
- mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20
+  shape = Normal(3, 0.5),
+  rate = Normal(2, 0.5),
+  max = 20
 )
 mean(dist2)
 
diff --git a/man/plot.dist_spec.Rd b/man/plot.dist_spec.Rd
index 27698649a..7e1a8d19e 100644
--- a/man/plot.dist_spec.Rd
+++ b/man/plot.dist_spec.Rd
@@ -33,9 +33,12 @@ plot(dist1)
 # Plot discretised distribution with 0.01 day discretisation window
 plot(dist1, res = 0.01, cumulative = FALSE)
 
-# An uncertain gamma distribution with mean 3 and sd 2
+# An uncertain gamma distribution with shape and rate normally distributed
+# as Normal(3, 0.5) and Normal(2, 0.5) respectively
 dist2 <- Gamma(
-  mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20
+  shape = Normal(3, 0.5),
+  rate = Normal(2, 0.5),
+  max = 20
 )
 plot(dist2)
 
diff --git a/man/plus-.dist_spec.Rd b/man/plus-.dist_spec.Rd
index 16e967faa..091b7e536 100644
--- a/man/plus-.dist_spec.Rd
+++ b/man/plus-.dist_spec.Rd
@@ -26,9 +26,12 @@ dist1 <- LogNormal(
 )
 dist1 + dist1
 
-# An uncertain gamma distribution with mean 3 and sd 2
+# An uncertain gamma distribution with shape and rate normally distributed
+# as Normal(3, 0.5) and Normal(2, 0.5) respectively
 dist2 <- Gamma(
-  mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20
+  shape = Normal(3, 0.5),
+  rate = Normal(2, 0.5),
+  max = 20
 )
 dist1 + dist2
 }
diff --git a/man/print.dist_spec.Rd b/man/print.dist_spec.Rd
index c3ea43cfb..b6925af9f 100644
--- a/man/print.dist_spec.Rd
+++ b/man/print.dist_spec.Rd
@@ -24,9 +24,10 @@ functions of fixed delay distributions combined in the passed <dist_spec>.
 dist1 <- LogNormal(mean = 1.5, sd = 0.5, max = 20)
 print(dist1)
 
-# An uncertain gamma distribution with mean 3 and sd 2
+# An uncertain gamma distribution with shape and rate normally distributed
+# as Normal(3, 0.5) and Normal(2, 0.5) respectively
 dist2 <- Gamma(
-  mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20
+  shape = Normal(3, 0.5), rate = Normal(2, 0.5), max = 20
 )
 print(dist2)
 }
diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R
index d6614247f..697d991f9 100644
--- a/tests/testthat/setup.R
+++ b/tests/testthat/setup.R
@@ -16,4 +16,6 @@ if (identical(Sys.getenv("NOT_CRAN"), "true")) {
   }
 }
 
-withr::defer(future::plan("sequential"), teardown_env())
+if (requireNamespace("future", quietly = TRUE)) {
+  withr::defer(future::plan("sequential"), teardown_env())
+}
diff --git a/tests/testthat/test-setup_future.R b/tests/testthat/test-setup_future.R
index aeb1b522f..a86011b9f 100644
--- a/tests/testthat/test-setup_future.R
+++ b/tests/testthat/test-setup_future.R
@@ -1,3 +1,4 @@
+skip_on_cran()
 reported_cases <- data.frame(region = c("test", "boo"))
 futile.logger::flog.threshold("FATAL")
 
diff --git a/vignettes/estimate_infections_options-bp-1.png b/vignettes/estimate_infections_options-bp-1.png
index 68fc9157c..c9321d83e 100644
Binary files a/vignettes/estimate_infections_options-bp-1.png and b/vignettes/estimate_infections_options-bp-1.png differ
diff --git a/vignettes/estimate_infections_options-fixed-1.png b/vignettes/estimate_infections_options-fixed-1.png
index c56ed5dc9..8e82ea94e 100644
Binary files a/vignettes/estimate_infections_options-fixed-1.png and b/vignettes/estimate_infections_options-fixed-1.png differ
diff --git a/vignettes/estimate_infections_options-gp_projection-1.png b/vignettes/estimate_infections_options-gp_projection-1.png
index cf18cf4ce..2137d64b6 100644
Binary files a/vignettes/estimate_infections_options-gp_projection-1.png and b/vignettes/estimate_infections_options-gp_projection-1.png differ
diff --git a/vignettes/estimate_infections_options-no_delays-1.png b/vignettes/estimate_infections_options-no_delays-1.png
index 2081b6b3b..d3792e1d0 100644
Binary files a/vignettes/estimate_infections_options-no_delays-1.png and b/vignettes/estimate_infections_options-no_delays-1.png differ
diff --git a/vignettes/estimate_infections_options-nonparametric-1.png b/vignettes/estimate_infections_options-nonparametric-1.png
index 31b938af6..9586e7a32 100644
Binary files a/vignettes/estimate_infections_options-nonparametric-1.png and b/vignettes/estimate_infections_options-nonparametric-1.png differ
diff --git a/vignettes/estimate_infections_options-truncation-1.png b/vignettes/estimate_infections_options-truncation-1.png
index 1f6c31cf6..a1848d3b4 100644
Binary files a/vignettes/estimate_infections_options-truncation-1.png and b/vignettes/estimate_infections_options-truncation-1.png differ
diff --git a/vignettes/estimate_infections_options-weekly_rw-1.png b/vignettes/estimate_infections_options-weekly_rw-1.png
index bd7c7a6c0..069381639 100644
Binary files a/vignettes/estimate_infections_options-weekly_rw-1.png and b/vignettes/estimate_infections_options-weekly_rw-1.png differ
diff --git a/vignettes/estimate_infections_options.Rmd b/vignettes/estimate_infections_options.Rmd
index 28cabeeba..748760b2f 100644
--- a/vignettes/estimate_infections_options.Rmd
+++ b/vignettes/estimate_infections_options.Rmd
@@ -198,10 +198,10 @@ summary(def)
 # elapsed time (in seconds)
 get_elapsed_time(def$fit)
 #>         warmup sample
-#> chain:1 23.301 18.698
-#> chain:2 29.287 28.181
-#> chain:3 37.962 17.803
-#> chain:4 37.761 14.090
+#> chain:1 23.529 18.991
+#> chain:2 29.477 28.414
+#> chain:3 37.718 18.235
+#> chain:4 38.292 14.125
 # summary plot
 plot(def)
 ```
@@ -232,10 +232,10 @@ summary(agp)
 # elapsed time (in seconds)
 get_elapsed_time(agp$fit)
 #>         warmup sample
-#> chain:1 27.748 19.571
-#> chain:2 23.701 18.866
-#> chain:3 22.463 28.347
-#> chain:4 26.436 19.798
+#> chain:1 28.106 18.768
+#> chain:2 24.271 19.359
+#> chain:3 22.602 30.107
+#> chain:4 26.954 20.210
 # summary plot
 plot(agp)
 ```
@@ -270,10 +270,10 @@ summary(dep)
 # elapsed time (in seconds)
 get_elapsed_time(dep$fit)
 #>         warmup sample
-#> chain:1 42.083 23.377
-#> chain:2 31.143 20.577
-#> chain:3 39.918 19.875
-#> chain:4 38.332 20.637
+#> chain:1 42.625 23.562
+#> chain:2 31.544 20.825
+#> chain:3 40.299 20.096
+#> chain:4 38.885 20.701
 # summary plot
 plot(dep)
 ```
@@ -288,30 +288,16 @@ Here, instead of doing so we assume that we know about truncation with mean of 1
 
 ``` r
 trunc_dist <- LogNormal(
-  mean = Normal(0.5, 0.1),
-  sd = Normal(0.5, 0.1),
+  mean = 0.5,
+  sd = 0.5,
   max = 3
 )
-#> Warning: ! Uncertain lognormal distribution specified in terms of parameters that are
-#>   not the "natural" parameters of the distribution meanlog and sdlog.
-#> ℹ Converting using a crude and very approximate method that is likely to
-#>   produce biased results.
-#> ℹ If possible it is preferable to specify the distribution directly in terms of
-#>   the natural parameters.
 trunc_dist
 #> - lognormal distribution (max: 3):
 #>   meanlog:
-#>     - normal distribution:
-#>       mean:
-#>         -1
-#>       sd:
-#>         0.14
+#>     -1
 #>   sdlog:
-#>     - normal distribution:
-#>       mean:
-#>         0.83
-#>       sd:
-#>         0.13
+#>     0.83
 ```
 
 We can then use this in the `esimtate_infections()` function using the `truncation` option.
@@ -324,14 +310,22 @@ trunc <- estimate_infections(reported_cases,
   truncation = trunc_opts(trunc_dist),
   rt = rt_opts(prior = rt_prior)
 )
-#> Error in vapply(delays[parametric], attr, "weight_prior", FUN.VALUE = logical(1)): values must be length 1,
-#>  but FUN(X[[3]]) result is length 0
 # summarise results
 summary(trunc)
-#> Error in object[[i]]: object of type 'builtin' is not subsettable
+#>                             measure               estimate
+#>                              <char>                 <char>
+#> 1:           New infections per day    4273 (2401 -- 8111)
+#> 2: Expected change in daily reports      Likely increasing
+#> 3:       Effective reproduction no.      1.1 (0.87 -- 1.4)
+#> 4:                   Rate of growth 0.032 (-0.043 -- 0.13)
+#> 5:     Doubling/halving time (days)        21 (5.5 -- -16)
 # elapsed time (in seconds)
 get_elapsed_time(trunc$fit)
-#> Error in (function (cond) : error in evaluating the argument 'object' in selecting a method for function 'get_elapsed_time': object of type 'builtin' is not subsettable
+#>         warmup sample
+#> chain:1 32.460 19.479
+#> chain:2 36.538 18.114
+#> chain:3 24.322 31.656
+#> chain:4 30.328 20.364
 # summary plot
 plot(trunc)
 ```
@@ -356,18 +350,18 @@ project_rt <- estimate_infections(reported_cases,
 summary(project_rt)
 #>                             measure                 estimate
 #>                              <char>                   <char>
-#> 1:           New infections per day      2265 (1291 -- 4223)
+#> 1:           New infections per day      2277 (1300 -- 4159)
 #> 2: Expected change in daily reports        Likely decreasing
-#> 3:       Effective reproduction no.        0.9 (0.71 -- 1.2)
-#> 4:                   Rate of growth -0.026 (-0.099 -- 0.057)
-#> 5:     Doubling/halving time (days)           -27 (12 -- -7)
+#> 3:       Effective reproduction no.        0.9 (0.71 -- 1.1)
+#> 4:                   Rate of growth -0.028 (-0.098 -- 0.052)
+#> 5:     Doubling/halving time (days)           -25 (13 -- -7)
 # elapsed time (in seconds)
 get_elapsed_time(project_rt$fit)
 #>         warmup sample
-#> chain:1 32.436 20.433
-#> chain:2 30.657 30.856
-#> chain:3 38.811 29.771
-#> chain:4 35.218 21.061
+#> chain:1 41.792 18.973
+#> chain:2 34.302 20.610
+#> chain:3 38.729 18.627
+#> chain:4 33.687 20.569
 # summary plot
 plot(project_rt)
 ```
@@ -389,18 +383,18 @@ fixed <- estimate_infections(reported_cases,
 summary(fixed)
 #>                             measure               estimate
 #>                              <char>                 <char>
-#> 1:           New infections per day  15985 (9499 -- 27668)
+#> 1:           New infections per day  15719 (9135 -- 27752)
 #> 2: Expected change in daily reports             Increasing
 #> 3:       Effective reproduction no.       1.2 (1.1 -- 1.3)
-#> 4:                   Rate of growth 0.048 (0.034 -- 0.063)
-#> 5:     Doubling/halving time (days)          14 (11 -- 20)
+#> 4:                   Rate of growth 0.048 (0.033 -- 0.062)
+#> 5:     Doubling/halving time (days)          14 (11 -- 21)
 # elapsed time (in seconds)
 get_elapsed_time(fixed$fit)
 #>         warmup sample
-#> chain:1  2.584  1.781
-#> chain:2  2.518  1.999
-#> chain:3  2.782  1.998
-#> chain:4  3.113  1.698
+#> chain:1  3.135  2.078
+#> chain:2  3.071  1.926
+#> chain:3  2.708  1.865
+#> chain:4  3.069  1.897
 # summary plot
 plot(fixed)
 ```
@@ -435,18 +429,18 @@ bkp <- estimate_infections(bp_cases,
 summary(bkp)
 #>                             measure                 estimate
 #>                              <char>                   <char>
-#> 1:           New infections per day      2343 (1927 -- 2880)
+#> 1:           New infections per day      2363 (1947 -- 2891)
 #> 2: Expected change in daily reports               Decreasing
 #> 3:       Effective reproduction no.      0.89 (0.86 -- 0.92)
-#> 4:                   Rate of growth -0.027 (-0.035 -- -0.02)
-#> 5:     Doubling/halving time (days)         -25 (-35 -- -20)
+#> 4:                   Rate of growth -0.027 (-0.034 -- -0.02)
+#> 5:     Doubling/halving time (days)         -26 (-35 -- -20)
 # elapsed time (in seconds)
 get_elapsed_time(bkp$fit)
 #>         warmup sample
-#> chain:1  4.169  4.144
-#> chain:2  4.123  4.556
-#> chain:3  4.067  4.346
-#> chain:4  3.840  4.380
+#> chain:1  4.777  4.458
+#> chain:2  4.579  4.862
+#> chain:3  4.953  5.284
+#> chain:4  4.385  5.380
 # summary plot
 plot(bkp)
 ```
@@ -470,18 +464,18 @@ rw <- estimate_infections(reported_cases,
 summary(rw)
 #>                             measure                estimate
 #>                              <char>                  <char>
-#> 1:           New infections per day     2139 (1069 -- 4280)
+#> 1:           New infections per day     2079 (1032 -- 4439)
 #> 2: Expected change in daily reports       Likely decreasing
-#> 3:       Effective reproduction no.      0.87 (0.63 -- 1.2)
-#> 4:                   Rate of growth -0.036 (-0.11 -- 0.044)
-#> 5:     Doubling/halving time (days)        -19 (16 -- -6.4)
+#> 3:       Effective reproduction no.      0.86 (0.62 -- 1.2)
+#> 4:                   Rate of growth -0.038 (-0.12 -- 0.047)
+#> 5:     Doubling/halving time (days)          -18 (15 -- -6)
 # elapsed time (in seconds)
 get_elapsed_time(rw$fit)
 #>         warmup sample
-#> chain:1  9.343 13.485
-#> chain:2  8.843 15.550
-#> chain:3  8.556 10.431
-#> chain:4  9.442 13.621
+#> chain:1  9.655 15.971
+#> chain:2  9.614 16.314
+#> chain:3  9.230 12.979
+#> chain:4 10.482 16.692
 # summary plot
 plot(rw)
 ```
@@ -502,18 +496,18 @@ no_delay <- estimate_infections(
 summary(no_delay)
 #>                             measure                  estimate
 #>                              <char>                    <char>
-#> 1:           New infections per day       2804 (2403 -- 3257)
+#> 1:           New infections per day       2801 (2409 -- 3289)
 #> 2: Expected change in daily reports                Decreasing
-#> 3:       Effective reproduction no.       0.89 (0.81 -- 0.97)
-#> 4:                   Rate of growth -0.031 (-0.062 -- 0.0012)
-#> 5:     Doubling/halving time (days)          -22 (570 -- -11)
+#> 3:       Effective reproduction no.        0.89 (0.8 -- 0.98)
+#> 4:                   Rate of growth -0.03 (-0.064 -- 0.00074)
+#> 5:     Doubling/halving time (days)          -23 (940 -- -11)
 # elapsed time (in seconds)
 get_elapsed_time(no_delay$fit)
 #>         warmup sample
-#> chain:1 33.675 32.734
-#> chain:2 43.810 35.023
-#> chain:3 47.974 42.682
-#> chain:4 48.048 32.326
+#> chain:1 40.605 34.244
+#> chain:2 44.485 36.481
+#> chain:3 41.615 37.638
+#> chain:4 38.434 40.516
 # summary plot
 plot(no_delay)
 ```
@@ -540,7 +534,7 @@ non_parametric <- estimate_infections(reported_cases,
 summary(non_parametric)
 #>                             measure                  estimate
 #>                              <char>                    <char>
-#> 1:           New infections per day       2730 (2689 -- 2768)
+#> 1:           New infections per day       2730 (2688 -- 2774)
 #> 2: Expected change in daily reports                Decreasing
 #> 3:       Effective reproduction no.       0.92 (0.86 -- 0.96)
 #> 4:                   Rate of growth -0.024 (-0.025 -- -0.022)
@@ -548,10 +542,10 @@ summary(non_parametric)
 # elapsed time (in seconds)
 get_elapsed_time(non_parametric$fit)
 #>         warmup sample
-#> chain:1  4.114  0.720
-#> chain:2  4.542  0.649
-#> chain:3  5.133  0.419
-#> chain:4  4.258  0.474
+#> chain:1  4.475  0.516
+#> chain:2  4.663  0.641
+#> chain:3  4.827  0.549
+#> chain:4  4.169  0.725
 # summary plot
 plot(non_parametric)
 ```
diff --git a/vignettes/estimate_infections_options.Rmd.orig b/vignettes/estimate_infections_options.Rmd.orig
index 25706703e..0cd532013 100644
--- a/vignettes/estimate_infections_options.Rmd.orig
+++ b/vignettes/estimate_infections_options.Rmd.orig
@@ -172,8 +172,8 @@ Here, instead of doing so we assume that we know about truncation with mean of 1
 
 ```{r define_truncation}
 trunc_dist <- LogNormal(
-  mean = Normal(0.5, 0.1),
-  sd = Normal(0.5, 0.1),
+  mean = 0.5,
+  sd = 0.5,
   max = 3
 )
 trunc_dist
diff --git a/vignettes/estimate_infections_workflow-plot_combined_delay-1.png b/vignettes/estimate_infections_workflow-plot_combined_delay-1.png
index 44e97d6eb..b95b0a171 100644
Binary files a/vignettes/estimate_infections_workflow-plot_combined_delay-1.png and b/vignettes/estimate_infections_workflow-plot_combined_delay-1.png differ
diff --git a/vignettes/estimate_infections_workflow-plot_uncertain_gamma-1.png b/vignettes/estimate_infections_workflow-plot_uncertain_gamma-1.png
index 0f050dd93..ae321f653 100644
Binary files a/vignettes/estimate_infections_workflow-plot_uncertain_gamma-1.png and b/vignettes/estimate_infections_workflow-plot_uncertain_gamma-1.png differ
diff --git a/vignettes/estimate_infections_workflow-results-1.png b/vignettes/estimate_infections_workflow-results-1.png
index 6ae795e8b..64afcc351 100644
Binary files a/vignettes/estimate_infections_workflow-results-1.png and b/vignettes/estimate_infections_workflow-results-1.png differ
diff --git a/vignettes/estimate_infections_workflow.Rmd b/vignettes/estimate_infections_workflow.Rmd
index d1167a8d5..50dab0fe1 100644
--- a/vignettes/estimate_infections_workflow.Rmd
+++ b/vignettes/estimate_infections_workflow.Rmd
@@ -107,31 +107,25 @@ plot(fixed_gamma)
 ![plot of chunk plot_fixed_gamma](estimate_infections_workflow-plot_fixed_gamma-1.png)
 
 If distributions are variable, the values with uncertainty are treated as [prior probability densities](https://en.wikipedia.org/wiki/Prior_probability) in the Bayesian inference framework used by _EpiNow2_, i.e. they are estimated as part of the inference.
-For example, to define a variable gamma distribution where uncertainty in the mean is given by a normal distribution with mean 3 and sd 2, and uncertainty in the standard deviation is given by a normal distribution with mean 1 and sd 0.1, with a maximum value 10, you would write
+For example, to define a variable gamma distribution where uncertainty in the shape is given by a normal distribution with mean 3 and sd 2, and uncertainty in the rate is given by a normal distribution with mean 1 and sd 0.1, with a maximum value 10, you would write
 
 
 ``` r
-uncertain_gamma <- Gamma(mean = Normal(3, 2), sd = Normal(1, 0.1), max = 10)
-#> Warning: ! Uncertain gamma distribution specified in terms of parameters that are not
-#>   the "natural" parameters of the distribution shape and rate.
-#> ℹ Converting using a crude and very approximate method that is likely to
-#>   produce biased results.
-#> ℹ If possible it is preferable to specify the distribution directly in terms of
-#>   the natural parameters.
+uncertain_gamma <- Gamma(shape = Normal(3, 2), rate = Normal(1, 0.1), max = 10)
 uncertain_gamma
 #> - gamma distribution (max: 10):
 #>   shape:
 #>     - normal distribution:
 #>       mean:
-#>         9
+#>         3
 #>       sd:
-#>         2.5
+#>         2
 #>   rate:
 #>     - normal distribution:
 #>       mean:
-#>         3
+#>         1
 #>       sd:
-#>         1.4
+#>         0.1
 ```
 
 which looks like this when plotted
@@ -142,10 +136,6 @@ plot(uncertain_gamma)
 ```
 
 ![plot of chunk plot_uncertain_gamma](estimate_infections_workflow-plot_uncertain_gamma-1.png)
-
-Note the warning about parameters.
-We used the mean and standard deviation to define this distribution with uncertain parameters, but it would be better to use the "natural" parameters of the gamma distribution, shape and rate, for example using the values estimate and reported back after calling the previous command.
-
 There are various ways the specific delay distributions mentioned below might be obtained.
 Often, they will come directly from the existing literature reviewed by the user and studies conducted elsewhere.
 Sometimes it might be possible to obtain them from existing databases, e.g. using the [epiparameter](https://github.com/epiverse-trace/epiparameter) R package.
@@ -278,7 +268,7 @@ For example, if the user believes that at the very start of the data the reprodu
 
 
 ``` r
-rt_prior <- LogNormal(mean =  2, sd = 1)
+rt_prior <- LogNormal(mean = 2, sd = 1)
 rt_opts(prior = rt_prior)
 ```
 
diff --git a/vignettes/estimate_infections_workflow.Rmd.orig b/vignettes/estimate_infections_workflow.Rmd.orig
index 2c8f48a19..6f28d6cec 100644
--- a/vignettes/estimate_infections_workflow.Rmd.orig
+++ b/vignettes/estimate_infections_workflow.Rmd.orig
@@ -90,10 +90,10 @@ plot(fixed_gamma)
 ```
 
 If distributions are variable, the values with uncertainty are treated as [prior probability densities](https://en.wikipedia.org/wiki/Prior_probability) in the Bayesian inference framework used by _EpiNow2_, i.e. they are estimated as part of the inference.
-For example, to define a variable gamma distribution where uncertainty in the mean is given by a normal distribution with mean 3 and sd 2, and uncertainty in the standard deviation is given by a normal distribution with mean 1 and sd 0.1, with a maximum value 10, you would write
+For example, to define a variable gamma distribution where uncertainty in the shape is given by a normal distribution with mean 3 and sd 2, and uncertainty in the rate is given by a normal distribution with mean 1 and sd 0.1, with a maximum value 10, you would write
 
 ```{r}
-uncertain_gamma <- Gamma(mean = Normal(3, 2), sd = Normal(1, 0.1), max = 10)
+uncertain_gamma <- Gamma(shape = Normal(3, 2), rate = Normal(1, 0.1), max = 10)
 uncertain_gamma
 ```
 
@@ -102,10 +102,6 @@ which looks like this when plotted
 ```{r plot_uncertain_gamma}
 plot(uncertain_gamma)
 ```
-
-Note the warning about parameters.
-We used the mean and standard deviation to define this distribution with uncertain parameters, but it would be better to use the "natural" parameters of the gamma distribution, shape and rate, for example using the values estimate and reported back after calling the previous command.
-
 There are various ways the specific delay distributions mentioned below might be obtained.
 Often, they will come directly from the existing literature reviewed by the user and studies conducted elsewhere.
 Sometimes it might be possible to obtain them from existing databases, e.g. using the [epiparameter](https://github.com/epiverse-trace/epiparameter) R package.
@@ -209,7 +205,7 @@ It can be changed using the `rt_opts()` function.
 For example, if the user believes that at the very start of the data the reproduction number was 2, with uncertainty in this belief represented by a standard deviation of 1, they would use
 
 ```{r results = 'hide'}
-rt_prior <- LogNormal(mean =  2, sd = 1)
+rt_prior <- LogNormal(mean = 2, sd = 1)
 rt_opts(prior = rt_prior)
 ```