diff --git a/NAMESPACE b/NAMESPACE
index c31d339a5..8a7e92267 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -38,7 +38,7 @@ export(convert_summary_stats_to_params)
 export(create_citation)
 export(create_metadata)
 export(create_method_assess)
-export(create_prob_dist)
+export(create_prob_distribution)
 export(create_region)
 export(create_summary_stats)
 export(create_uncertainty)
diff --git a/R/accessors.R b/R/accessors.R
index eb8054fcb..b90e38105 100644
--- a/R/accessors.R
+++ b/R/accessors.R
@@ -36,10 +36,10 @@ get_parameters <- function(x, ...) {
 get_parameters.epiparameter <- function(x, ...) {
   chkDots(...)
   # extract parameters depending on prob distribution class
-  if (inherits(x$prob_dist, "distcrete")) {
-    params <- unlist(x$prob_dist$parameters)
-  } else if (inherits(x$prob_dist, "distribution")) {
-    params <- unlist(distributional::parameters(x$prob_dist))
+  if (inherits(x$prob_distribution, "distcrete")) {
+    params <- unlist(x$prob_distribution$parameters)
+  } else if (inherits(x$prob_distribution, "distribution")) {
+    params <- unlist(distributional::parameters(x$prob_distribution))
 
     # if dist is truncated clean names
     if (is_truncated(x)) {
@@ -53,8 +53,8 @@ get_parameters.epiparameter <- function(x, ...) {
 
     # convert to meanlog and sdlog names
     params <- .clean_params(
-      prob_dist = family(x),
-      prob_dist_params = params
+      prob_distribution = family(x),
+      prob_distribution_params = params
     )
   } else {
     return(NA)
diff --git a/R/calc_dist_params.R b/R/calc_dist_params.R
index ad88f15e5..9d43b85cd 100644
--- a/R/calc_dist_params.R
+++ b/R/calc_dist_params.R
@@ -29,14 +29,13 @@
 #'
 #' @return A named `numeric` vector with parameters.
 #' @keywords internal
-.calc_dist_params <- function(prob_dist, # nolint cyclocomp
-                              prob_dist_params,
+.calc_dist_params <- function(prob_distribution, # nolint cyclocomp
                               summary_stats,
                               sample_size) {
-  if (is.na(prob_dist)) {
+  if (is.na(prob_distribution)) {
     message(
       "No adequate summary statistics available to calculate the parameters ",
-      "of the ", prob_dist, " distribution"
+      "of the ", prob_distribution, " distribution"
     )
     return(NA)
   }
@@ -60,19 +59,21 @@
     unlist(summary_stats$range)
   )
 
-  # extract dispersion
-  disp <- unname(prob_dist_params[names(prob_dist_params) == "dispersion"])
-  median_disp <- c(median = summary_stats$median, dispersion = disp)
-
   # extract mean and sd
   mean_sd <- c(summary_stats$mean, summary_stats$sd)
 
   is_mean_sd <- checkmate::test_numeric(
     mean_sd, any.missing = FALSE, len = 2, finite = TRUE
   )
-  is_median_disp <- checkmate::test_numeric(
-    median_disp, len = 2, finite = TRUE
+
+  is_median_disp <- checkmate::test_number(
+    summary_stats$median,
+    finite = TRUE
+  ) && checkmate::test_number(
+    summary_stats$dispersion,
+    finite = TRUE
   )
+
   is_median_range <- checkmate::test_numeric(
     median_range, len = 3, finite = TRUE
   ) && checkmate::test_count(sample_size, positive = TRUE)
@@ -93,40 +94,49 @@
     )
     summary_stats_ <- summary_stats_[idx]
     # create flat list structure to be passed to ... in conversion
-    args <- unlist(list(prob_dist, as.list(summary_stats_)), recursive = FALSE)
-    prob_dist_params <- unlist(do.call(
+    args <- unlist(
+      list(prob_distribution, as.list(summary_stats_)),
+      recursive = FALSE
+    )
+    prob_distribution_params <- unlist(do.call(
       convert_summary_stats_to_params,
       args = args
     ))
   } else if (is_median_disp) {
-    med <- summary_stats$median
-    meanlog <- log(med / sqrt(1 + disp^2))
-    sdlog <- sqrt(log(1 + disp^2))
-    prob_dist_params <- c(meanlog = meanlog, sdlog = sdlog)
+    args <- list(
+      x = prob_distribution,
+      median = summary_stats$median,
+      dispersion = summary_stats$dispersion
+    )
+    prob_distribution_params <- do.call(
+      convert_summary_stats_to_params,
+      args = args
+    )
+    prob_distribution_params <- unlist(prob_distribution_params)
   } else if (!anyNA(percentiles)) {
     # calculate the parameters from the percentiles
     # percentiles required to be [0, 1] so divide by 100
-    prob_dist_params <- extract_param(
+    prob_distribution_params <- extract_param(
       type = "percentile",
       values = percentiles,
-      distribution = prob_dist,
+      distribution = prob_distribution,
       percentiles = as.numeric(names(percentiles)) / 100
     )
   } else if (is_median_range) {
-    prob_dist_params <- extract_param(
+    prob_distribution_params <- extract_param(
       type = "range",
       values = median_range,
-      distribution = prob_dist,
+      distribution = prob_distribution,
       samples = sample_size
     )
   } else {
     message(
       "No adequate summary statistics available to calculate the parameters ",
-      "of the ", prob_dist, " distribution"
+      "of the ", prob_distribution, " distribution"
     )
     return(NA)
   }
 
   # return params
-  prob_dist_params
+  prob_distribution_params
 }
diff --git a/R/checkers.R b/R/checkers.R
index c527078e7..2de76462b 100644
--- a/R/checkers.R
+++ b/R/checkers.R
@@ -16,8 +16,10 @@
 #' ep <- epiparameter(
 #'   disease = "ebola",
 #'   epi_dist = "incubation",
-#'   prob_distribution = "gamma",
-#'   prob_distribution_params = c(shape = 1, scale = 1)
+#'   prob_distribution = create_prob_distribution(
+#'     prob_distribution = "gamma",
+#'     prob_distribution_params = c(shape = 1, scale = 1)
+#'   )
 #' )
 #' is_parameterised(ep)
 #'
@@ -43,7 +45,7 @@ is_parameterized <- is_parameterised
 is_parameterised.epiparameter <- function(x, ...) {
   chkDots(...)
   # probability distribution object
-  return(is.object(x$prob_dist))
+  return(is.object(x$prob_distribution))
 }
 
 #' @export
diff --git a/R/coercion.R b/R/coercion.R
index 75821967e..098f61a8a 100644
--- a/R/coercion.R
+++ b/R/coercion.R
@@ -62,7 +62,7 @@ as.data.frame.epiparameter <- function(x, ...) {
     disease = x$disease,
     pathogen = x$pathogen,
     epi_distribution = x$epi_dist,
-    prob_distribution = I(list(prob_dist = x$prob_dist)),
+    prob_distribution = I(list(prob_dist = x$prob_distribution)),
     uncertainty = I(list(uncertainty = x$uncertainty)),
     summary_stats = I(list(summary_stats = x$summary_stats)),
     citation = I(x$citation),
@@ -220,15 +220,17 @@ is_epiparameter_df <- function(x) {
     disease = x$disease,
     pathogen = x$pathogen,
     epi_dist = x$epi_distribution,
-    prob_distribution = family(x$prob_distribution[[1]]),
-    prob_distribution_params = prob_distribution_params,
+    prob_distribution = create_prob_distribution(
+      prob_distribution = family(x$prob_distribution[[1]]),
+        prob_distribution_params = prob_distribution_params,
+      discretise = inherits(x$prob_distribution[[1]], "distcrete"),
+      truncation = truncation
+    ),
     uncertainty = uncertainty,
     summary_stats = x$summary_stats$summary_stats,
     citation = x$citation,
     metadata = x$metadata,
     method_assess = x$method_assess,
-    discretise = inherits(x$prob_distribution[[1]], "distcrete"),
-    truncation = truncation,
     notes = x$notes
   )
 }
@@ -268,7 +270,7 @@ is_epiparameter_df <- function(x) {
   sd_ <- NULL
   if (rlang::is_na(prob_dist)) {
     prob_dist_params <- NA_real_
-    uncertainty <- create_uncertainty()
+    uncertainty <- list(uncertainty = create_uncertainty())
   } else {
     prob_dist <- switch(
       prob_dist,
@@ -309,6 +311,14 @@ is_epiparameter_df <- function(x) {
     )
     prob_dist_params_names <- .clean_string(prob_dist_params_names)
     names(prob_dist_params) <- prob_dist_params_names
+
+    # overwrite prob_dist with user specified if given to make use of conversion
+    if (!is.null(prob_dist_in)) {
+      prob_dist <- prob_dist_in
+      # remove user specified to not trigger overwriting prob_dist below
+      prob_dist_in <- NULL
+    }
+
     if (all(c("mean", "sd") %in% names(prob_dist_params))) {
       sd_ <- prob_dist_params[["sd"]]
       prob_dist_params <- do.call(
@@ -422,8 +432,10 @@ is_epiparameter_df <- function(x) {
     disease = disease,
     pathogen = pathogen,
     epi_dist = epi_dist,
-    prob_distribution = prob_dist,
-    prob_distribution_params = prob_dist_params,
+    prob_distribution = create_prob_distribution(
+      prob_distribution = prob_dist,
+      prob_distribution_params = prob_dist_params
+    ),
     uncertainty = uncertainty,
     summary_stats = summary_stats,
     citation = citation,
diff --git a/R/convert_params.R b/R/convert_params.R
index a5fe41b20..f4a2beea0 100644
--- a/R/convert_params.R
+++ b/R/convert_params.R
@@ -230,7 +230,7 @@ convert_params_to_summary_stats.epiparameter <- function(x, ...) {
   if (!is_parameterised(x)) {
     if (length(dots) > 0) {
       # unparameterised with parameters supplied through `...`
-      return(convert_params_to_summary_stats(x$prob_dist, ...))
+      return(convert_params_to_summary_stats(x$prob_distribution, ...))
     }
     # unparameterised with no parameters supplied through `...`
     stop(
diff --git a/R/create_prob_dist.R b/R/create_prob_distribution.R
similarity index 50%
rename from R/create_prob_dist.R
rename to R/create_prob_distribution.R
index 282c71e0c..563e80652 100644
--- a/R/create_prob_dist.R
+++ b/R/create_prob_distribution.R
@@ -22,7 +22,17 @@
 #' for the discretisation is 1. `w` can be `[0,1]`. For more information please
 #' see [distcrete::distcrete()].
 #'
-#' @inheritParams new_epiparameter
+#' @param prob_distribution A `character` string specifying the probability
+#' distribution. This should match the \R naming convention of probability
+#' distributions (e.g. lognormal is `lnorm`, negative binomial is `nbinom`, and
+#' geometric is `geom`).
+#' @param prob_distribution_params A named vector of probability distribution
+#' parameters.
+#' @param discretise A boolean `logical` whether the distribution is
+#' discretised. Default is `FALSE` which assumes a continuous probability
+#' distribution.
+#' @param truncation A `numeric` specifying the truncation point if the inferred
+#' distribution was truncated, `NA` if not or unknown.
 #' @param ... [dots] Extra arguments to be passed to
 #' \pkg{distributional} or \pkg{distcrete} functions that construct the S3
 #' distribution objects. To see which arguments can be adjusted for discretised
@@ -31,74 +41,99 @@
 #' constructor function, e.g. for the Gamma distribution see
 #' [distributional::dist_gamma()].
 #'
-#' @return An S3 class containing the probability distribution.
+#' @return An S3 class containing the probability distribution or a `character`
+#' string if the parameters of the probability distribution are unknown.
 #' @export
 #'
 #' @examples
 #' # example with continuous distribution without truncation
-#' create_prob_dist(
-#'   prob_dist = "gamma",
-#'   prob_dist_params = c(shape = 1, scale = 1),
+#' create_prob_distribution(
+#'   prob_distribution = "gamma",
+#'   prob_distribution_params = c(shape = 1, scale = 1),
 #'   discretise = FALSE,
 #'   truncation = NA
 #' )
 #'
 #' # example with continuous distribution with truncation
-#' create_prob_dist(
-#'   prob_dist = "gamma",
-#'   prob_dist_params = c(shape = 1, scale = 1),
+#' create_prob_distribution(
+#'   prob_distribution = "gamma",
+#'   prob_distribution_params = c(shape = 1, scale = 1),
 #'   discretise = FALSE,
 #'   truncation = 10
 #' )
 #'
 #' # example with discrete distribution
-#' create_prob_dist(
-#'   prob_dist = "gamma",
-#'   prob_dist_params = c(shape = 1, scale = 1),
+#' create_prob_distribution(
+#'   prob_distribution = "gamma",
+#'   prob_distribution_params = c(shape = 1, scale = 1),
 #'   discretise = TRUE,
 #'   truncation = NA
 #' )
 #'
 #' # example passing extra arguments to distcrete
-#' create_prob_dist(
-#'   prob_dist = "gamma",
-#'   prob_dist_params = c(shape = 1, scale = 1),
+#' create_prob_distribution(
+#'   prob_distribution = "gamma",
+#'   prob_distribution_params = c(shape = 1, scale = 1),
 #'   discretise = TRUE,
 #'   truncation = NA,
 #'   w = 0.5
 #' )
-create_prob_dist <- function(prob_dist,
-                             prob_dist_params,
-                             discretise = FALSE,
-                             truncation = NA,
-                             ...) {
+create_prob_distribution <- function(prob_distribution,
+                                     prob_distribution_params,
+                                     discretise = FALSE,
+                                     truncation = NA,
+                                     ...) {
   checkmate::assert_character(
-    prob_dist,
+    prob_distribution,
     min.chars = 1,
     min.len = 1,
     max.len = 2
   )
-  checkmate::assert_numeric(prob_dist_params, names = "unique")
+
+  # when only the type of probability distribution is known return string
+  if (missing(prob_distribution_params) || anyNA(prob_distribution_params)) {
+    return(prob_distribution)
+  }
+
+  # NA parameters will be caught by if above so only need to check numeric
+  stopifnot(
+    "`prob_distribution_params` must be a named vector of numerics or NA" =
+      checkmate::test_numeric(prob_distribution_params, names = "unique")
+  )
   checkmate::assert_logical(discretise, len = 1)
   checkmate::assert_number(truncation, na.ok = TRUE)
 
+  # set prob_distribution to lowercase for downstream case sensitive matching
+  prob_distribution <- .clean_string(prob_distribution)
+
+  stopifnot(
+    "Incorrect parameters provided for probability distribution." =
+      is_epiparameter_params(prob_distribution, prob_distribution_params)
+  )
+
+  # standardise common distribution parameters
+  prob_distribution_params <- .clean_params(
+    prob_distribution = prob_distribution,
+    prob_distribution_params = prob_distribution_params
+  )
+
   dots <- list(...)
   if (discretise) {
-    prob_dist <- match.arg(
-      prob_dist,
+    prob_distribution <- match.arg(
+      prob_distribution,
       choices = c("gamma", "lnorm", "weibull", "norm")
     )
     # create default list of args to construct <distcrete>
     distcrete_args <- c(
-      name = prob_dist,
+      name = prob_distribution,
       interval = 1,
-      as.list(prob_dist_params),
+      as.list(prob_distribution_params),
       w = 1
     )
     # modify <distcrete> args if provided in dots
     distcrete_args <- utils::modifyList(distcrete_args, dots)
     # create discretised probability distribution object
-    prob_dist <- do.call(
+    prob_distribution <- do.call(
       distcrete::distcrete,
       distcrete_args
     )
@@ -106,39 +141,39 @@ create_prob_dist <- function(prob_dist,
     # currently dots not used to construct <distribution>
     chkDots(...)
     # create non-discretised probability distribution object
-    prob_dist <- switch(prob_dist,
+    prob_distribution <- switch(prob_distribution,
       gamma = distributional::dist_gamma(
-        shape = prob_dist_params[["shape"]],
-        rate = 1 / prob_dist_params[["scale"]]
+        shape = prob_distribution_params[["shape"]],
+        rate = 1 / prob_distribution_params[["scale"]]
       ),
       lnorm = distributional::dist_lognormal(
-        mu = prob_dist_params[["meanlog"]],
-        sigma = prob_dist_params[["sdlog"]]
+        mu = prob_distribution_params[["meanlog"]],
+        sigma = prob_distribution_params[["sdlog"]]
       ),
       weibull = distributional::dist_weibull(
-        shape = prob_dist_params[["shape"]],
-        scale = prob_dist_params[["scale"]]
+        shape = prob_distribution_params[["shape"]],
+        scale = prob_distribution_params[["scale"]]
       ),
       nbinom = distributional::dist_negative_binomial(
-        size = prob_dist_params[["dispersion"]],
+        size = prob_distribution_params[["dispersion"]],
         prob = convert_summary_stats_to_params(
           "nbinom",
-          mean = prob_dist_params[["mean"]],
-          dispersion = prob_dist_params[["dispersion"]]
+          mean = prob_distribution_params[["mean"]],
+          dispersion = prob_distribution_params[["dispersion"]]
         )$prob
       ),
       geom = distributional::dist_geometric(
-        prob = unname(prob_dist_params)
+        prob = unname(prob_distribution_params)
       ),
       pois = distributional::dist_poisson(
-        lambda = unname(prob_dist_params)
+        lambda = unname(prob_distribution_params)
       ),
       norm = distributional::dist_normal(
-        mu = prob_dist_params[["mean"]],
-        sigma = prob_dist_params[["sd"]]
+        mu = prob_distribution_params[["mean"]],
+        sigma = prob_distribution_params[["sd"]]
       ),
       exp = distributional::dist_exponential(
-        rate = prob_dist_params[["rate"]]
+        rate = prob_distribution_params[["rate"]]
       ),
       stop("Did not recognise distribution name", call. = FALSE)
     )
@@ -152,13 +187,13 @@ create_prob_dist <- function(prob_dist,
         call. = FALSE
       )
     } else {
-      prob_dist <- distributional::dist_truncated(
-        prob_dist,
+      prob_distribution <- distributional::dist_truncated(
+        prob_distribution,
         upper = truncation
       )
     }
   }
 
-  # return prob_dist object
-  prob_dist
+  # return prob_distribution object
+  prob_distribution
 }
diff --git a/R/epiparameter-utils.R b/R/epiparameter-utils.R
index eb95b5a5a..7896fb5a0 100644
--- a/R/epiparameter-utils.R
+++ b/R/epiparameter-utils.R
@@ -484,32 +484,33 @@ create_method_assess <- function(censored = NA,
 #' This check for valid parameters is independent of whether the distribution
 #' is truncated or discretised.
 #'
-#' @inheritParams new_epiparameter
+#' @inheritParams create_prob_distribution
 #'
 #' @return A boolean `logical`.
 #' @export
 #'
 #' @examples
 #' is_epiparameter_params(
-#'   prob_dist = "gamma",
-#'   prob_dist_params = c(shape = 2, scale = 1)
+#'   prob_distribution = "gamma",
+#'   prob_distribution_params = c(shape = 2, scale = 1)
 #' )
-is_epiparameter_params <- function(prob_dist, prob_dist_params) {
-  if (is.na(prob_dist) || anyNA(prob_dist_params)) {
+is_epiparameter_params <- function(prob_distribution,
+                                   prob_distribution_params) {
+  if (is.na(prob_distribution) || anyNA(prob_distribution_params)) {
     return(FALSE)
   }
 
   # check input
   checkmate::assert_numeric(
-    prob_dist_params,
+    prob_distribution_params,
     min.len = 1,
     names = "unique"
   )
 
   # remove truncation parameters if truncated
-  if ("upper" %in% names(prob_dist_params)) {
-    prob_dist_params <- prob_dist_params[
-      names(prob_dist_params) != c("lower", "upper")
+  if ("upper" %in% names(prob_distribution_params)) {
+    prob_distribution_params <- prob_distribution_params[
+      names(prob_distribution_params) != c("lower", "upper")
     ]
   }
 
@@ -524,13 +525,13 @@ is_epiparameter_params <- function(prob_dist, prob_dist_params) {
     norm = list(c("mean", "sd"), c("mu", "sigma")),
     exp = list("rate", "lambda", "mean")
   )
-  possible_params <- possible_params[[prob_dist]]
+  possible_params <- possible_params[[prob_distribution]]
 
   # check whether any combinations are valid
   matches <- vapply(
     possible_params,
     function(x, y) all(names(y) %in% x) && identical(length(y), length(x)),
-    y = prob_dist_params,
+    y = prob_distribution_params,
     FUN.VALUE = logical(1)
   )
 
@@ -552,36 +553,36 @@ is_epiparameter_params <- function(prob_dist, prob_dist_params) {
 #' This means that the distribution specific parameter cleaning functions do
 #' not need to check and error for incorrect parameterisation.
 #'
-#' @inheritParams new_epiparameter
+#' @inheritParams create_prob_distribution
 #'
 #' @name .clean_params
 #'
 #' @return Named `numeric` vector of parameters.
 #' @keywords internal
-.clean_params <- function(prob_dist, prob_dist_params) {
+.clean_params <- function(prob_distribution, prob_distribution_params) {
   valid_params <- is_epiparameter_params(
-    prob_dist = prob_dist,
-    prob_dist_params = prob_dist_params
+    prob_distribution = prob_distribution,
+    prob_distribution_params = prob_distribution_params
   )
   if (!valid_params) {
     stop(
-      "Invalid parameterisation for ", prob_dist, " distribution",
+      "Invalid parameterisation for ", prob_distribution, " distribution",
       call. = FALSE
     )
   }
-  is_trunc <- "upper" %in% names(prob_dist_params)
+  is_trunc <- "upper" %in% names(prob_distribution_params)
   # remove truncation parameters if truncated
   if (is_trunc) {
-    trunc_params <- prob_dist_params[
-      names(prob_dist_params) == c("lower", "upper")
+    trunc_params <- prob_distribution_params[
+      names(prob_distribution_params) == c("lower", "upper")
     ]
-    prob_dist_params <- prob_dist_params[
-      names(prob_dist_params) != c("lower", "upper")
+    prob_distribution_params <- prob_distribution_params[
+      names(prob_distribution_params) != c("lower", "upper")
     ]
   }
   # weibull only has one parameterisation so does not need cleaning
   clean_func <- switch(
-    prob_dist,
+    prob_distribution,
     gamma = .clean_params_gamma,
     lnorm = .clean_params_lnorm,
     weibull = function(x) x,
@@ -592,7 +593,7 @@ is_epiparameter_params <- function(prob_dist, prob_dist_params) {
     exp = .clean_params_exp,
     stop("Probability distribution not recognised", call. = FALSE)
   )
-  clean_params <- do.call(clean_func, list(prob_dist_params))
+  clean_params <- do.call(clean_func, list(prob_distribution_params))
   # reappend truncation parameter if truncated
   if (is_trunc) {
     clean_params <- append(clean_params, trunc_params)
@@ -736,3 +737,52 @@ is_epiparameter_params <- function(prob_dist, prob_dist_params) {
   checkmate::assert_character(x)
   gsub(pattern = "_|-", replacement = " ", x = trimws(tolower(x)))
 }
+
+#' Standardise distribution parameter uncertainty
+#'
+#' @param x An `<epiparameter>` object.
+#' @inheritParams epiparameter
+#' @param uncertainty_missing A boolean `logical` as to whether `uncertainty`
+#' is missing (see [missing()]) from the parent function.
+#'
+#' @return An uncertainty list for an `<epiparameter>` object.
+#' @keywords internal
+.clean_uncertainty <- function(x, prob_distribution, uncertainty_missing) {
+  param_names <- names(get_parameters(x))
+  param_names <- param_names %||% NA_character_
+  if (uncertainty_missing ||
+      !identical(prob_distribution, x$prob_distribution)) {
+    # create uncertainty for each parameter if not provided or auto calculated
+    x$uncertainty <- lapply(
+      param_names,
+      function(xx) list(uncertainty = create_uncertainty())
+    )
+    if (!anyNA(param_names)) names(x$uncertainty) <- param_names
+  } else {
+    ci_limits_ <- lapply(x$uncertainty, `[[`, 1)
+    # if uncertainty is unspecified then it only needs renaming
+    if (anyNA(ci_limits_, recursive = TRUE)) {
+      if (!anyNA(param_names)) names(x$uncertainty) <- param_names
+    } else {
+      # standardise parameter uncertainty in to match parameters
+      dist <- family(x)
+      params_ <- vector(mode = "list", length = length(x$uncertainty))
+      names(params_) <- param_names
+      for (i in seq_along(ci_limits_)) {
+        params <- vapply(ci_limits_, `[[`, FUN.VALUE = numeric(1), i)
+        temp <-  as.list(.clean_params(
+          prob_distribution = dist,
+          prob_distribution_params = params
+        ))
+        params_ <- Map(c, params_, temp)
+      }
+      ci_limits_ <- lapply(params_, function(x) sort(unname(x)))
+      for (i in seq_along(x$uncertainty)) {
+        x$uncertainty[[i]]$ci_limits <- ci_limits_[[i]]
+      }
+      names(x$uncertainty) <- param_names
+    }
+  }
+  # return <epiparameter> uncertainty
+  x$uncertainty
+}
diff --git a/R/epiparameter.R b/R/epiparameter.R
index e89d7f927..9dd2031cc 100644
--- a/R/epiparameter.R
+++ b/R/epiparameter.R
@@ -9,12 +9,6 @@
 #' `{distributional}` when `discretise = FALSE`, or a `distcrete` object from
 #' `{distcrete}` when `discretise = TRUE`.
 #'
-#' @param prob_dist A character string specifying the probability
-#' distribution. This should match the R naming convention of probability
-#' distributions (e.g. lognormal is `lnorm`, negative binomial is `nbinom`, and
-#' geometric is `geom`).
-#' @param prob_dist_params A named vector of probability distribution
-#' parameters.
 #' @inheritParams epiparameter
 #'
 #' @inherit epiparameter return
@@ -22,69 +16,37 @@
 new_epiparameter <- function(disease = character(),
                              pathogen = character(),
                              epi_dist = character(),
-                             prob_dist = list(),
-                             prob_dist_params = numeric(),
+                             prob_distribution = list(),
                              uncertainty = list(),
                              summary_stats = list(),
-                             auto_calc_params = logical(),
                              citation = character(),
                              metadata = list(),
                              method_assess = list(),
-                             discretise = logical(),
-                             truncation = numeric(),
                              notes = character(),
+                             auto_calc_params = logical(),
                              ...) {
-  # check ci has been given for each param and param and uncertainty names match
-  stopifnot(
-    "uncertainty must be provided for each parameter" =
-      anyNA(uncertainty) ||
-      length(prob_dist_params) == length(uncertainty),
-    "parameters and uncertainty must be named and match" =
-      anyNA(uncertainty) ||
-      identical(names(prob_dist_params), names(uncertainty))
-  )
-
-  # set string to lowercase for downstream case sensitive matching of prob_dist
-  prob_dist <- ifelse(
-    test = is.character(prob_dist),
-    yes = .clean_string(prob_dist),
-    no = prob_dist
-  )
-
-  # TODO: formalise if statement below or remove
-  # include mean in prob_dist_params
-  if (!is.null(summary_stats$mean) && !is.na(summary_stats$mean) &&
-      prob_dist %in% c("nbinom", "geom", "pois", "norm")) {
-    prob_dist_params["mean"] <- summary_stats$mean
-  }
-
-  if (is_epiparameter_params(prob_dist, prob_dist_params)) {
-    # standardise common distribution parameters
-    prob_dist_params <- .clean_params(
-      prob_dist = prob_dist,
-      prob_dist_params = prob_dist_params
-    )
-  } else if (auto_calc_params) {
+  if (auto_calc_params && is.character(prob_distribution)) {
     # calculate parameters if not provided
-    prob_dist_params <- .calc_dist_params(
-      prob_dist = prob_dist,
-      prob_dist_params = prob_dist_params,
+    prob_distribution_params <- .calc_dist_params(
+      prob_distribution = prob_distribution,
       summary_stats = summary_stats,
       sample_size = metadata$sample_size
     )
-  }
-
-  if (anyNA(prob_dist_params)) {
-    message("Unparameterised <epiparameter> object")
-  } else {
-    # create a S3 object holding the probability distribution
-    prob_dist <- create_prob_dist(
-      prob_dist = prob_dist,
-      prob_dist_params = prob_dist_params,
-      discretise = discretise,
-      truncation = truncation,
-      ...
-    )
+    if (!anyNA(prob_distribution_params)) {
+      prob_distribution <- create_prob_distribution(
+        prob_distribution = prob_distribution,
+        prob_distribution_params = prob_distribution_params,
+        ...
+      )
+      message(
+        "Parameterising the probability distribution with the summary ",
+        "statistics.\n Probability distribution is assumed not to be ",
+        "discretised or truncated."
+      )
+    }
+    if (!inherits(prob_distribution, c("distribution", "distcrete"))) {
+      message("Unparameterised <epiparameter> object")
+    }
   }
 
   if (epi_dist == "offspring_distribution") {
@@ -104,7 +66,7 @@ new_epiparameter <- function(disease = character(),
       disease = disease,
       pathogen = pathogen,
       epi_dist = epi_dist,
-      prob_dist = prob_dist,
+      prob_distribution = prob_distribution,
       uncertainty = uncertainty,
       summary_stats = summary_stats,
       citation = citation,
@@ -152,12 +114,11 @@ new_epiparameter <- function(disease = character(),
 #' disease, or `NA` if not known.
 #' @param epi_dist A `character` string with the name of the
 #' epidemiological distribution type.
-#' @param prob_distribution A `character` string specifying the probability
-#' distribution. This should match the \R naming convention of probability
-#' distributions (e.g. lognormal is `lnorm`, negative binomial is `nbinom`, and
-#' geometric is `geom`).
-#' @param prob_distribution_params A named vector of probability distribution
-#' parameters.
+#' @param prob_distribution An S3 class containing the probability
+#' distribution or a character string if the parameters of the probability
+#' distribution are unknown but the name of the distribution is known, or `NA`
+#' if the distribution name and parameters are unknown. Use
+#' [create_prob_distribution()] to create `prob_distribution`.
 #' @param uncertainty A list of named vectors with the uncertainty around
 #' the probability distribution parameters. If uncertainty around the parameter
 #' estimates is unknown use [create_uncertainty()] (which is the
@@ -189,11 +150,6 @@ new_epiparameter <- function(disease = character(),
 #' @param method_assess A list of methodological aspects used when fitting
 #' the distribution, use [create_method_assess()] to create method
 #' assessment.
-#' @param discretise A boolean `logical` whether the distribution is
-#' discretised.
-#' Default is FALSE which assumes a continuous probability distribution
-#' @param truncation A `numeric` specifying the truncation point if the inferred
-#' distribution was truncated, `NA` if not or unknown.
 #' @param notes A `character` string with any additional information about the
 #' data, inference method or disease.
 #' @param ... [dots] Extra arguments to be passed to internal functions.
@@ -211,17 +167,21 @@ new_epiparameter <- function(disease = character(),
 #' ebola_incubation <- epiparameter(
 #'   disease = "ebola",
 #'   epi_dist = "incubation_period",
-#'   prob_distribution = "gamma",
-#'   prob_distribution_params = c(shape = 1, scale = 1)
+#'   prob_distribution = create_prob_distribution(
+#'     prob_distribution = "gamma",
+#'     prob_distribution_params = c(shape = 1, scale = 1)
+#'   )
 #' )
 #'
 #' # minimal input required for discrete `epiparameter`
 #' ebola_incubation <- epiparameter(
 #'   disease = "ebola",
 #'   epi_dist = "incubation_period",
-#'   prob_distribution = "gamma",
-#'   prob_distribution_params = c(shape = 1, scale = 1),
-#'   discretise = TRUE
+#'   prob_distribution = create_prob_distribution(
+#'     prob_distribution = "gamma",
+#'     prob_distribution_params = c(shape = 1, scale = 1),
+#'     discretise = TRUE
+#'   )
 #' )
 #'
 #' # example with more fields filled in
@@ -229,9 +189,16 @@ new_epiparameter <- function(disease = character(),
 #'   disease = "ebola",
 #'   pathogen = "ebola_virus",
 #'   epi_dist = "incubation",
-#'   prob_distribution = "gamma",
-#'   prob_distribution_params = c(shape = 1, scale = 1),
-#'   uncertainty = create_uncertainty(),
+#'   prob_distribution = create_prob_distribution(
+#'     prob_distribution = "gamma",
+#'     prob_distribution_params = c(shape = 1, scale = 1),
+#'     discretise = FALSE,
+#'     truncation = NA
+#'   ),
+#'   uncertainty = list(
+#'     shape = create_uncertainty(),
+#'     scale = create_uncertainty()
+#'   ),
 #'   summary_stats = create_summary_stats(mean = 2, sd = 1),
 #'   citation = create_citation(
 #'     author = person(given = "John", family = "Smith"),
@@ -249,34 +216,30 @@ new_epiparameter <- function(disease = character(),
 #'   method_assess = create_method_assess(
 #'     censored = TRUE
 #'   ),
-#'   discretise = FALSE,
-#'   truncation = NA,
 #'   notes = "No notes"
 #' )
 epiparameter <- function(disease,
                          pathogen = NA_character_,
                          epi_dist,
-                         prob_distribution = NA_character_,
-                         prob_distribution_params = NA_real_,
+                         prob_distribution = create_prob_distribution(
+                           prob_distribution = NA_character_
+                         ),
                          uncertainty = create_uncertainty(),
                          summary_stats = create_summary_stats(),
-                         auto_calc_params = TRUE,
                          citation = create_citation(),
                          metadata = create_metadata(),
                          method_assess = create_method_assess(),
-                         discretise = FALSE,
-                         truncation = NA_real_,
                          notes = NULL,
+                         auto_calc_params = TRUE,
                          ...) {
   # check input
   checkmate::assert_string(disease)
   checkmate::assert_character(pathogen)
   checkmate::assert_string(epi_dist)
-  checkmate::assert_character(
-    prob_distribution,
-    min.chars = 1,
-    min.len = 1,
-    max.len = 2
+  stopifnot(
+    "Probability distribution must be a distribution object or a character" =
+      !inherits(prob_distribution, c("distribution", "distcrete")) ||
+      !is.character(prob_distribution)
   )
   checkmate::assert_list(uncertainty, names = "unique")
   checkmate::assert_list(
@@ -288,40 +251,41 @@ epiparameter <- function(disease,
   checkmate::assert_class(citation, classes = "bibentry")
   checkmate::assert_list(metadata)
   checkmate::assert_list(method_assess)
-  checkmate::assert_number(truncation, na.ok = TRUE)
-  checkmate::assert_logical(discretise, len = 1)
   checkmate::assert_character(notes, null.ok = TRUE)
 
-  # check whether ci has been provided for each parameter
-  # check whether probability params are named or na
-  stopifnot(
-    "uncertainty must be provided for each parameter" =
-      anyNA(uncertainty, recursive = TRUE) ||
-        length(prob_distribution_params) == length(uncertainty),
-    "probability distribution params must be a named vector or NA" =
-      anyNA(prob_distribution_params, recursive = TRUE) ||
-        !is.null(names(prob_distribution_params))
-  )
-
   # call epiparameter constructor
   epiparameter <- new_epiparameter(
     disease = disease,
     pathogen = pathogen,
     epi_dist = epi_dist,
-    prob_dist = prob_distribution,
-    prob_dist_params = prob_distribution_params,
+    prob_distribution = prob_distribution,
     uncertainty = uncertainty,
     summary_stats = summary_stats,
     auto_calc_params = auto_calc_params,
     citation = citation,
     metadata = metadata,
     method_assess = method_assess,
-    discretise = discretise,
-    truncation = truncation,
     notes = notes,
     ...
   )
 
+  # uncertainty is checked after new_epiparameter to use <epiparameter> methods
+  epiparameter$uncertainty <- .clean_uncertainty(
+    epiparameter,
+    prob_distribution = prob_distribution,
+    uncertainty_missing = missing(uncertainty)
+  )
+
+  param_names <- names(get_parameters(epiparameter)) %||% NA_character_
+
+  stopifnot(
+    "uncertainty must be provided for each parameter" =
+      length(param_names) == length(epiparameter$uncertainty),
+    "parameters and uncertainty must be named and match" =
+      identical(param_names, names(epiparameter$uncertainty)) ||
+      is.na(param_names)
+  )
+
   # call epiparameter validator
   assert_epiparameter(epiparameter)
 
@@ -343,7 +307,7 @@ assert_epiparameter <- function(x) {
   }
 
   list_names <- c(
-    "disease", "pathogen", "epi_dist", "prob_dist", "uncertainty",
+    "disease", "pathogen", "epi_dist", "prob_distribution", "uncertainty",
     "summary_stats", "citation", "metadata", "method_assess", "notes"
   )
   missing_list_names <- list_names[!list_names %in% attributes(x)$names]
@@ -361,8 +325,8 @@ assert_epiparameter <- function(x) {
       checkmate::test_string(x$epi_dist),
     "epiparameter must contain a <distribution> or <distcrete> or NA" =
     checkmate::test_multi_class(
-      x$prob_dist, classes = c("distribution", "distcrete")
-    ) || checkmate::test_string(x$prob_dist, na.ok = TRUE),
+      x$prob_distribution, classes = c("distribution", "distcrete")
+    ) || checkmate::test_string(x$prob_distribution, na.ok = TRUE),
     "epidisit must contain uncertainty, summary stats and metadata" =
       all(
         is.list(x$uncertainty), is.list(x$summary_stats), is.list(x$metadata)
@@ -387,7 +351,7 @@ test_epiparameter <- function(x) { # nolint cyclocomp_linter
   if (!is_epiparameter(x)) return(FALSE)
 
   list_names <- c(
-    "disease", "pathogen", "epi_dist", "prob_dist", "uncertainty",
+    "disease", "pathogen", "epi_dist", "prob_distribution", "uncertainty",
     "summary_stats", "citation", "metadata", "method_assess", "notes"
   )
   missing_list_names <- list_names[!list_names %in% attributes(x)$names]
@@ -396,8 +360,8 @@ test_epiparameter <- function(x) { # nolint cyclocomp_linter
   valid_elements <- checkmate::test_string(x$disease) &&
     checkmate::test_string(x$epi_dist) &&
     (checkmate::test_multi_class(
-      x$prob_dist, classes = c("distribution", "distcrete")
-    ) || checkmate::test_string(x$prob_dist, na.ok = TRUE)) &&
+      x$prob_distribution, classes = c("distribution", "distcrete")
+    ) || checkmate::test_string(x$prob_distribution, na.ok = TRUE)) &&
     all(
       is.list(x$uncertainty), is.list(x$summary_stats), is.list(x$metadata)
     ) &&
@@ -420,8 +384,10 @@ test_epiparameter <- function(x) { # nolint cyclocomp_linter
 #' epiparameter <- epiparameter(
 #'   disease = "ebola",
 #'   epi_dist = "incubation_period",
-#'   prob_distribution = "gamma",
-#'   prob_distribution_params = c(shape = 1, scale = 1)
+#'   prob_distribution = create_prob_distribution(
+#'     prob_distribution = "gamma",
+#'     prob_distribution_params = c(shape = 1, scale = 1)
+#'   )
 #' )
 #' epiparameter
 print.epiparameter <- function(x, ...) {
@@ -440,8 +406,10 @@ print.epiparameter <- function(x, ...) {
 #' epiparameter <- epiparameter(
 #'   disease = "ebola",
 #'   epi_dist = "incubation_period",
-#'   prob_distribution = "gamma",
-#'   prob_distribution_params = c(shape = 1, scale = 1)
+#'   prob_distribution = create_prob_distribution(
+#'     prob_distribution = "gamma",
+#'     prob_distribution_params = c(shape = 1, scale = 1)
+#'   )
 #' )
 #' format(epiparameter)
 format.epiparameter <- function(x, ...) {
@@ -454,9 +422,9 @@ format.epiparameter <- function(x, ...) {
       )
     )
 
-  if (is.object(x$prob_dist) || is.character(x$prob_dist)) {
+  if (is.object(x$prob_distribution) || is.character(x$prob_distribution)) {
     dist_string <- ifelse(
-      test = inherits(x$prob_dist, "distcrete"),
+      test = inherits(x$prob_distribution, "distcrete"),
       yes = tr_("Distribution: discrete %s"),
       no = tr_("Distribution: %s")
     )
@@ -465,7 +433,7 @@ format.epiparameter <- function(x, ...) {
     writeLines(tr_("Parameters: <no parameters>"))
   }
 
-  if (is.object(x$prob_dist)) {
+  if (is.object(x$prob_distribution)) {
     params <- get_parameters(x)
 
     # decide on parameter format from magnitude of number
@@ -502,8 +470,10 @@ format.epiparameter <- function(x, ...) {
 #' ep <- epiparameter(
 #'   disease = "ebola",
 #'   epi_dist = "serial_interval",
-#'   prob_distribution = "gamma",
-#'   prob_distribution_params = c(shape = 1, scale = 1)
+#'   prob_distribution = create_prob_distribution(
+#'     prob_distribution = "gamma",
+#'     prob_distribution_params = c(shape = 1, scale = 1)
+#'   )
 #' )
 #'
 #' is_epiparameter(ep)
@@ -546,8 +516,10 @@ is_epiparameter <- function(x) {
 #' ep <- epiparameter(
 #'   disease = "ebola",
 #'   epi_dist = "incubation_period",
-#'   prob_distribution = "gamma",
-#'   prob_distribution_params = c(shape = 1, scale = 1)
+#'   prob_distribution = create_prob_distribution(
+#'     prob_distribution = "gamma",
+#'     prob_distribution_params = c(shape = 1, scale = 1)
+#'   )
 #' )
 #'
 #' # example of each distribution method for an `epiparameter` object
@@ -564,11 +536,11 @@ density.epiparameter <- function(x, at, ...) {
   if (isFALSE(is_parameterised(x))) {
     stop("<epiparameter> is unparameterised", call. = FALSE)
   }
-  unlist <- length(x$prob_dist) == 1
-  if (inherits(x$prob_dist, "distcrete")) {
-    out <- x$prob_dist$d(at)
+  unlist <- length(x$prob_distribution) == 1
+  if (inherits(x$prob_distribution, "distcrete")) {
+    out <- x$prob_distribution$d(at)
   } else {
-    out <- stats::density(x$prob_dist, at = at)
+    out <- stats::density(x$prob_distribution, at = at)
   }
   out <- if (unlist) unlist(out, recursive = FALSE) else out
   out
@@ -585,12 +557,12 @@ cdf.epiparameter <- function(x, q, ..., log = FALSE) {
   if (isFALSE(is_parameterised(x))) {
     stop("<epiparameter> is unparameterised", call. = FALSE)
   }
-  unlist <- length(x$prob_dist) == 1
-  if (inherits(x$prob_dist, "distcrete")) {
-    out <- x$prob_dist$p(q)
+  unlist <- length(x$prob_distribution) == 1
+  if (inherits(x$prob_distribution, "distcrete")) {
+    out <- x$prob_distribution$p(q)
     if (log) out <- log(out)
   } else {
-    out <- distributional::cdf(x$prob_dist, q = q, ..., log = log)
+    out <- distributional::cdf(x$prob_distribution, q = q, ..., log = log)
   }
   out <- if (unlist) unlist(out, recursive = FALSE) else out
   out
@@ -603,11 +575,11 @@ quantile.epiparameter <- function(x, p, ...) {
   if (isFALSE(is_parameterised(x))) {
     stop("<epiparameter> is unparameterised", call. = FALSE)
   }
-  unlist <- length(x$prob_dist) == 1
-  if (inherits(x$prob_dist, "distcrete")) {
-    out <- x$prob_dist$q(p)
+  unlist <- length(x$prob_distribution) == 1
+  if (inherits(x$prob_distribution, "distcrete")) {
+    out <- x$prob_distribution$q(p)
   } else {
-    out <- stats::quantile(x$prob_dist, p = p)
+    out <- stats::quantile(x$prob_distribution, p = p)
   }
   out <- if (unlist) unlist(out, recursive = FALSE) else out
   out
@@ -626,13 +598,13 @@ generate.epiparameter <- function(x, times, ...) {
   }
   # check times is a single number for consistent behaviour
   checkmate::assert_number(times)
-  if (inherits(x$prob_dist, "distcrete")) {
-    unlist <- length(x$prob_dist) == 1
-    out <- x$prob_dist$r(n = times)
+  if (inherits(x$prob_distribution, "distcrete")) {
+    unlist <- length(x$prob_distribution) == 1
+    out <- x$prob_distribution$r(n = times)
     out <- if (unlist) unlist(out, recursive = FALSE) else out
   } else {
-    recursive <- length(x$prob_dist) == 1
-    out <- distributional::generate(x$prob_dist, times = times)
+    recursive <- length(x$prob_distribution) == 1
+    out <- distributional::generate(x$prob_distribution, times = times)
     out <- unlist(out, recursive = recursive)
   }
   out
@@ -654,8 +626,10 @@ generate.epiparameter <- function(x, times, ...) {
 #' ebola_incubation <- epiparameter(
 #'   disease = "ebola",
 #'   epi_dist = "incubation_period",
-#'   prob_distribution = "gamma",
-#'   prob_distribution_params = c(shape = 1, scale = 1)
+#'   prob_distribution = create_prob_distribution(
+#'     prob_distribution = "gamma",
+#'     prob_distribution_params = c(shape = 1, scale = 1)
+#'   )
 #' )
 #' discretise(ebola_incubation)
 discretise <- function(x, ...) {
@@ -670,7 +644,7 @@ discretize <- discretise
 #' @export
 discretise.epiparameter <- function(x, ...) {
   # check if distribution is already discretised if so return early
-  if (inherits(x$prob_dist, "distcrete")) {
+  if (inherits(x$prob_distribution, "distcrete")) {
     message("Distribution in `epiparameter` is already discretised")
     return(x)
   } else {
@@ -695,7 +669,7 @@ discretise.epiparameter <- function(x, ...) {
 
       # trunc dist family is truncated so get prob dist by unclassing dist and
       # extracting name
-      list_dist <- unclass(x$prob_dist)
+      list_dist <- unclass(x$prob_distribution)
       prob_dist <- gsub(
         pattern = "dist_",
         replacement = "",
@@ -706,14 +680,14 @@ discretise.epiparameter <- function(x, ...) {
 
     # standardise distribution parameter names
     prob_dist_params <- .clean_params(
-      prob_dist = prob_dist,
-      prob_dist_params = prob_dist_params
+      prob_distribution = prob_dist,
+      prob_distribution_params = prob_dist_params
     )
 
     # create a new discretised probability distribution
-    x$prob_dist <- create_prob_dist(
-      prob_dist = prob_dist,
-      prob_dist_params = prob_dist_params,
+    x$prob_distribution <- create_prob_distribution(
+      prob_distribution = prob_dist,
+      prob_distribution_params = prob_dist_params,
       discretise = TRUE,
       truncation = NA
     )
@@ -749,8 +723,10 @@ discretise.default <- function(x, ...) {
 #' ep <- epiparameter(
 #'   disease = "ebola",
 #'   epi_dist = "incubation_period",
-#'   prob_distribution = "gamma",
-#'   prob_distribution_params = c(shape = 1, scale = 1)
+#'   prob_distribution = create_prob_distribution(
+#'     prob_distribution = "gamma",
+#'     prob_distribution_params = c(shape = 1, scale = 1)
+#'   )
 #' )
 #' family(ep)
 #'
@@ -758,27 +734,29 @@ discretise.default <- function(x, ...) {
 #' ep <- epiparameter(
 #'   disease = "ebola",
 #'   epi_dist = "incubation_period",
-#'   prob_distribution = "lnorm",
-#'   prob_distribution_params = c(meanlog = 1, sdlog = 1),
-#'   discretise = TRUE
+#'   prob_distribution = create_prob_distribution(
+#'     prob_distribution = "lnorm",
+#'     prob_distribution_params = c(meanlog = 1, sdlog = 1),
+#'     discretise = TRUE
+#'   )
 #' )
 #' family(ep)
 family.epiparameter <- function(object, ...) {
-  if (inherits(object$prob_dist, "distcrete")) {
-    prob_dist <- object$prob_dist$name
-  } else if (inherits(object$prob_dist, "distribution")) {
+  if (inherits(object$prob_distribution, "distcrete")) {
+    prob_dist <- object$prob_distribution$name
+  } else if (inherits(object$prob_distribution, "distribution")) {
     if (is_truncated(object)) {
       prob_dist <- gsub(
         pattern = "dist_",
         replacement = "",
-        x = class(unclass(unclass(object$prob_dist)[[1]])[[1]])[1],
+        x = class(unclass(unclass(object$prob_distribution)[[1]])[[1]])[1],
         fixed = TRUE
       )
     } else {
-      prob_dist <- stats::family(object$prob_dist)
+      prob_dist <- stats::family(object$prob_distribution)
     }
-  } else if (is.character(object$prob_dist)) {
-    prob_dist <- object$prob_dist
+  } else if (is.character(object$prob_distribution)) {
+    prob_dist <- object$prob_distribution
   } else {
     return(NA)
   }
@@ -814,17 +792,21 @@ family.epiparameter <- function(object, ...) {
 #' ep <- epiparameter(
 #'   disease = "ebola",
 #'   epi_dist = "incubation_period",
-#'   prob_distribution = "lnorm",
-#'   prob_distribution_params = c(meanlog = 1, sdlog = 1)
+#'   prob_distribution = create_prob_distribution(
+#'     prob_distribution = "lnorm",
+#'     prob_distribution_params = c(meanlog = 1, sdlog = 1)
+#'   )
 #' )
 #' is_truncated(ep)
 #'
 #' ep <- epiparameter(
 #'   disease = "ebola",
 #'   epi_dist = "incubation_period",
-#'   prob_distribution = "lnorm",
-#'   prob_distribution_params = c(meanlog = 1, sdlog = 1),
-#'   truncation = 10
+#'   prob_distribution = create_prob_distribution(
+#'     prob_distribution = "lnorm",
+#'     prob_distribution_params = c(meanlog = 1, sdlog = 1),
+#'     truncation = 10
+#'   )
 #' )
 #' is_truncated(ep)
 is_truncated <- function(x) {
@@ -834,18 +816,18 @@ is_truncated <- function(x) {
   )
 
   # distcrete distributions cannot be truncated
-  if (inherits(x$prob_dist, "distcrete")) {
+  if (inherits(x$prob_distribution, "distcrete")) {
     return(FALSE)
   }
 
   # unparameterised objects cannot be truncated
   # dont use is_parameterised due to infinite recursion
-  if (is.na(x$prob_dist) || is.character(x$prob_dist)) {
+  if (is.na(x$prob_distribution) || is.character(x$prob_distribution)) {
     return(FALSE)
   }
 
   # use stats::family instead of epiparameter::family to check truncated
-  if (identical(stats::family(x$prob_dist), "truncated")) {
+  if (identical(stats::family(x$prob_distribution), "truncated")) {
     return(TRUE)
   } else {
     return(FALSE)
@@ -869,8 +851,10 @@ is_truncated <- function(x) {
 #' ep <- epiparameter(
 #'   disease = "ebola",
 #'   epi_dist = "incubation_period",
-#'   prob_distribution = "lnorm",
-#'   prob_distribution_params = c(meanlog = 1, sdlog = 1)
+#'   prob_distribution = create_prob_distribution(
+#'     prob_distribution = "lnorm",
+#'     prob_distribution_params = c(meanlog = 1, sdlog = 1)
+#'   )
 #' )
 #' is_continuous(ep)
 #' is_continuous(discretise(ep))
@@ -878,8 +862,10 @@ is_truncated <- function(x) {
 #' ep <- epiparameter(
 #'   disease = "ebola",
 #'   epi_dist = "offspring distribution",
-#'   prob_distribution = "nbinom",
-#'   prob_distribution_params = c(mean = 2, dispersion = 0.5)
+#'   prob_distribution = create_prob_distribution(
+#'     prob_distribution = "nbinom",
+#'     prob_distribution_params = c(mean = 2, dispersion = 0.5)
+#'   )
 #' )
 #' is_continuous(ep)
 is_continuous <- function(x) {
@@ -888,7 +874,7 @@ is_continuous <- function(x) {
       is_epiparameter(x)
   )
   family(x) %in% c("gamma", "lnorm", "weibull", "normal") &&
-    !inherits(x$prob_dist, "distcrete")
+    !inherits(x$prob_distribution, "distcrete")
 }
 
 #' Mean method for `<epiparameter>` class
diff --git a/R/epiparameter_db.R b/R/epiparameter_db.R
index c75b312a1..86a4555b9 100644
--- a/R/epiparameter_db.R
+++ b/R/epiparameter_db.R
@@ -559,6 +559,17 @@ epidist_db <- function(disease = "all",
   params <- params_uncertainty$params
   uncertainty <- params_uncertainty$uncertainty
 
+  # TODO: remove this and move precision from parameter to summary statistic
+  if ("precision" %in% names(params)) {
+    if (length(params) == 1) {
+      params <- NA_real_
+      uncertainty <- list(uncertainty = create_uncertainty())
+    } else {
+      params <- params[-which(names(params) == "precision")]
+      uncertainty <- uncertainty[-which(names(uncertainty) == "precision")]
+    }
+  }
+
   # format summary statistics
   ss <- lapply(x$summary_statistics, unlist)
   names(ss$quantile_values) <- ss$quantile_names
@@ -604,16 +615,18 @@ epidist_db <- function(disease = "all",
     disease = x$disease,
     pathogen = x$pathogen,
     epi_dist = x$epi_distribution,
-    prob_distribution = x$probability_distribution$prob_distribution,
-    prob_distribution_params = params,
+    prob_distribution = create_prob_distribution(
+      prob_distribution = x$probability_distribution$prob_distribution,
+      prob_distribution_params = params,
+      discretise = discretised,
+      truncation = truncation
+    ),
     uncertainty = uncertainty,
     summary_stats = ss,
     auto_calc_params = TRUE,
     citation = cit,
     metadata = meta,
     method_assess = method,
-    discretise = discretised,
-    truncation = truncation,
     notes = x$notes
   )
 }
@@ -633,7 +646,7 @@ epidist_db <- function(disease = "all",
     return(
       list(
         params = NA_real_,
-        uncertainty = create_uncertainty()
+        uncertainty = list(uncertainty = create_uncertainty())
       )
     )
   }
@@ -717,7 +730,8 @@ epidist_db <- function(disease = "all",
   stopifnot(is_epiparameter(lst))
   if (nse_subject == "prob_dist") {
     # special case to subset by dist as name is extracted first
-    if (is.object(lst$prob_dist) || is.character(lst$prob_dist)) {
+    if (is.object(lst$prob_distribution) ||
+        is.character(lst$prob_distribution)) {
       prob_dist <- family(lst) # nolint prob_dist used in eval
       eval(expr = condition)
     } else {
diff --git a/R/parameter_tbl.R b/R/parameter_tbl.R
index 01e04c696..62012cc0c 100644
--- a/R/parameter_tbl.R
+++ b/R/parameter_tbl.R
@@ -58,10 +58,10 @@ parameter_tbl <- function(multi_epiparameter,
   )
   prob_dist <- vapply(
     multi_epiparameter, function(x) {
-      switch(class(x$prob_dist)[1],
+      switch(class(x$prob_distribution)[1],
         distcrete = family(x),
         distribution = family(x),
-        character = x$prob_dist,
+        character = x$prob_distribution,
         logical = NA_character_
       )
     },
diff --git a/R/plot.R b/R/plot.R
index 38836a0f9..bea286975 100644
--- a/R/plot.R
+++ b/R/plot.R
@@ -26,8 +26,10 @@
 #' ep <- epiparameter(
 #'   disease = "ebola",
 #'   epi_dist = "incubation_period",
-#'   prob_distribution = "gamma",
-#'   prob_distribution_params = c(shape = 2, scale = 1)
+#'   prob_distribution = create_prob_distribution(
+#'     prob_distribution = "gamma",
+#'     prob_distribution_params = c(shape = 2, scale = 1)
+#'   )
 #' )
 #' plot(ep)
 #'
diff --git a/R/sysdata.rda b/R/sysdata.rda
index 8699d8125..203610a8f 100644
Binary files a/R/sysdata.rda and b/R/sysdata.rda differ
diff --git a/R/utils.R b/R/utils.R
index d3dc5f656..616260503 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -68,3 +68,5 @@ calc_disc_dist_quantile <- function(prob, days, quantile) {
   cit <- paste0(cit, " (", x$year, ")")
   cit
 }
+
+`%||%` <- function(x, y) if (is.null(x)) y else x
diff --git a/man/create_prob_dist.Rd b/man/create_prob_distribution.Rd
similarity index 69%
rename from man/create_prob_dist.Rd
rename to man/create_prob_distribution.Rd
index a180d724a..3c9adb308 100644
--- a/man/create_prob_dist.Rd
+++ b/man/create_prob_distribution.Rd
@@ -1,29 +1,29 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/create_prob_dist.R
-\name{create_prob_dist}
-\alias{create_prob_dist}
+% Please edit documentation in R/create_prob_distribution.R
+\name{create_prob_distribution}
+\alias{create_prob_distribution}
 \title{Create a distribution object}
 \usage{
-create_prob_dist(
-  prob_dist,
-  prob_dist_params,
+create_prob_distribution(
+  prob_distribution,
+  prob_distribution_params,
   discretise = FALSE,
   truncation = NA,
   ...
 )
 }
 \arguments{
-\item{prob_dist}{A character string specifying the probability
-distribution. This should match the R naming convention of probability
+\item{prob_distribution}{A \code{character} string specifying the probability
+distribution. This should match the \R naming convention of probability
 distributions (e.g. lognormal is \code{lnorm}, negative binomial is \code{nbinom}, and
 geometric is \code{geom}).}
 
-\item{prob_dist_params}{A named vector of probability distribution
+\item{prob_distribution_params}{A named vector of probability distribution
 parameters.}
 
 \item{discretise}{A boolean \code{logical} whether the distribution is
-discretised.
-Default is FALSE which assumes a continuous probability distribution}
+discretised. Default is \code{FALSE} which assumes a continuous probability
+distribution.}
 
 \item{truncation}{A \code{numeric} specifying the truncation point if the inferred
 distribution was truncated, \code{NA} if not or unknown.}
@@ -37,7 +37,8 @@ constructor function, e.g. for the Gamma distribution see
 \code{\link[distributional:dist_gamma]{distributional::dist_gamma()}}.}
 }
 \value{
-An S3 class containing the probability distribution.
+An S3 class containing the probability distribution or a \code{character}
+string if the parameters of the probability distribution are unknown.
 }
 \description{
 Creates an S3 class holding the distribution and parameters
@@ -65,33 +66,33 @@ see \code{\link[distcrete:distcrete]{distcrete::distcrete()}}.
 }
 \examples{
 # example with continuous distribution without truncation
-create_prob_dist(
-  prob_dist = "gamma",
-  prob_dist_params = c(shape = 1, scale = 1),
+create_prob_distribution(
+  prob_distribution = "gamma",
+  prob_distribution_params = c(shape = 1, scale = 1),
   discretise = FALSE,
   truncation = NA
 )
 
 # example with continuous distribution with truncation
-create_prob_dist(
-  prob_dist = "gamma",
-  prob_dist_params = c(shape = 1, scale = 1),
+create_prob_distribution(
+  prob_distribution = "gamma",
+  prob_distribution_params = c(shape = 1, scale = 1),
   discretise = FALSE,
   truncation = 10
 )
 
 # example with discrete distribution
-create_prob_dist(
-  prob_dist = "gamma",
-  prob_dist_params = c(shape = 1, scale = 1),
+create_prob_distribution(
+  prob_distribution = "gamma",
+  prob_distribution_params = c(shape = 1, scale = 1),
   discretise = TRUE,
   truncation = NA
 )
 
 # example passing extra arguments to distcrete
-create_prob_dist(
-  prob_dist = "gamma",
-  prob_dist_params = c(shape = 1, scale = 1),
+create_prob_distribution(
+  prob_distribution = "gamma",
+  prob_distribution_params = c(shape = 1, scale = 1),
   discretise = TRUE,
   truncation = NA,
   w = 0.5
diff --git a/man/discretise.Rd b/man/discretise.Rd
index bb0b25ae1..878728996 100644
--- a/man/discretise.Rd
+++ b/man/discretise.Rd
@@ -35,8 +35,10 @@ discretised distribution (using an object from the \code{{distcrete}} package).
 ebola_incubation <- epiparameter(
   disease = "ebola",
   epi_dist = "incubation_period",
-  prob_distribution = "gamma",
-  prob_distribution_params = c(shape = 1, scale = 1)
+  prob_distribution = create_prob_distribution(
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 1, scale = 1)
+  )
 )
 discretise(ebola_incubation)
 }
diff --git a/man/dot-calc_dist_params.Rd b/man/dot-calc_dist_params.Rd
index bd7cf8783..ea6f8f14f 100644
--- a/man/dot-calc_dist_params.Rd
+++ b/man/dot-calc_dist_params.Rd
@@ -5,16 +5,14 @@
 \title{Calculate the parameters of a probability distribution from a list of
 summary statistics}
 \usage{
-.calc_dist_params(prob_dist, prob_dist_params, summary_stats, sample_size)
+.calc_dist_params(prob_distribution, summary_stats, sample_size)
 }
 \arguments{
-\item{prob_dist}{A character string specifying the probability
-distribution. This should match the R naming convention of probability
-distributions (e.g. lognormal is \code{lnorm}, negative binomial is \code{nbinom}, and
-geometric is \code{geom}).}
-
-\item{prob_dist_params}{A named vector of probability distribution
-parameters.}
+\item{prob_distribution}{An S3 class containing the probability
+distribution or a character string if the parameters of the probability
+distribution are unknown but the name of the distribution is known, or \code{NA}
+if the distribution name and parameters are unknown. Use
+\code{\link[=create_prob_distribution]{create_prob_distribution()}} to create \code{prob_distribution}.}
 
 \item{summary_stats}{A list of summary statistics, use
 \code{\link[=create_summary_stats]{create_summary_stats()}} to create list. This list can include
diff --git a/man/dot-clean_params.Rd b/man/dot-clean_params.Rd
index c024d6245..409aeb3d6 100644
--- a/man/dot-clean_params.Rd
+++ b/man/dot-clean_params.Rd
@@ -11,7 +11,7 @@
 \alias{.clean_params_exp}
 \title{Standardise distribution parameters}
 \usage{
-.clean_params(prob_dist, prob_dist_params)
+.clean_params(prob_distribution, prob_distribution_params)
 
 .clean_params_gamma(prob_dist_params)
 
@@ -28,12 +28,12 @@
 .clean_params_exp(prob_dist_params)
 }
 \arguments{
-\item{prob_dist}{A character string specifying the probability
-distribution. This should match the R naming convention of probability
+\item{prob_distribution}{A \code{character} string specifying the probability
+distribution. This should match the \R naming convention of probability
 distributions (e.g. lognormal is \code{lnorm}, negative binomial is \code{nbinom}, and
 geometric is \code{geom}).}
 
-\item{prob_dist_params}{A named vector of probability distribution
+\item{prob_distribution_params}{A named vector of probability distribution
 parameters.}
 }
 \value{
diff --git a/man/dot-clean_uncertainty.Rd b/man/dot-clean_uncertainty.Rd
new file mode 100644
index 000000000..d3a540656
--- /dev/null
+++ b/man/dot-clean_uncertainty.Rd
@@ -0,0 +1,27 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/epiparameter-utils.R
+\name{.clean_uncertainty}
+\alias{.clean_uncertainty}
+\title{Standardise distribution parameter uncertainty}
+\usage{
+.clean_uncertainty(x, prob_distribution, uncertainty_missing)
+}
+\arguments{
+\item{x}{An \verb{<epiparameter>} object.}
+
+\item{prob_distribution}{An S3 class containing the probability
+distribution or a character string if the parameters of the probability
+distribution are unknown but the name of the distribution is known, or \code{NA}
+if the distribution name and parameters are unknown. Use
+\code{\link[=create_prob_distribution]{create_prob_distribution()}} to create \code{prob_distribution}.}
+
+\item{uncertainty_missing}{A boolean \code{logical} as to whether \code{uncertainty}
+is missing (see \code{\link[=missing]{missing()}}) from the parent function.}
+}
+\value{
+An uncertainty list for an \verb{<epiparameter>} object.
+}
+\description{
+Standardise distribution parameter uncertainty
+}
+\keyword{internal}
diff --git a/man/epiparameter.Rd b/man/epiparameter.Rd
index 67c23cd50..3474b0fe4 100644
--- a/man/epiparameter.Rd
+++ b/man/epiparameter.Rd
@@ -8,17 +8,14 @@ epiparameter(
   disease,
   pathogen = NA_character_,
   epi_dist,
-  prob_distribution = NA_character_,
-  prob_distribution_params = NA_real_,
+  prob_distribution = create_prob_distribution(prob_distribution = NA_character_),
   uncertainty = create_uncertainty(),
   summary_stats = create_summary_stats(),
-  auto_calc_params = TRUE,
   citation = create_citation(),
   metadata = create_metadata(),
   method_assess = create_method_assess(),
-  discretise = FALSE,
-  truncation = NA_real_,
   notes = NULL,
+  auto_calc_params = TRUE,
   ...
 )
 }
@@ -31,13 +28,11 @@ disease, or \code{NA} if not known.}
 \item{epi_dist}{A \code{character} string with the name of the
 epidemiological distribution type.}
 
-\item{prob_distribution}{A \code{character} string specifying the probability
-distribution. This should match the \R naming convention of probability
-distributions (e.g. lognormal is \code{lnorm}, negative binomial is \code{nbinom}, and
-geometric is \code{geom}).}
-
-\item{prob_distribution_params}{A named vector of probability distribution
-parameters.}
+\item{prob_distribution}{An S3 class containing the probability
+distribution or a character string if the parameters of the probability
+distribution are unknown but the name of the distribution is known, or \code{NA}
+if the distribution name and parameters are unknown. Use
+\code{\link[=create_prob_distribution]{create_prob_distribution()}} to create \code{prob_distribution}.}
 
 \item{uncertainty}{A list of named vectors with the uncertainty around
 the probability distribution parameters. If uncertainty around the parameter
@@ -53,13 +48,6 @@ data used to fit the distribution such as lower and upper range. The summary
 statistics can also include uncertainty around metrics such as confidence
 interval around mean and standard deviation.}
 
-\item{auto_calc_params}{A boolean \code{logical} determining whether to try and
-calculate the probability distribution parameters from summary statistics if
-distribution parameters are not provided. Default is \code{TRUE}. In the case when
-sufficient summary statistics are provided and the parameter(s) of the
-distribution are not, the \code{\link[=.calc_dist_params]{.calc_dist_params()}} function is called to
-calculate the parameters and add them to the \code{epiparameter} object created.}
-
 \item{citation}{A \verb{<bibentry>} with the citation of the source of the
 data or the paper that inferred the distribution parameters, use
 \code{\link[=create_citation]{create_citation()}} to create citation.}
@@ -76,16 +64,16 @@ delay distribution such as extrinsic incubation period) unless
 the distribution, use \code{\link[=create_method_assess]{create_method_assess()}} to create method
 assessment.}
 
-\item{discretise}{A boolean \code{logical} whether the distribution is
-discretised.
-Default is FALSE which assumes a continuous probability distribution}
-
-\item{truncation}{A \code{numeric} specifying the truncation point if the inferred
-distribution was truncated, \code{NA} if not or unknown.}
-
 \item{notes}{A \code{character} string with any additional information about the
 data, inference method or disease.}
 
+\item{auto_calc_params}{A boolean \code{logical} determining whether to try and
+calculate the probability distribution parameters from summary statistics if
+distribution parameters are not provided. Default is \code{TRUE}. In the case when
+sufficient summary statistics are provided and the parameter(s) of the
+distribution are not, the \code{\link[=.calc_dist_params]{.calc_dist_params()}} function is called to
+calculate the parameters and add them to the \code{epiparameter} object created.}
+
 \item{...}{\link{dots} Extra arguments to be passed to internal functions.
 
 This is most commonly used to pass arguments to \code{\link[distcrete:distcrete]{distcrete::distcrete()}}
@@ -131,17 +119,21 @@ Accepted \verb{<epiparameter>} distribution parameterisations are:
 ebola_incubation <- epiparameter(
   disease = "ebola",
   epi_dist = "incubation_period",
-  prob_distribution = "gamma",
-  prob_distribution_params = c(shape = 1, scale = 1)
+  prob_distribution = create_prob_distribution(
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 1, scale = 1)
+  )
 )
 
 # minimal input required for discrete `epiparameter`
 ebola_incubation <- epiparameter(
   disease = "ebola",
   epi_dist = "incubation_period",
-  prob_distribution = "gamma",
-  prob_distribution_params = c(shape = 1, scale = 1),
-  discretise = TRUE
+  prob_distribution = create_prob_distribution(
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 1, scale = 1),
+    discretise = TRUE
+  )
 )
 
 # example with more fields filled in
@@ -149,9 +141,16 @@ ebola_incubation <- epiparameter(
   disease = "ebola",
   pathogen = "ebola_virus",
   epi_dist = "incubation",
-  prob_distribution = "gamma",
-  prob_distribution_params = c(shape = 1, scale = 1),
-  uncertainty = create_uncertainty(),
+  prob_distribution = create_prob_distribution(
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 1, scale = 1),
+    discretise = FALSE,
+    truncation = NA
+  ),
+  uncertainty = list(
+    shape = create_uncertainty(),
+    scale = create_uncertainty()
+  ),
   summary_stats = create_summary_stats(mean = 2, sd = 1),
   citation = create_citation(
     author = person(given = "John", family = "Smith"),
@@ -169,8 +168,6 @@ ebola_incubation <- epiparameter(
   method_assess = create_method_assess(
     censored = TRUE
   ),
-  discretise = FALSE,
-  truncation = NA,
   notes = "No notes"
 )
 }
diff --git a/man/epiparameter_distribution_functions.Rd b/man/epiparameter_distribution_functions.Rd
index 89b5c3296..d2c4e06c9 100644
--- a/man/epiparameter_distribution_functions.Rd
+++ b/man/epiparameter_distribution_functions.Rd
@@ -46,8 +46,10 @@ functions. These operate on any distribution that can be included in an
 ep <- epiparameter(
   disease = "ebola",
   epi_dist = "incubation_period",
-  prob_distribution = "gamma",
-  prob_distribution_params = c(shape = 1, scale = 1)
+  prob_distribution = create_prob_distribution(
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 1, scale = 1)
+  )
 )
 
 # example of each distribution method for an `epiparameter` object
diff --git a/man/family.epiparameter.Rd b/man/family.epiparameter.Rd
index 0417b2251..ab14c2e7e 100644
--- a/man/family.epiparameter.Rd
+++ b/man/family.epiparameter.Rd
@@ -26,8 +26,10 @@ output irrespective of the internal distribution class.
 ep <- epiparameter(
   disease = "ebola",
   epi_dist = "incubation_period",
-  prob_distribution = "gamma",
-  prob_distribution_params = c(shape = 1, scale = 1)
+  prob_distribution = create_prob_distribution(
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 1, scale = 1)
+  )
 )
 family(ep)
 
@@ -35,9 +37,11 @@ family(ep)
 ep <- epiparameter(
   disease = "ebola",
   epi_dist = "incubation_period",
-  prob_distribution = "lnorm",
-  prob_distribution_params = c(meanlog = 1, sdlog = 1),
-  discretise = TRUE
+  prob_distribution = create_prob_distribution(
+    prob_distribution = "lnorm",
+    prob_distribution_params = c(meanlog = 1, sdlog = 1),
+    discretise = TRUE
+  )
 )
 family(ep)
 }
diff --git a/man/format.epiparameter.Rd b/man/format.epiparameter.Rd
index 930f74db3..07bf8e618 100644
--- a/man/format.epiparameter.Rd
+++ b/man/format.epiparameter.Rd
@@ -22,8 +22,10 @@ Format method for \verb{<epiparameter>} class
 epiparameter <- epiparameter(
   disease = "ebola",
   epi_dist = "incubation_period",
-  prob_distribution = "gamma",
-  prob_distribution_params = c(shape = 1, scale = 1)
+  prob_distribution = create_prob_distribution(
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 1, scale = 1)
+  )
 )
 format(epiparameter)
 }
diff --git a/man/is_continuous.Rd b/man/is_continuous.Rd
index 0ce293433..8743e0f5a 100644
--- a/man/is_continuous.Rd
+++ b/man/is_continuous.Rd
@@ -26,8 +26,10 @@ binomial), and all \verb{<distcrete>} objects are discrete.
 ep <- epiparameter(
   disease = "ebola",
   epi_dist = "incubation_period",
-  prob_distribution = "lnorm",
-  prob_distribution_params = c(meanlog = 1, sdlog = 1)
+  prob_distribution = create_prob_distribution(
+    prob_distribution = "lnorm",
+    prob_distribution_params = c(meanlog = 1, sdlog = 1)
+  )
 )
 is_continuous(ep)
 is_continuous(discretise(ep))
@@ -35,8 +37,10 @@ is_continuous(discretise(ep))
 ep <- epiparameter(
   disease = "ebola",
   epi_dist = "offspring distribution",
-  prob_distribution = "nbinom",
-  prob_distribution_params = c(mean = 2, dispersion = 0.5)
+  prob_distribution = create_prob_distribution(
+    prob_distribution = "nbinom",
+    prob_distribution_params = c(mean = 2, dispersion = 0.5)
+  )
 )
 is_continuous(ep)
 }
diff --git a/man/is_epiparameter.Rd b/man/is_epiparameter.Rd
index 489a357c0..37bca1ec2 100644
--- a/man/is_epiparameter.Rd
+++ b/man/is_epiparameter.Rd
@@ -20,8 +20,10 @@ Check object is an \verb{<epiparameter>}
 ep <- epiparameter(
   disease = "ebola",
   epi_dist = "serial_interval",
-  prob_distribution = "gamma",
-  prob_distribution_params = c(shape = 1, scale = 1)
+  prob_distribution = create_prob_distribution(
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 1, scale = 1)
+  )
 )
 
 is_epiparameter(ep)
diff --git a/man/is_epiparameter_params.Rd b/man/is_epiparameter_params.Rd
index 34efa40af..817607b10 100644
--- a/man/is_epiparameter_params.Rd
+++ b/man/is_epiparameter_params.Rd
@@ -5,15 +5,15 @@
 \title{Check whether the vector of parameters for the probability distribution
 are in the set of possible parameters used in the epiparameter package}
 \usage{
-is_epiparameter_params(prob_dist, prob_dist_params)
+is_epiparameter_params(prob_distribution, prob_distribution_params)
 }
 \arguments{
-\item{prob_dist}{A character string specifying the probability
-distribution. This should match the R naming convention of probability
+\item{prob_distribution}{A \code{character} string specifying the probability
+distribution. This should match the \R naming convention of probability
 distributions (e.g. lognormal is \code{lnorm}, negative binomial is \code{nbinom}, and
 geometric is \code{geom}).}
 
-\item{prob_dist_params}{A named vector of probability distribution
+\item{prob_distribution_params}{A named vector of probability distribution
 parameters.}
 }
 \value{
@@ -29,7 +29,7 @@ is truncated or discretised.
 }
 \examples{
 is_epiparameter_params(
-  prob_dist = "gamma",
-  prob_dist_params = c(shape = 2, scale = 1)
+  prob_distribution = "gamma",
+  prob_distribution_params = c(shape = 2, scale = 1)
 )
 }
diff --git a/man/is_parameterised.Rd b/man/is_parameterised.Rd
index 3ca83db4c..c34eaaa06 100644
--- a/man/is_parameterised.Rd
+++ b/man/is_parameterised.Rd
@@ -31,8 +31,10 @@ distribution and distribution parameters
 ep <- epiparameter(
   disease = "ebola",
   epi_dist = "incubation",
-  prob_distribution = "gamma",
-  prob_distribution_params = c(shape = 1, scale = 1)
+  prob_distribution = create_prob_distribution(
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 1, scale = 1)
+  )
 )
 is_parameterised(ep)
 
diff --git a/man/is_truncated.Rd b/man/is_truncated.Rd
index d86c81f36..546784b00 100644
--- a/man/is_truncated.Rd
+++ b/man/is_truncated.Rd
@@ -26,17 +26,21 @@ return \code{FALSE} by default.
 ep <- epiparameter(
   disease = "ebola",
   epi_dist = "incubation_period",
-  prob_distribution = "lnorm",
-  prob_distribution_params = c(meanlog = 1, sdlog = 1)
+  prob_distribution = create_prob_distribution(
+    prob_distribution = "lnorm",
+    prob_distribution_params = c(meanlog = 1, sdlog = 1)
+  )
 )
 is_truncated(ep)
 
 ep <- epiparameter(
   disease = "ebola",
   epi_dist = "incubation_period",
-  prob_distribution = "lnorm",
-  prob_distribution_params = c(meanlog = 1, sdlog = 1),
-  truncation = 10
+  prob_distribution = create_prob_distribution(
+    prob_distribution = "lnorm",
+    prob_distribution_params = c(meanlog = 1, sdlog = 1),
+    truncation = 10
+  )
 )
 is_truncated(ep)
 }
diff --git a/man/new_epiparameter.Rd b/man/new_epiparameter.Rd
index 7719ff079..fbac046fc 100644
--- a/man/new_epiparameter.Rd
+++ b/man/new_epiparameter.Rd
@@ -8,17 +8,14 @@ new_epiparameter(
   disease = character(),
   pathogen = character(),
   epi_dist = character(),
-  prob_dist = list(),
-  prob_dist_params = numeric(),
+  prob_distribution = list(),
   uncertainty = list(),
   summary_stats = list(),
-  auto_calc_params = logical(),
   citation = character(),
   metadata = list(),
   method_assess = list(),
-  discretise = logical(),
-  truncation = numeric(),
   notes = character(),
+  auto_calc_params = logical(),
   ...
 )
 }
@@ -31,13 +28,11 @@ disease, or \code{NA} if not known.}
 \item{epi_dist}{A \code{character} string with the name of the
 epidemiological distribution type.}
 
-\item{prob_dist}{A character string specifying the probability
-distribution. This should match the R naming convention of probability
-distributions (e.g. lognormal is \code{lnorm}, negative binomial is \code{nbinom}, and
-geometric is \code{geom}).}
-
-\item{prob_dist_params}{A named vector of probability distribution
-parameters.}
+\item{prob_distribution}{An S3 class containing the probability
+distribution or a character string if the parameters of the probability
+distribution are unknown but the name of the distribution is known, or \code{NA}
+if the distribution name and parameters are unknown. Use
+\code{\link[=create_prob_distribution]{create_prob_distribution()}} to create \code{prob_distribution}.}
 
 \item{uncertainty}{A list of named vectors with the uncertainty around
 the probability distribution parameters. If uncertainty around the parameter
@@ -53,13 +48,6 @@ data used to fit the distribution such as lower and upper range. The summary
 statistics can also include uncertainty around metrics such as confidence
 interval around mean and standard deviation.}
 
-\item{auto_calc_params}{A boolean \code{logical} determining whether to try and
-calculate the probability distribution parameters from summary statistics if
-distribution parameters are not provided. Default is \code{TRUE}. In the case when
-sufficient summary statistics are provided and the parameter(s) of the
-distribution are not, the \code{\link[=.calc_dist_params]{.calc_dist_params()}} function is called to
-calculate the parameters and add them to the \code{epiparameter} object created.}
-
 \item{citation}{A \verb{<bibentry>} with the citation of the source of the
 data or the paper that inferred the distribution parameters, use
 \code{\link[=create_citation]{create_citation()}} to create citation.}
@@ -76,16 +64,16 @@ delay distribution such as extrinsic incubation period) unless
 the distribution, use \code{\link[=create_method_assess]{create_method_assess()}} to create method
 assessment.}
 
-\item{discretise}{A boolean \code{logical} whether the distribution is
-discretised.
-Default is FALSE which assumes a continuous probability distribution}
-
-\item{truncation}{A \code{numeric} specifying the truncation point if the inferred
-distribution was truncated, \code{NA} if not or unknown.}
-
 \item{notes}{A \code{character} string with any additional information about the
 data, inference method or disease.}
 
+\item{auto_calc_params}{A boolean \code{logical} determining whether to try and
+calculate the probability distribution parameters from summary statistics if
+distribution parameters are not provided. Default is \code{TRUE}. In the case when
+sufficient summary statistics are provided and the parameter(s) of the
+distribution are not, the \code{\link[=.calc_dist_params]{.calc_dist_params()}} function is called to
+calculate the parameters and add them to the \code{epiparameter} object created.}
+
 \item{...}{\link{dots} Extra arguments to be passed to internal functions.
 
 This is most commonly used to pass arguments to \code{\link[distcrete:distcrete]{distcrete::distcrete()}}
diff --git a/man/plot.epiparameter.Rd b/man/plot.epiparameter.Rd
index d9f32bfe8..49d8ac80a 100644
--- a/man/plot.epiparameter.Rd
+++ b/man/plot.epiparameter.Rd
@@ -32,8 +32,10 @@ first and last day to plot on the x-axis can be supplied to \code{xlim}
 ep <- epiparameter(
   disease = "ebola",
   epi_dist = "incubation_period",
-  prob_distribution = "gamma",
-  prob_distribution_params = c(shape = 2, scale = 1)
+  prob_distribution = create_prob_distribution(
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 2, scale = 1)
+  )
 )
 plot(ep)
 
diff --git a/man/print.epiparameter.Rd b/man/print.epiparameter.Rd
index ac3b04e75..263d4af00 100644
--- a/man/print.epiparameter.Rd
+++ b/man/print.epiparameter.Rd
@@ -21,8 +21,10 @@ Print method for \verb{<epiparameter>} class
 epiparameter <- epiparameter(
   disease = "ebola",
   epi_dist = "incubation_period",
-  prob_distribution = "gamma",
-  prob_distribution_params = c(shape = 1, scale = 1)
+  prob_distribution = create_prob_distribution(
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 1, scale = 1)
+  )
 )
 epiparameter
 }
diff --git a/tests/testthat/_snaps/epiparameter.md b/tests/testthat/_snaps/epiparameter.md
index 47936b5fb..a42605a8b 100644
--- a/tests/testthat/_snaps/epiparameter.md
+++ b/tests/testthat/_snaps/epiparameter.md
@@ -1,8 +1,8 @@
 # epiparameter.print works as expected
 
     Code
-      epiparameter(disease = "ebola", epi_dist = "incubation", prob_distribution = "gamma",
-        prob_distribution_params = c(shape = 1, scale = 1))
+      epiparameter(disease = "ebola", epi_dist = "incubation", prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma", prob_distribution_params = c(shape = 1, scale = 1)))
     Message
       Citation cannot be created as author, year, journal or title is missing
     Output
@@ -18,8 +18,9 @@
 ---
 
     Code
-      epiparameter(disease = "ebola", epi_dist = "incubation", prob_distribution = "gamma",
-        prob_distribution_params = c(shape = 1, scale = 1), discretise = TRUE)
+      epiparameter(disease = "ebola", epi_dist = "incubation", prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma", prob_distribution_params = c(shape = 1, scale = 1),
+        discretise = TRUE))
     Message
       Citation cannot be created as author, year, journal or title is missing
     Output
diff --git a/tests/testthat/_snaps/print.md b/tests/testthat/_snaps/print.md
index e1f8c93bb..98820ff0b 100644
--- a/tests/testthat/_snaps/print.md
+++ b/tests/testthat/_snaps/print.md
@@ -23,8 +23,8 @@
       <https://doi.org/10.1016/S1473-3099%2809%2970069-6>.
       Distribution: lnorm
       Parameters:
-        meanlog: 1.247
-        sdlog: 0.975
+        meanlog: 1.723
+        sdlog: 0.231
       
       [[2]]
       Disease: Human Coronavirus
@@ -37,8 +37,8 @@
       <https://doi.org/10.1016/S1473-3099%2809%2970069-7>.
       Distribution: lnorm
       Parameters:
-        meanlog: 0.742
-        sdlog: 0.918
+        meanlog: 1.163
+        sdlog: 0.140
       
       [[3]]
       Disease: SARS
@@ -51,8 +51,8 @@
       <https://doi.org/10.1016/S1473-3099%2809%2970069-8>.
       Distribution: lnorm
       Parameters:
-        meanlog: 0.660
-        sdlog: 1.205
+        meanlog: 1.386
+        sdlog: 0.593
       
       # i 122 more elements
       # i Use `print(n = ...)` to see more elements.
diff --git a/tests/testthat/test-accessors.R b/tests/testthat/test-accessors.R
index 431dcbd1a..5c7dd94b4 100644
--- a/tests/testthat/test-accessors.R
+++ b/tests/testthat/test-accessors.R
@@ -11,7 +11,7 @@ test_that("get_parameters works as expected for unparameterised epiparameter", {
     epiparameter(
       disease = "Ebola",
       epi_dist = "incubation period",
-      prob_distribution = "gamma"
+      prob_distribution = create_prob_distribution(prob_distribution = "gamma")
     )
   )
   expect_true(is.na(get_parameters(ep)))
@@ -22,8 +22,10 @@ test_that("get_parameters works as expected for continuous epiparameter", {
     epiparameter(
       disease = "Ebola",
       epi_dist = "incubation period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1)
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1)
+      )
     )
   )
   params <- get_parameters(ep)
@@ -36,9 +38,11 @@ test_that("get_parameters works as expected for discretised epiparameter", {
     epiparameter(
       disease = "Ebola",
       epi_dist = "incubation period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1),
-      discretise = TRUE
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1),
+        discretise = TRUE
+      )
     )
   )
   params <- get_parameters(ep)
@@ -66,8 +70,10 @@ test_that("get_citation works as expected for manual epiparameter", {
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1),
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1)
+      ),
       citation = create_citation(
         author = person(given = "John F.", family = "Smith"),
         year = 2000,
@@ -86,8 +92,10 @@ test_that("get_citation works as expected for epiparameter missing citation", {
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1)
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1)
+      )
     )
   )
   citation <- get_citation(ep)
diff --git a/tests/testthat/test-calc_dist_params.R b/tests/testthat/test-calc_dist_params.R
index 72e0c6588..5b9187b69 100644
--- a/tests/testthat/test-calc_dist_params.R
+++ b/tests/testthat/test-calc_dist_params.R
@@ -1,7 +1,6 @@
 test_that(".calc_dist_params works as expected converting from mean and sd", {
   params <- .calc_dist_params(
-    prob_dist = "gamma",
-    prob_dist_params = NA,
+    prob_distribution = "gamma",
     summary_stats = create_summary_stats(mean = 5, sd = 2),
     sample_size = NA
   )
@@ -12,8 +11,7 @@ test_that(".calc_dist_params works as expected converting from mean and sd", {
 
 test_that(".calc_dist_params works as expected converting for different dist", {
   params <- .calc_dist_params(
-    prob_dist = "lnorm",
-    prob_dist_params = NA,
+    prob_distribution = "lnorm",
     summary_stats = create_summary_stats(mean = 5, sd = 2),
     sample_size = NA
   )
@@ -25,8 +23,7 @@ test_that(".calc_dist_params works as expected converting for different dist", {
 test_that(".calc_dist_params works as expected extracting from percentiles", {
   # messages for numerical optimisation suppressed
   params <- suppressMessages(.calc_dist_params(
-    prob_dist = "gamma",
-    prob_dist_params = NA,
+    prob_distribution = "gamma",
     summary_stats = create_summary_stats(
       quantiles = c("25" = 10, "75" = 20)
     ),
@@ -38,8 +35,7 @@ test_that(".calc_dist_params works as expected extracting from percentiles", {
 
   # messages for numerical optimisation suppressed
   params <- suppressMessages(.calc_dist_params(
-    prob_dist = "lnorm",
-    prob_dist_params = NA,
+    prob_distribution = "lnorm",
     summary_stats = create_summary_stats(
       quantiles = c("25" = 10, "75" = 20)
     ),
@@ -53,8 +49,7 @@ test_that(".calc_dist_params works as expected extracting from percentiles", {
 test_that(".calc_dist_params works as expected extracting from median & range", {
   # messages for numerical optimisation suppressed
   params <- suppressMessages(.calc_dist_params(
-    prob_dist = "gamma",
-    prob_dist_params = NA,
+    prob_distribution = "gamma",
     summary_stats = create_summary_stats(
       median = 10,
       lower_range = 5,
@@ -70,8 +65,7 @@ test_that(".calc_dist_params works as expected extracting from median & range",
 test_that(".calc_dist_params fails as expected extracting without sample size", {
   expect_message(
     params <- .calc_dist_params(
-      prob_dist = "gamma",
-      prob_dist_params = NA,
+      prob_distribution = "gamma",
       summary_stats = create_summary_stats(
         median = 10,
         lower_range = 5,
@@ -88,8 +82,7 @@ test_that(".calc_dist_params fails as expected extracting without sample size",
 test_that(".calc_dist_params messages as expected without summary stats", {
   expect_message(
     params <- .calc_dist_params(
-      prob_dist = "gamma",
-      prob_dist_params = NA,
+      prob_distribution = "gamma",
       summary_stats = create_summary_stats(mean = 5, median = 5),
       sample_size = NA
     ),
diff --git a/tests/testthat/test-checkers.R b/tests/testthat/test-checkers.R
index c8d19117e..8b0329c80 100644
--- a/tests/testthat/test-checkers.R
+++ b/tests/testthat/test-checkers.R
@@ -3,8 +3,10 @@ test_that("is_parameterised works as expected with epiparameter parameters", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation_period",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1)
+    )
   ))
   expect_true(is_parameterised(ep))
 })
diff --git a/tests/testthat/test-coercion.R b/tests/testthat/test-coercion.R
index d2405cd7f..354968cff 100644
--- a/tests/testthat/test-coercion.R
+++ b/tests/testthat/test-coercion.R
@@ -27,7 +27,7 @@ test_that("as_epiparameter works for ebola infectious period (issue #327 & #306)
   expect_s3_class(ebola_infectiousness_epiparameter, class = "epiparameter")
   # Lau 2017 has information to parameterise an exponential distribution
   expect_s3_class(
-    ebola_infectiousness_epiparameter$prob_dist,
+    ebola_infectiousness_epiparameter$prob_distribution,
     class = "distribution"
   )
   # populate mean summary statistics
@@ -52,7 +52,7 @@ test_that("as_epiparameter works for lassa incubation period (issue #306)", {
   expect_s3_class(lassa_incub_epiparameter, class = "epiparameter")
   # Akhmetzhanov 2019 has information to parameterise a gamma distribution
   expect_s3_class(
-    lassa_incub_epiparameter$prob_dist,
+    lassa_incub_epiparameter$prob_distribution,
     class = "distribution"
   )
   # populate mean and sd summary statistics without uncertainty
@@ -77,7 +77,7 @@ test_that("as_epiparameter works for ebola serial interval (issue #303)", {
   expect_s3_class(ebola_serial_epiparameter, class = "epiparameter")
   # Chan 2020 has information to parameterise a gamma distribution
   expect_s3_class(
-    ebola_serial_epiparameter$prob_dist,
+    ebola_serial_epiparameter$prob_distribution,
     class = "distribution"
   )
   # populate mean and sd summary statistics without uncertainty
@@ -103,7 +103,7 @@ test_that("as_epiparameter works for ebola SI assumed prob_dist (issue #310)", {
   expect_s3_class(ebola_serial_epiparameter, class = "epiparameter")
   # Faye 2015 parameterise an assumed gamma distribution
   expect_s3_class(
-    ebola_serial_epiparameter$prob_dist,
+    ebola_serial_epiparameter$prob_distribution,
     class = "distribution"
   )
   # populate mean and sd summary statistics without uncertainty
@@ -128,7 +128,7 @@ test_that("as_epiparameter works for lassa incubation overwritten prob_dist", {
   expect_s3_class(lassa_incub_epiparameter, class = "epiparameter")
   # Akhmetzhanov 2019 has information to parameterise a gamma distribution
   expect_s3_class(
-    lassa_incub_epiparameter$prob_dist,
+    lassa_incub_epiparameter$prob_distribution,
     class = "distribution"
   )
   expect_identical(family(lassa_incub_epiparameter), "lnorm")
@@ -176,13 +176,12 @@ test_that("as_epiparameter fails as expected with overwritten prob_dist", {
         ebola_params$article_label == "Fallah 2015 (1)"
     ),
   ]
-  # suppress warning and message about citation
-  ebola_si_lnorm <- suppressWarnings(
-    suppressMessages(
-      as_epiparameter(ebola_si, prob_dist = "lnorm")
-    )
+  expect_error(
+    suppressWarnings(
+      suppressMessages(
+         as_epiparameter(ebola_si, prob_dist = "lnorm")
+      )
+    ),
+    regexp = "Incorrect parameters provided for probability distribution."
   )
-  expect_s3_class(ebola_si_lnorm, class = "epiparameter")
-  expect_false(is_parameterised(ebola_si_lnorm))
-  expect_identical(family(ebola_si_lnorm), "lnorm")
 })
diff --git a/tests/testthat/test-convert_params.R b/tests/testthat/test-convert_params.R
index f968dfad4..02e673da4 100644
--- a/tests/testthat/test-convert_params.R
+++ b/tests/testthat/test-convert_params.R
@@ -2,8 +2,10 @@ ep <- suppressMessages(
   epiparameter(
     disease = "Ebola",
     epi_dist = "incubation period",
-    prob_distribution = "lnorm",
-    prob_distribution_params = c(meanlog = 1, sdlog = 1),
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "lnorm",
+      prob_distribution_params = c(meanlog = 1, sdlog = 1)
+    ),
     summary_stats = create_summary_stats(mean = 1, sd = 1)
   )
 )
@@ -196,7 +198,7 @@ test_that("convert_params_to_summary_stats.character fails as expected", {
 })
 
 test_that("convert_params_to_summary_stats.epiparameter fails as expected", {
-  ep$prob_dist <- "lnorm"
+  ep$prob_distribution <- "lnorm"
   expect_error(
     convert_params_to_summary_stats(ep),
     regexp = "<epiparameter> supplied has no parameters and none are suppled"
diff --git a/tests/testthat/test-create_epiparameter_prob_dist.R b/tests/testthat/test-create_epiparameter_prob_dist.R
deleted file mode 100644
index d737755c7..000000000
--- a/tests/testthat/test-create_epiparameter_prob_dist.R
+++ /dev/null
@@ -1,243 +0,0 @@
-test_that("create_prob_dist works for continuous gamma", {
-  res <- create_prob_dist(
-    prob_dist = "gamma",
-    prob_dist_params = c(shape = 1, scale = 1),
-    discretise = FALSE,
-    truncation = NA
-  )
-
-  expect_s3_class(res, "distribution")
-  expect_identical(family(res), "gamma")
-  expect_identical(
-    distributional::parameters(res),
-    data.frame(shape = 1, rate = 1)
-  )
-})
-
-test_that("create_prob_dist works for continuous lognormal", {
-  res <- create_prob_dist(
-    prob_dist = "lnorm",
-    prob_dist_params = c(meanlog = 1, sdlog = 1),
-    discretise = FALSE,
-    truncation = NA
-  )
-
-  expect_s3_class(res, "distribution")
-  expect_identical(family(res), "lognormal")
-  expect_identical(
-    distributional::parameters(res),
-    data.frame(mu = 1, sigma = 1)
-  )
-})
-
-test_that("create_prob_dist works for continuous Weibull", {
-  res <- create_prob_dist(
-    prob_dist = "weibull",
-    prob_dist_params = c(shape = 1, scale = 1),
-    discretise = FALSE,
-    truncation = NA
-  )
-
-  expect_s3_class(res, "distribution")
-  expect_identical(family(res), "weibull")
-  expect_identical(
-    distributional::parameters(res),
-    data.frame(shape = 1, scale = 1)
-  )
-})
-
-test_that("create_prob_dist works for negative binomial", {
-  res <- create_prob_dist(
-    prob_dist = "nbinom",
-    prob_dist_params = c(mean = 1, dispersion = 1),
-    discretise = FALSE,
-    truncation = NA
-  )
-
-  expect_s3_class(res, "distribution")
-  expect_identical(family(res), "negbin")
-  expect_identical(
-    distributional::parameters(res),
-    data.frame(n = 1, p = 0.5)
-  )
-})
-
-test_that("create_prob_dist works for geometric", {
-  res <- create_prob_dist(
-    prob_dist = "geom",
-    prob_dist_params = c(mean = 1),
-    discretise = FALSE,
-    truncation = NA
-  )
-
-  expect_s3_class(res, "distribution")
-  expect_identical(family(res), "geometric")
-  expect_identical(
-    distributional::parameters(res),
-    data.frame(p = 1)
-  )
-})
-
-test_that("create_prob_dist works for poisson", {
-  res <- create_prob_dist(
-    prob_dist = "pois",
-    prob_dist_params = c(mean = 1),
-    discretise = FALSE,
-    truncation = NA
-  )
-
-  expect_s3_class(res, "distribution")
-  expect_identical(family(res), "poisson")
-  expect_identical(
-    distributional::parameters(res),
-    data.frame(l = 1)
-  )
-})
-
-test_that("create_prob_dist works for exponential", {
-  res <- create_prob_dist(
-    prob_dist = "exp",
-    prob_dist_params = c(rate = 2),
-    discretise = FALSE,
-    truncation = NA
-  )
-
-  expect_s3_class(res, "distribution")
-  expect_identical(family(res), "exponential")
-  expect_identical(
-    distributional::parameters(res),
-    data.frame(rate = 2)
-  )
-})
-
-test_that("create_prob_dist works for discrete gamma", {
-  res <- create_prob_dist(
-    prob_dist = "gamma",
-    prob_dist_params = c(shape = 1, scale = 1),
-    discretise = TRUE,
-    truncation = NA
-  )
-
-  expect_s3_class(res, "distcrete")
-  expect_identical(res$name, "gamma")
-  expect_identical(
-    res$parameters,
-    list(shape = 1, scale = 1)
-  )
-})
-
-test_that("create_prob_dist works for discrete lognormal", {
-  res <- create_prob_dist(
-    prob_dist = "lnorm",
-    prob_dist_params = c(meanlog = 1, sdlog = 1),
-    discretise = TRUE,
-    truncation = NA
-  )
-
-  expect_s3_class(res, "distcrete")
-  expect_identical(res$name, "lnorm")
-  expect_identical(
-    res$parameters,
-    list(meanlog = 1, sdlog = 1)
-  )
-})
-
-test_that("create_prob_dist works for discrete Weibull", {
-  res <- create_prob_dist(
-    prob_dist = "weibull",
-    prob_dist_params = c(shape = 1, scale = 1),
-    discretise = TRUE,
-    truncation = NA
-  )
-
-  expect_s3_class(res, "distcrete")
-  expect_identical(res$name, "weibull")
-  expect_identical(
-    res$parameters,
-    list(shape = 1, scale = 1)
-  )
-})
-
-test_that("create_prob_dist works for discrete normal", {
-  res <- create_prob_dist(
-    prob_dist = "norm",
-    prob_dist_params = c(mean = 1, sd = 1),
-    discretise = TRUE,
-    truncation = NA
-  )
-
-  expect_s3_class(res, "distcrete")
-  expect_identical(res$name, "norm")
-  expect_identical(
-    res$parameters,
-    list(mean = 1, sd = 1)
-  )
-})
-
-test_that("create_prob_dist works for truncated continuous", {
-  res <- create_prob_dist(
-    prob_dist = "gamma",
-    prob_dist_params = c(shape = 1, scale = 1),
-    discretise = FALSE,
-    truncation = 10
-  )
-
-  expect_s3_class(res, "distribution")
-  expect_identical(family(res), "truncated")
-})
-
-test_that("create_prob_dist works passing args via ...", {
-  dist1 <- create_prob_dist(
-    prob_dist = "gamma",
-    prob_dist_params = c(shape = 1, scale = 1),
-    discretise = TRUE,
-    truncation = NA,
-    w = 1
-  )
-  dist2 <- create_prob_dist(
-    prob_dist = "gamma",
-    prob_dist_params = c(shape = 1, scale = 1),
-    discretise = TRUE,
-    truncation = NA,
-    w = 0.5
-  )
-  expect_false(identical(dist1, dist2))
-  expect_identical(dist1$w, 1)
-  expect_identical(dist2$w, 0.5)
-})
-
-test_that("create_prob_dist fails for unrecognised discretised", {
-  expect_error(
-    create_prob_dist(
-      prob_dist = "distribution",
-      prob_dist_params = c(shape = 1, scale = 1),
-      discretise = TRUE,
-      truncation = NA
-    ),
-    regexp = "(arg)*(should be one of)*(gamma)*(lnorm)*(weibull)"
-  )
-})
-
-test_that("create_prob_dist fails for unrecognised non-discretised", {
-  expect_error(
-    create_prob_dist(
-      prob_dist = "distribution",
-      prob_dist_params = c(shape = 1, scale = 1),
-      discretise = FALSE,
-      truncation = NA
-    ),
-    regexp = "Did not recognise distribution name"
-  )
-})
-
-test_that("create_prob_dist errors for discrete truncation", {
-  expect_error(
-    create_prob_dist(
-      prob_dist = "gamma",
-      prob_dist_params = c(shape = 1, scale = 1),
-      discretise = TRUE,
-      truncation = 10
-    ),
-    regexp = "Truncation is not yet implemented for discrete distributions"
-  )
-})
diff --git a/tests/testthat/test-create_prob_distribution.R b/tests/testthat/test-create_prob_distribution.R
new file mode 100644
index 000000000..f7e4eb980
--- /dev/null
+++ b/tests/testthat/test-create_prob_distribution.R
@@ -0,0 +1,243 @@
+test_that("create_prob_distribution works for continuous gamma", {
+  res <- create_prob_distribution(
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 1, scale = 1),
+    discretise = FALSE,
+    truncation = NA
+  )
+
+  expect_s3_class(res, "distribution")
+  expect_identical(family(res), "gamma")
+  expect_identical(
+    distributional::parameters(res),
+    data.frame(shape = 1, rate = 1)
+  )
+})
+
+test_that("create_prob_distribution works for continuous lognormal", {
+  res <- create_prob_distribution(
+    prob_distribution = "lnorm",
+    prob_distribution_params = c(meanlog = 1, sdlog = 1),
+    discretise = FALSE,
+    truncation = NA
+  )
+
+  expect_s3_class(res, "distribution")
+  expect_identical(family(res), "lognormal")
+  expect_identical(
+    distributional::parameters(res),
+    data.frame(mu = 1, sigma = 1)
+  )
+})
+
+test_that("create_prob_distribution works for continuous Weibull", {
+  res <- create_prob_distribution(
+    prob_distribution = "weibull",
+    prob_distribution_params = c(shape = 1, scale = 1),
+    discretise = FALSE,
+    truncation = NA
+  )
+
+  expect_s3_class(res, "distribution")
+  expect_identical(family(res), "weibull")
+  expect_identical(
+    distributional::parameters(res),
+    data.frame(shape = 1, scale = 1)
+  )
+})
+
+test_that("create_prob_distribution works for negative binomial", {
+  res <- create_prob_distribution(
+    prob_distribution = "nbinom",
+    prob_distribution_params = c(mean = 1, dispersion = 1),
+    discretise = FALSE,
+    truncation = NA
+  )
+
+  expect_s3_class(res, "distribution")
+  expect_identical(family(res), "negbin")
+  expect_identical(
+    distributional::parameters(res),
+    data.frame(n = 1, p = 0.5)
+  )
+})
+
+test_that("create_prob_distribution works for geometric", {
+  res <- create_prob_distribution(
+    prob_distribution = "geom",
+    prob_distribution_params = c(mean = 1),
+    discretise = FALSE,
+    truncation = NA
+  )
+
+  expect_s3_class(res, "distribution")
+  expect_identical(family(res), "geometric")
+  expect_identical(
+    distributional::parameters(res),
+    data.frame(p = 1)
+  )
+})
+
+test_that("create_prob_distribution works for poisson", {
+  res <- create_prob_distribution(
+    prob_distribution = "pois",
+    prob_distribution_params = c(mean = 1),
+    discretise = FALSE,
+    truncation = NA
+  )
+
+  expect_s3_class(res, "distribution")
+  expect_identical(family(res), "poisson")
+  expect_identical(
+    distributional::parameters(res),
+    data.frame(l = 1)
+  )
+})
+
+test_that("create_prob_distribution works for exponential", {
+  res <- create_prob_distribution(
+    prob_distribution = "exp",
+    prob_distribution_params = c(rate = 2),
+    discretise = FALSE,
+    truncation = NA
+  )
+
+  expect_s3_class(res, "distribution")
+  expect_identical(family(res), "exponential")
+  expect_identical(
+    distributional::parameters(res),
+    data.frame(rate = 2)
+  )
+})
+
+test_that("create_prob_distribution works for discrete gamma", {
+  res <- create_prob_distribution(
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 1, scale = 1),
+    discretise = TRUE,
+    truncation = NA
+  )
+
+  expect_s3_class(res, "distcrete")
+  expect_identical(res$name, "gamma")
+  expect_identical(
+    res$parameters,
+    list(shape = 1, scale = 1)
+  )
+})
+
+test_that("create_prob_distribution works for discrete lognormal", {
+  res <- create_prob_distribution(
+    prob_distribution = "lnorm",
+    prob_distribution_params = c(meanlog = 1, sdlog = 1),
+    discretise = TRUE,
+    truncation = NA
+  )
+
+  expect_s3_class(res, "distcrete")
+  expect_identical(res$name, "lnorm")
+  expect_identical(
+    res$parameters,
+    list(meanlog = 1, sdlog = 1)
+  )
+})
+
+test_that("create_prob_distribution works for discrete Weibull", {
+  res <- create_prob_distribution(
+    prob_distribution = "weibull",
+    prob_distribution_params = c(shape = 1, scale = 1),
+    discretise = TRUE,
+    truncation = NA
+  )
+
+  expect_s3_class(res, "distcrete")
+  expect_identical(res$name, "weibull")
+  expect_identical(
+    res$parameters,
+    list(shape = 1, scale = 1)
+  )
+})
+
+test_that("create_prob_distribution works for discrete normal", {
+  res <- create_prob_distribution(
+    prob_distribution = "norm",
+    prob_distribution_params = c(mean = 1, sd = 1),
+    discretise = TRUE,
+    truncation = NA
+  )
+
+  expect_s3_class(res, "distcrete")
+  expect_identical(res$name, "norm")
+  expect_identical(
+    res$parameters,
+    list(mean = 1, sd = 1)
+  )
+})
+
+test_that("create_prob_distribution works for truncated continuous", {
+  res <- create_prob_distribution(
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 1, scale = 1),
+    discretise = FALSE,
+    truncation = 10
+  )
+
+  expect_s3_class(res, "distribution")
+  expect_identical(family(res), "truncated")
+})
+
+test_that("create_prob_distribution works passing args via ...", {
+  dist1 <- create_prob_distribution(
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 1, scale = 1),
+    discretise = TRUE,
+    truncation = NA,
+    w = 1
+  )
+  dist2 <- create_prob_distribution(
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 1, scale = 1),
+    discretise = TRUE,
+    truncation = NA,
+    w = 0.5
+  )
+  expect_false(identical(dist1, dist2))
+  expect_identical(dist1$w, 1)
+  expect_identical(dist2$w, 0.5)
+})
+
+test_that("create_prob_distribution fails for unrecognised discretised", {
+  expect_error(
+    create_prob_distribution(
+      prob_distribution = "distribution",
+      prob_distribution_params = c(shape = 1, scale = 1),
+      discretise = TRUE,
+      truncation = NA
+    ),
+    regexp = "Incorrect parameters provided for probability distribution."
+  )
+})
+
+test_that("create_prob_distribution fails for unrecognised non-discretised", {
+  expect_error(
+    create_prob_distribution(
+      prob_distribution = "distribution",
+      prob_distribution_params = c(shape = 1, scale = 1),
+      discretise = FALSE,
+      truncation = NA
+    ),
+    regexp = "Incorrect parameters provided for probability distribution."
+  )
+})
+
+test_that("create_prob_distribution errors for discrete truncation", {
+  expect_error(
+    create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1),
+      discretise = TRUE,
+      truncation = 10
+    ),
+    regexp = "Truncation is not yet implemented for discrete distributions"
+  )
+})
diff --git a/tests/testthat/test-epiparameter-utils.R b/tests/testthat/test-epiparameter-utils.R
index c16973f73..1413ec2c7 100644
--- a/tests/testthat/test-epiparameter-utils.R
+++ b/tests/testthat/test-epiparameter-utils.R
@@ -91,14 +91,14 @@ test_that("create_citation works with PMID", {
 
 test_that(".clean_params works as expected for gamma", {
   params <- .clean_params(
-    prob_dist = "gamma",
-    prob_dist_params = c(shape = 1, scale = 1)
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 1, scale = 1)
   )
   expect_identical(params, c(shape = 1, scale = 1))
 
   params <- .clean_params(
-    prob_dist = "gamma",
-    prob_dist_params = c(shape = 1, rate = 0.5)
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 1, rate = 0.5)
   )
   expect_identical(params, c(shape = 1, scale = 2))
 })
@@ -106,8 +106,8 @@ test_that(".clean_params works as expected for gamma", {
 test_that(".clean_params fails when gamma parameters are incorrect", {
   expect_error(
     .clean_params(
-      prob_dist = "gamma",
-      prob_dist_params = c(meanlog = 1, sdlog = 1)
+      prob_distribution = "gamma",
+      prob_distribution_params = c(meanlog = 1, sdlog = 1)
     ),
     regexp = "Invalid parameterisation for gamma distribution"
   )
@@ -115,14 +115,14 @@ test_that(".clean_params fails when gamma parameters are incorrect", {
 
 test_that(".clean_params works as expected for lnorm", {
   params <- .clean_params(
-    prob_dist = "lnorm",
-    prob_dist_params = c(meanlog = 1, sdlog = 1)
+    prob_distribution = "lnorm",
+    prob_distribution_params = c(meanlog = 1, sdlog = 1)
   )
   expect_identical(params, c(meanlog = 1, sdlog = 1))
 
   params <- .clean_params(
-    prob_dist = "lnorm",
-    prob_dist_params = c(mu = 2, sigma = 2)
+    prob_distribution = "lnorm",
+    prob_distribution_params = c(mu = 2, sigma = 2)
   )
   expect_identical(params, c(meanlog = 2, sdlog = 2))
 })
@@ -130,8 +130,8 @@ test_that(".clean_params works as expected for lnorm", {
 test_that(".clean_params fails when lnorm parameters are incorrect", {
   expect_error(
     .clean_params(
-      prob_dist = "lnorm",
-      prob_dist_params = c(shape = 1, scale = 1)
+      prob_distribution = "lnorm",
+      prob_distribution_params = c(shape = 1, scale = 1)
     ),
     regexp = "Invalid parameterisation for lnorm distribution"
   )
@@ -139,8 +139,8 @@ test_that(".clean_params fails when lnorm parameters are incorrect", {
 
 test_that(".clean_params works as expected for weibull", {
   params <- .clean_params(
-    prob_dist = "weibull",
-    prob_dist_params = c(shape = 1, scale = 1)
+    prob_distribution = "weibull",
+    prob_distribution_params = c(shape = 1, scale = 1)
   )
   expect_identical(params, c(shape = 1, scale = 1))
 })
@@ -148,8 +148,8 @@ test_that(".clean_params works as expected for weibull", {
 test_that(".clean_params fails when weibull parameters are incorrect", {
   expect_error(
     .clean_params(
-      prob_dist = "weibull",
-      prob_dist_params = c(meanlog = 1, sdlog = 1)
+      prob_distribution = "weibull",
+      prob_distribution_params = c(meanlog = 1, sdlog = 1)
     ),
     regexp = "Invalid parameterisation for weibull distribution"
   )
@@ -157,14 +157,14 @@ test_that(".clean_params fails when weibull parameters are incorrect", {
 
 test_that(".clean_params works as expected for nbinom", {
   params <- .clean_params(
-    prob_dist = "nbinom",
-    prob_dist_params = c(n = 2, p = 0.5)
+    prob_distribution = "nbinom",
+    prob_distribution_params = c(n = 2, p = 0.5)
   )
   expect_identical(params, c(mean = 2, dispersion = 2))
 
   params <- .clean_params(
-    prob_dist = "nbinom",
-    prob_dist_params = c(mean = 1, dispersion = 1)
+    prob_distribution = "nbinom",
+    prob_distribution_params = c(mean = 1, dispersion = 1)
   )
   expect_identical(params, c(mean = 1, dispersion = 1))
 })
@@ -172,8 +172,8 @@ test_that(".clean_params works as expected for nbinom", {
 test_that(".clean_params fails when nbinom parameters are incorrect", {
   expect_error(
     .clean_params(
-      prob_dist = "nbinom",
-      prob_dist_params = c(meanlog = 1, sdlog = 1)
+      prob_distribution = "nbinom",
+      prob_distribution_params = c(meanlog = 1, sdlog = 1)
     ),
     regexp = "Invalid parameterisation for nbinom distribution"
   )
@@ -181,20 +181,20 @@ test_that(".clean_params fails when nbinom parameters are incorrect", {
 
 test_that(".clean_params works as expected for geom", {
   params <- .clean_params(
-    prob_dist = "geom",
-    prob_dist_params = c(prob = 0.5)
+    prob_distribution = "geom",
+    prob_distribution_params = c(prob = 0.5)
   )
   expect_identical(params, c(prob = 0.5))
 
   params <- .clean_params(
-    prob_dist = "geom",
-    prob_dist_params = c(p = 0.5)
+    prob_distribution = "geom",
+    prob_distribution_params = c(p = 0.5)
   )
   expect_identical(params, c(prob = 0.5))
 
   params <- .clean_params(
-    prob_dist = "geom",
-    prob_dist_params = c(mean = 2)
+    prob_distribution = "geom",
+    prob_distribution_params = c(mean = 2)
   )
   expect_identical(params, c(prob = 0.5))
 })
@@ -202,8 +202,8 @@ test_that(".clean_params works as expected for geom", {
 test_that(".clean_params fails when geom parameters are incorrect", {
   expect_error(
     .clean_params(
-      prob_dist = "geom",
-      prob_dist_params = c(meanlog = 1, sdlog = 1)
+      prob_distribution = "geom",
+      prob_distribution_params = c(meanlog = 1, sdlog = 1)
     ),
     regexp = "Invalid parameterisation for geom distribution"
   )
@@ -211,20 +211,20 @@ test_that(".clean_params fails when geom parameters are incorrect", {
 
 test_that(".clean_params works as expected for pois", {
   params <- .clean_params(
-    prob_dist = "pois",
-    prob_dist_params = c(mean = 0.5)
+    prob_distribution = "pois",
+    prob_distribution_params = c(mean = 0.5)
   )
   expect_identical(params, c(mean = 0.5))
 
   params <- .clean_params(
-    prob_dist = "pois",
-    prob_dist_params = c(l = 0.5)
+    prob_distribution = "pois",
+    prob_distribution_params = c(l = 0.5)
   )
   expect_identical(params, c(mean = 0.5))
 
   params <- .clean_params(
-    prob_dist = "pois",
-    prob_dist_params = c(lambda = 0.5)
+    prob_distribution = "pois",
+    prob_distribution_params = c(lambda = 0.5)
   )
   expect_identical(params, c(mean = 0.5))
 })
@@ -232,8 +232,8 @@ test_that(".clean_params works as expected for pois", {
 test_that(".clean_params fails when pois parameters are incorrect", {
   expect_error(
     .clean_params(
-      prob_dist = "pois",
-      prob_dist_params = c(means = 1)
+      prob_distribution = "pois",
+      prob_distribution_params = c(means = 1)
     ),
     regexp = "Invalid parameterisation for pois distribution"
   )
@@ -241,20 +241,20 @@ test_that(".clean_params fails when pois parameters are incorrect", {
 
 test_that(".clean_params works as expected for exp", {
   params <- .clean_params(
-    prob_dist = "exp",
-    prob_dist_params = c(rate = 2)
+    prob_distribution = "exp",
+    prob_distribution_params = c(rate = 2)
   )
   expect_identical(params, c(rate = 2))
 
   params <- .clean_params(
-    prob_dist = "exp",
-    prob_dist_params = c(lambda = 2)
+    prob_distribution = "exp",
+    prob_distribution_params = c(lambda = 2)
   )
   expect_identical(params, c(rate = 2))
 
   params <- .clean_params(
-    prob_dist = "exp",
-    prob_dist_params = c(mean = 0.5)
+    prob_distribution = "exp",
+    prob_distribution_params = c(mean = 0.5)
   )
   expect_identical(params, c(rate = 2))
 })
@@ -262,8 +262,8 @@ test_that(".clean_params works as expected for exp", {
 test_that(".clean_params fails when exp parameters are incorrect", {
   expect_error(
     .clean_params(
-      prob_dist = "exp",
-      prob_dist_params = c(means = 1)
+      prob_distribution = "exp",
+      prob_distribution_params = c(means = 1)
     ),
     regexp = "Invalid parameterisation for exp distribution"
   )
@@ -272,8 +272,8 @@ test_that(".clean_params fails when exp parameters are incorrect", {
 test_that(".clean_params fails as expected", {
   expect_error(
     .clean_params(
-      prob_dist = "distribution",
-      prob_dist_params = c(meanlog = 1, sdlog = 1)
+      prob_distribution = "distribution",
+      prob_distribution_params = c(meanlog = 1, sdlog = 1)
     ),
     regexp = "Invalid parameterisation for distribution distribution"
   )
diff --git a/tests/testthat/test-epiparameter.R b/tests/testthat/test-epiparameter.R
index ac5d91aff..b374d9e1f 100644
--- a/tests/testthat/test-epiparameter.R
+++ b/tests/testthat/test-epiparameter.R
@@ -3,8 +3,10 @@ test_that("epiparameter works with minimal viable input", {
   ebola_dist <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1)
+    )
   ))
 
   expect_s3_class(ebola_dist, class = "epiparameter")
@@ -18,8 +20,10 @@ test_that("epiparameter works with all arguments set", {
       disease = "MERS",
       pathogen = "MERS_CoV",
       epi_dist = "serial_interval",
-      prob_distribution = "lnorm",
-      prob_distribution_params = c(meanlog = 2, sdlog = 1),
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "lnorm",
+        prob_distribution_params = c(meanlog = 2, sdlog = 1)
+      ),
       uncertainty = list(
         meanlog = create_uncertainty(
           ci_limits = c(1, 3),
@@ -84,15 +88,20 @@ test_that("epiparameter works with default helper functions", {
     disease = "SARS",
     pathogen = "SARS_CoV",
     epi_dist = "onset_to_death",
-    prob_distribution = "lnorm",
-    prob_distribution_params = c(meanlog = 2, sdlog = 1),
-    uncertainty = create_uncertainty(),
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "lnorm",
+      prob_distribution_params = c(meanlog = 2, sdlog = 1),
+      discretise = FALSE,
+      truncation = NA
+    ),
+    uncertainty = list(
+      meanlog = create_uncertainty(),
+      sdlog = create_uncertainty()
+    ),
     summary_stats = create_summary_stats(),
     citation = create_citation(),
     metadata = create_metadata(),
     method_assess = create_method_assess(),
-    discretise = FALSE,
-    truncation = NA,
     notes = "No notes"
   ))
 
@@ -105,8 +114,10 @@ test_that("epiparameter fails as expected", {
     epiparameter(
       disease = 1,
       epi_dist = "incubation",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1)
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1)
+      )
     ),
     regexp = paste0(
       "Assertion on 'disease' failed: Must be of type ",
@@ -118,8 +129,10 @@ test_that("epiparameter fails as expected", {
     epiparameter(
       disease = "ebola",
       epi_dist = 1,
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1)
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1)
+      )
     ),
     regexp = paste0(
       "Assertion on 'epi_dist' failed: Must be of type ",
@@ -131,12 +144,10 @@ test_that("epiparameter fails as expected", {
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation",
-      prob_distribution = 1,
-      prob_distribution_params = c(shape = 1, scale = 1)
+      prob_distribution = 1
     ),
     regexp = paste0(
-      "(Assertion on 'prob_distribution' failed)*(Must be of type)*",
-      "(character)*(NULL)*(double)"
+      "epiparameter must contain a <distribution> or <distcrete> or NA"
     )
   )
 
@@ -161,8 +172,10 @@ test_that("epiparameter.print works as expected", {
   expect_snapshot(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1)
+    )
   ))
 })
 
@@ -170,9 +183,11 @@ test_that("epiparameter.print works as expected", {
   expect_snapshot(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1),
-    discretise = TRUE
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1),
+      discretise = TRUE
+    )
   ))
 })
 
@@ -180,8 +195,10 @@ test_that("epiparameter.plot does not produce an error", {
   ebola_dist <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1)
+    )
   ))
 
 
@@ -198,8 +215,10 @@ test_that("epiparameter.plot works with non-default x-axis", {
   ebola_dist <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1)
+    )
   ))
 
   expect_silent(
@@ -225,10 +244,12 @@ test_that("new_epiparameter works with minimal viable input", {
   epiparameter_obj <- suppressMessages(
     new_epiparameter(
       disease = "ebola",
-      pathogen = "ebola_virus",
+      pathogen = "ebola virus",
       epi_dist = "incubation",
-      prob_dist = "gamma",
-      prob_dist_params = c(shape = 1, scale = 1),
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1)
+      ),
       uncertainty = list(
         shape = create_uncertainty(
           ci_limits = c(0, 2),
@@ -247,8 +268,6 @@ test_that("new_epiparameter works with minimal viable input", {
         title = "Ebola incubation",
         journal = "Journal of Epi"
       ),
-      discretise = FALSE,
-      truncation = NA,
       notes = "No notes"
     )
   )
@@ -264,8 +283,12 @@ test_that("assert_epiparameter passes when expected", {
       disease = "ebola",
       pathogen = "ebola_virus",
       epi_dist = "incubation",
-      prob_dist = "gamma",
-      prob_dist_params = c(shape = 1, scale = 1),
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distributions_params = c(shape = 1, scale = 1),
+        discretise = FALSE,
+        truncation = NA
+      ),
       uncertainty = list(
         shape = create_uncertainty(
           ci_limits = c(0, 2),
@@ -285,9 +308,8 @@ test_that("assert_epiparameter passes when expected", {
         journal = "Journal of Epi",
         doi = "10.1872372hc"
       ),
-      discretise = FALSE,
-      truncation = NA,
-      notes = "No notes"
+      notes = "No notes",
+      auto_calc_params = FALSE
     )
   )
 
@@ -299,8 +321,12 @@ test_that("assert_epiparameter catches class faults when expected", {
     disease = "ebola",
     pathogen = "ebola_virus",
     epi_dist = "incubation",
-    prob_dist = "gamma",
-    prob_dist_params = c(shape = 1, scale = 1),
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_dist_params = c(shape = 1, scale = 1),
+      discretise = FALSE,
+      truncation = NA
+    ),
     uncertainty = list(
       shape = create_uncertainty(
         ci_limits = c(0, 2),
@@ -314,8 +340,7 @@ test_that("assert_epiparameter catches class faults when expected", {
       )
     ),
     citation = "Smith (2002) <10.128372837>",
-    discretise = FALSE,
-    truncation = NA
+    auto_calc_params = FALSE
   )
 
   epiparameter_obj$disease <- NULL
@@ -329,8 +354,12 @@ test_that("assert_epiparameter catches class faults when expected", {
     disease = "ebola",
     pathogen = "ebola_virus",
     epi_dist = "incubation",
-    prob_dist = "gamma",
-    prob_dist_params = c(shape = 1, scale = 1),
+    prob_distribution = create_prob_distribution(
+      prob_distribution =  "gamma",
+      prob_dist_params = c(shape = 1, scale = 1),
+      discretise = FALSE,
+      truncation = NA
+    ),
     uncertainty = list(
       shape = create_uncertainty(
         ci_limits = c(0, 2),
@@ -344,8 +373,7 @@ test_that("assert_epiparameter catches class faults when expected", {
       )
     ),
     citation = "Smith (2002) <10.128372837>",
-    discretise = FALSE,
-    truncation = NA
+    auto_calc_params = FALSE
   )
 
   epiparameter_obj$disease <- factor("disease")
@@ -359,8 +387,12 @@ test_that("assert_epiparameter catches class faults when expected", {
     disease = "ebola",
     pathogen = "ebola_virus",
     epi_dist = "incubation",
-    prob_dist = "gamma",
-    prob_dist_params = c(shape = 1, scale = 1),
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_dist_params = c(shape = 1, scale = 1),
+      discretise = FALSE,
+      truncation = NA
+    ),
     uncertainty = list(
       shape = create_uncertainty(
         ci_limits = c(0, 2),
@@ -374,8 +406,7 @@ test_that("assert_epiparameter catches class faults when expected", {
       )
     ),
     citation = "Smith (2002) <10.128372837>",
-    discretise = FALSE,
-    truncation = NA
+    auto_calc_params = FALSE
   )
 
   epiparameter_obj$epi_dist <- c("incubation", "period")
@@ -421,8 +452,10 @@ test_that("density works as expected on continuous epiparameter object", {
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1)
+      prob_distribution = create_prob_distribution(
+        prob_distribution =  "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1)
+      )
     )
   )
 
@@ -452,9 +485,11 @@ test_that("density works as expected on discrete epiparameter object", {
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1),
-      discretise = TRUE
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1),
+        discretise = TRUE
+      )
     )
   )
 
@@ -485,8 +520,10 @@ test_that("density works as expected on continuous epiparameter object with vect
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1)
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1)
+      )
     )
   )
 
@@ -502,9 +539,11 @@ test_that("density works as expected on discrete epiparameter object with vector
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1),
-      discretise = TRUE
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1),
+        discretise = TRUE
+      )
     )
   )
 
@@ -519,8 +558,10 @@ test_that("cdf works as expected on continuous epiparameter object", {
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1)
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1)
+      )
     )
   )
 
@@ -554,9 +595,11 @@ test_that("cdf works as expected on discrete epiparameter object", {
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1),
-      discretise = TRUE
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1),
+        discretise = TRUE
+      )
     )
   )
 
@@ -591,8 +634,10 @@ test_that("cdf works as expected on continuous epiparameter object with vector
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1)
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1)
+      )
     )
   )
 
@@ -609,9 +654,11 @@ test_that("cdf works as expected on discrete epiparameter object with vector
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1),
-      discretise = TRUE
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1),
+        discretise = TRUE
+      )
     )
   )
 
@@ -627,8 +674,10 @@ test_that("quantile works as expected on continuous epiparameter object", {
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1)
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1)
+      )
     )
   )
 
@@ -653,9 +702,11 @@ test_that("quantile works as expected on discrete epiparameter object", {
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1),
-      discretise = TRUE
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1),
+        discretise = TRUE
+      )
     )
   )
 
@@ -681,8 +732,10 @@ test_that("quantile works as expected on continuous epiparameter object with vec
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1)
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1)
+      )
     )
   )
 
@@ -698,9 +751,11 @@ test_that("quantile works as expected on discrete epiparameter object with vecto
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1),
-      discretise = TRUE
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1),
+        discretise = TRUE
+      )
     )
   )
 
@@ -715,8 +770,10 @@ test_that("generate works as expected on continuous epiparameter object", {
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1)
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1)
+      )
     )
   )
 
@@ -736,9 +793,11 @@ test_that("generate works as expected on discrete epiparameter object", {
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1),
-      discretise = TRUE
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1),
+        discretise = TRUE
+      )
     )
   )
 
@@ -759,8 +818,10 @@ test_that("generate fails as expected on continuous epiparameter object with vec
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1)
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1)
+      )
     )
   )
 
@@ -773,9 +834,11 @@ test_that("generate fails as expected on discrete epiparameter object with vecto
     epiparameter(
       disease = "ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1),
-      discretise = TRUE
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1),
+        discretise = TRUE
+      )
     )
   )
 
@@ -787,8 +850,10 @@ test_that("is_epiparameter returns TRUE when expected", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "serial_interval",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1)
+    )
   ))
 
   expect_true(is_epiparameter(ep))
@@ -798,8 +863,10 @@ test_that("is_epiparameter returns FALSE when expected", {
   false_ep <- list(
     disease = "ebola",
     epi_dist = "serial_interval",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1)
+    )
   )
 
   expect_false(is_epiparameter(false_ep))
@@ -810,14 +877,16 @@ test_that("discretise works as expected on continuous gamma", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1)
+    )
   ))
   ep <- discretise(ep)
 
-  expect_s3_class(ep$prob_dist, "distcrete")
-  expect_identical(ep$prob_dist$parameters, list(shape = 1, scale = 1))
-  expect_identical(ep$prob_dist$name, "gamma")
+  expect_s3_class(ep$prob_distribution, "distcrete")
+  expect_identical(ep$prob_distribution$parameters, list(shape = 1, scale = 1))
+  expect_identical(ep$prob_distribution$name, "gamma")
 })
 
 test_that("discretise works as expected on continuous lognormal", {
@@ -825,14 +894,16 @@ test_that("discretise works as expected on continuous lognormal", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "lnorm",
-    prob_distribution_params = c(meanlog = 1, sdlog = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "lnorm",
+      prob_distribution_params = c(meanlog = 1, sdlog = 1)
+    )
   ))
   ep <- discretise(ep)
 
-  expect_s3_class(ep$prob_dist, "distcrete")
-  expect_identical(ep$prob_dist$parameters, list(meanlog = 1, sdlog = 1))
-  expect_identical(ep$prob_dist$name, "lnorm")
+  expect_s3_class(ep$prob_distribution, "distcrete")
+  expect_identical(ep$prob_distribution$parameters, list(meanlog = 1, sdlog = 1))
+  expect_identical(ep$prob_distribution$name, "lnorm")
 })
 
 test_that("discretise works as expected on discretised dist", {
@@ -840,18 +911,20 @@ test_that("discretise works as expected on discretised dist", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1),
-    discretise = TRUE
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1),
+      discretise = TRUE
+    )
   ))
   expect_message(
     discretise(ep),
     regexp = "Distribution in `epiparameter` is already discretised"
   )
 
-  expect_s3_class(ep$prob_dist, "distcrete")
-  expect_identical(ep$prob_dist$parameters, list(shape = 1, scale = 1))
-  expect_identical(ep$prob_dist$name, "gamma")
+  expect_s3_class(ep$prob_distribution, "distcrete")
+  expect_identical(ep$prob_distribution$parameters, list(shape = 1, scale = 1))
+  expect_identical(ep$prob_distribution$name, "gamma")
 })
 
 test_that("discretise works as expected on truncated dist", {
@@ -859,9 +932,11 @@ test_that("discretise works as expected on truncated dist", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1),
-    truncation = 10
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1),
+      truncation = 10
+    )
   ))
 
   expect_warning(
@@ -874,9 +949,9 @@ test_that("discretise works as expected on truncated dist", {
 
   ep <- suppressWarnings(discretise(ep))
 
-  expect_s3_class(ep$prob_dist, "distcrete")
-  expect_identical(ep$prob_dist$parameters, list(shape = 1, scale = 1))
-  expect_identical(ep$prob_dist$name, "gamma")
+  expect_s3_class(ep$prob_distribution, "distcrete")
+  expect_identical(ep$prob_distribution$parameters, list(shape = 1, scale = 1))
+  expect_identical(ep$prob_distribution$name, "gamma")
 })
 
 test_that("discretise fails as expected on non-epiparameter object", {
@@ -896,8 +971,10 @@ test_that("parameters works as expected on continuous gamma", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1)
+    )
   ))
   params <- get_parameters(ep)
 
@@ -910,8 +987,10 @@ test_that("parameters works as expected on continuous lognormal", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "lnorm",
-    prob_distribution_params = c(meanlog = 1, sdlog = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "lnorm",
+      prob_distribution_params = c(meanlog = 1, sdlog = 1)
+    )
   ))
   params <- get_parameters(ep)
 
@@ -924,9 +1003,11 @@ test_that("parameters works as expected on discretised dist", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1),
-    discretise = TRUE
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1),
+      discretise = TRUE
+    )
   ))
   params <- get_parameters(ep)
 
@@ -939,9 +1020,11 @@ test_that("parameters works as expected on truncated dist", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1),
-    truncation = 10
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1),
+      truncation = 10
+    )
   ))
   params <- get_parameters(ep)
 
@@ -972,8 +1055,10 @@ test_that("family works as expected for distributional", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation_period",
-    prob_distribution = "lnorm",
-    prob_distribution_params = c(meanlog = 1, sdlog = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "lnorm",
+      prob_distribution_params = c(meanlog = 1, sdlog = 1)
+    )
   ))
   expect_identical(family(ep), "lnorm")
 })
@@ -983,9 +1068,11 @@ test_that("family works as expected for distcrete", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation_period",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1),
-    discretise = TRUE
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1),
+      discretise = TRUE
+    )
   ))
   expect_identical(family(ep), "gamma")
 })
@@ -995,9 +1082,11 @@ test_that("family works as expected for distributional truncated", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation_period",
-    prob_distribution = "weibull",
-    prob_distribution_params = c(shape = 1, scale = 1),
-    truncation = 10
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "weibull",
+      prob_distribution_params = c(shape = 1, scale = 1),
+      truncation = 10
+    )
   ))
   expect_identical(family(ep), "weibull")
 })
@@ -1007,8 +1096,10 @@ test_that("is_truncated works as expected for continuous distributions", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation_period",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1)
+    )
   ))
   expect_false(is_truncated(ep))
 })
@@ -1018,9 +1109,11 @@ test_that("is_truncated works as expected for discretised distributions", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation_period",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1),
-    discretise = TRUE
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1),
+      discretise = TRUE
+    )
   ))
   expect_false(is_truncated(ep))
 })
@@ -1030,9 +1123,11 @@ test_that("is_truncated works as expected for truncated distributions", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation_period",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1),
-    truncation = 10
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1),
+      truncation = 10
+    )
   ))
   expect_true(is_truncated(ep))
 })
@@ -1042,16 +1137,20 @@ test_that("is_continuous works as expected for continuous distributions", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation period",
-    prob_distribution = "lnorm",
-    prob_distribution_params = c(meanlog = 1, sdlog = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "lnorm",
+      prob_distribution_params = c(meanlog = 1, sdlog = 1)
+    )
   ))
   expect_true(is_continuous(ep))
 
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation period",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1)
+    )
   ))
   expect_true(is_continuous(ep))
 })
@@ -1061,17 +1160,21 @@ test_that("is_continuous works as expected for discrete distributions", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "offspring distribution",
-    prob_distribution = "nbinom",
-    prob_distribution_params = c(mean = 2, dispersion = 0.5)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "nbinom",
+      prob_distribution_params = c(mean = 2, dispersion = 0.5)
+    )
   ))
   expect_false(is_continuous(ep))
 
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation period",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1),
-    discretise = TRUE
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1),
+      discretise = TRUE
+    )
   ))
   expect_false(is_continuous(ep))
 })
@@ -1093,8 +1196,10 @@ test_that("mean works as expected with params and no mean", {
     epiparameter(
       disease = "Ebola",
       epi_dist = "incubation_period",
-      prob_distribution = "gamma",
-      prob_distribution_params = c(shape = 1, scale = 1)
+      prob_distribution = create_prob_distribution(
+        prob_distribution = "gamma",
+        prob_distribution_params = c(shape = 1, scale = 1)
+      )
     )
   )
   expect_identical(mean(ep), 1)
@@ -1128,8 +1233,10 @@ test_that("as.function works as expected for density", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1)
+    )
   ))
   ep_func <- as.function(ep, func_type = "density")
   expect_type(ep_func, type = "closure")
@@ -1141,8 +1248,10 @@ test_that("as.function works as expected for cdf", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1)
+    )
   ))
   ep_func <- as.function(ep, func_type = "cdf")
   expect_type(ep_func, type = "closure")
@@ -1154,8 +1263,10 @@ test_that("as.function works as expected for generate", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1)
+    )
   ))
   ep_func <- as.function(ep, func_type = "generate")
   expect_type(ep_func, type = "closure")
@@ -1167,8 +1278,10 @@ test_that("as.function works as expected for quantile", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1)
+    )
   ))
   ep_func <- as.function(ep, func_type = "quantile")
   expect_type(ep_func, type = "closure")
@@ -1180,8 +1293,10 @@ test_that("as.function fails as expected", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1)
+    )
   ))
   expect_error(
     as.function(ep, func_type = "random"),
@@ -1205,8 +1320,10 @@ test_that("as.data.frame works for <epiparameter>", {
   ep <- suppressMessages(epiparameter(
     disease = "ebola",
     epi_dist = "incubation",
-    prob_distribution = "gamma",
-    prob_distribution_params = c(shape = 1, scale = 1)
+    prob_distribution = create_prob_distribution(
+      prob_distribution = "gamma",
+      prob_distribution_params = c(shape = 1, scale = 1)
+    )
   ))
   df <- as.data.frame(ep)
   expect_s3_class(df, class = "data.frame")
diff --git a/vignettes/epiparameter.Rmd b/vignettes/epiparameter.Rmd
index 5a4ced39b..d3f5d5e39 100644
--- a/vignettes/epiparameter.Rmd
+++ b/vignettes/epiparameter.Rmd
@@ -93,8 +93,10 @@ covid_incubation <- epiparameter(
   disease = "COVID-19",
   pathogen = "SARS-CoV-2",
   epi_dist = "incubation period",
-  prob_distribution = "gamma",
-  prob_distribution_params = c(shape = 2, scale = 1),
+  prob_distribution = create_prob_distribution(
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 2, scale = 1)
+  ),
   summary_stats = create_summary_stats(mean = 2),
   citation = create_citation(
     author = person(
diff --git a/vignettes/extract_convert.Rmd b/vignettes/extract_convert.Rmd
index cc5160b87..85e8e1089 100644
--- a/vignettes/extract_convert.Rmd
+++ b/vignettes/extract_convert.Rmd
@@ -82,8 +82,10 @@ ep <- epiparameter(
   disease = "<Disease name>",
   pathogen = "<Pathogen name>",
   epi_dist = "<Epidemilogical Distribution name>",
-  prob_distribution = "gamma",
-  prob_distribution_params = c(shape = 2.5, scale = 1.5)
+  prob_distribution = create_prob_distribution(
+    prob_distribution = "gamma",
+    prob_distribution_params = c(shape = 2.5, scale = 1.5)
+  )
 )
 convert_params_to_summary_stats(ep)
 ```