diff --git a/NAMESPACE b/NAMESPACE
index 89f5cf4..895c4c1 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,7 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
 export(detect_extinct)
-export(dist_setup)
 export(extinct_prob)
 export(inf_fn)
 export(outbreak_model)
@@ -18,7 +17,6 @@ importFrom(future.apply,future_lapply)
 importFrom(purrr,map2)
 importFrom(purrr,map2_dbl)
 importFrom(purrr,map_lgl)
-importFrom(purrr,partial)
 importFrom(purrr,safely)
 importFrom(sn,rsn)
 importFrom(stats,as.formula)
diff --git a/R/aux_functions.R b/R/aux_functions.R
index 4739274..a88f54e 100644
--- a/R/aux_functions.R
+++ b/R/aux_functions.R
@@ -1,25 +1,3 @@
-#' Create partial function to sample from gamma distributions
-#' @author Joel Hellewell
-#' @param dist_shape a positive `numeric` scalar: shape parameter of Weibull
-#'   distribution
-#' @param dist_scale a positive `numeric` scalar: scale parameter of Weibull
-#'   distribution
-#'
-#' @return partial function that takes a numeric argument for number of samples
-#' @export
-#' @importFrom purrr partial
-#'
-#' @examples
-#' incfn <- dist_setup(dist_shape = 2.32, dist_scale = 6.49)
-#' incfn(5)
-dist_setup <- function(dist_shape = NULL, dist_scale = NULL) {
-  out <- purrr::partial(rweibull,
-                 shape = dist_shape,
-                 scale = dist_scale)
-  return(out)
-}
-
-
 #' Samples the serial interval for given incubation period samples
 #'
 #' @param inc_samp a positive `numeric` vector: samples from the incubation
@@ -63,8 +41,8 @@ inf_fn <- function(inc_samp = NULL, k = NULL) {
 #'   r0community = 2.5,
 #'   disp.com = 0.16,
 #'   disp.iso = 1,
-#'   delay_shape = 1.65,
-#'   delay_scale = 4.28,
+#'   onset_to_isolation = \(x) stats::rweibull(n = x, shape = 1.65, scale = 4.28),
+#'   incubation_period = \(x) stats::rweibull(n = x, shape = 2.322737, scale = 6.492272),
 #'   k = 0,
 #'   quarantine = FALSE
 #' )
@@ -107,8 +85,8 @@ extinct_prob <- function(outbreak_df_week = NULL, cap_cases  = NULL, week_range
 #'   r0community = 2.5,
 #'   disp.com = 0.16,
 #'   disp.iso = 1,
-#'   delay_shape = 1.65,
-#'   delay_scale = 4.28,
+#'   onset_to_isolation = \(x) stats::rweibull(n = x, shape = 1.65, scale = 4.28),
+#'   incubation_period = \(x) stats::rweibull(n = x, shape = 2.322737, scale = 6.492272),
 #'   k = 0,
 #'   quarantine = FALSE
 #' )
diff --git a/R/outbreak_model.R b/R/outbreak_model.R
index cf71a56..03e8979 100644
--- a/R/outbreak_model.R
+++ b/R/outbreak_model.R
@@ -21,10 +21,14 @@
 #'   `TRUE` then traced contacts are isolated before symptom onset
 #' @param prop.asym a nonnegative `numeric` scalar: proportion of cases that
 #'   are completely asymptomatic (sublinical) (between 0 and 1)
-#' @param delay_shape a positive `numeric` scalar: shape parameter of delay
-#'   distribution
-#' @param delay_scale a positive `numeric` scalar: scale parameter of delay
-#'   distribution
+#' @param onset_to_isolation a `function`: a random number generating
+#'   `function` that accepts a single `integer` argument specifying the
+#'   length of the `function` output.
+#' @param incubation_period a `function`: a random number generating
+#'   `function` that samples from incubation period distribution, the
+#'   `function` accepts a single `integer` argument specifying the number of
+#'   times to sample the incubation period (i.e. length of the `function`
+#'   output).
 #' @param num.initial.cases a nonnegative `integer` scalar: number of initial
 #'   or starting cases which are all assumed to be missed.
 #' @param cap_cases a positive `integer` scalar: number of cumulative cases at
@@ -59,8 +63,8 @@
 #'   disp.com = 0.16,
 #'   disp.subclin = 0.16,
 #'   k = 0,
-#'   delay_shape = 1.65,
-#'   delay_scale = 4.28,
+#'   onset_to_isolation = \(x) stats::rweibull(n = x, shape = 1.65, scale = 4.28),
+#'   incubation_period = \(x) stats::rweibull(n = x, shape = 2.32, scale = 6.49),
 #'   prop.asym = 0,
 #'   quarantine = FALSE
 #' )
@@ -70,18 +74,8 @@ outbreak_model <- function(num.initial.cases = NULL, prop.ascertain = NULL,
                            r0isolated = NULL, r0community = NULL,
                            r0subclin = NULL, disp.iso = NULL,
                            disp.com = NULL, disp.subclin = NULL,
-                           k, delay_shape = NULL,
-                           delay_scale = NULL, prop.asym = NULL,
-                           quarantine = NULL) {
-
-  # Set up functions to sample from distributions
-  # incubation period sampling function
-  incfn <- dist_setup(dist_shape = 2.322737,
-                      dist_scale = 6.492272)
-  # incfn <- dist_setup(dist_shape = 3.303525,dist_scale = 6.68849) # incubation function for ECDC run
-  # onset to isolation delay sampling function
-  delayfn <- dist_setup(delay_shape,
-                        delay_scale)
+                           k, onset_to_isolation, incubation_period,
+                           prop.asym = NULL, quarantine = NULL) {
 
   # Set initial values for loop indices
   total.cases <- num.initial.cases
@@ -90,9 +84,9 @@ outbreak_model <- function(num.initial.cases = NULL, prop.ascertain = NULL,
 
   # Initial setup
   case_data <- outbreak_setup(num.initial.cases = num.initial.cases,
-                            incfn = incfn,
+                            incubation_period = incubation_period,
                             prop.asym = prop.asym,
-                            delayfn = delayfn,
+                            onset_to_isolation = onset_to_isolation,
                             k = k)
 
   # Preallocate
@@ -110,8 +104,8 @@ outbreak_model <- function(num.initial.cases = NULL, prop.ascertain = NULL,
                              r0isolated = r0isolated,
                              r0community = r0community,
                              r0subclin = r0subclin,
-                             incfn = incfn,
-                             delayfn = delayfn,
+                             incubation_period = incubation_period,
+                             onset_to_isolation = onset_to_isolation,
                              prop.ascertain = prop.ascertain,
                              k = k,
                              quarantine = quarantine,
diff --git a/R/outbreak_setup.R b/R/outbreak_setup.R
index 69b3ba4..690d8ce 100644
--- a/R/outbreak_setup.R
+++ b/R/outbreak_setup.R
@@ -21,20 +21,20 @@
 #'
 #' @examples
 #' # incubation period sampling function
-#' incfn <- dist_setup(dist_shape = 2.32, dist_scale = 6.49)
+#' incubation_period <- \(x) stats::rweibull(n = x, shape = 2.32, scale = 6.49)
 #' # delay distribution sampling function
-#' delayfn <- dist_setup(dist_shape = 1.65, dist_scale = 4.28)
+#' onset_to_isolation <- \(x) stats::rweibull(n = x, shape = 1.65, scale = 4.28)
 #' out <- outbreak_setup(
 #'   num.initial.cases = 1,
-#'   incfn = incfn,
-#'   delayfn = delayfn,
+#'   incubation_period = incubation_period,
+#'   onset_to_isolation = onset_to_isolation,
 #'   k = 1.95,
 #'   prop.asym = 0
 #' )
 #' out
-outbreak_setup <- function(num.initial.cases, incfn, delayfn, k, prop.asym) {
+outbreak_setup <- function(num.initial.cases, incubation_period, onset_to_isolation, k, prop.asym) {
   # Set up table of initial cases
-  inc_samples <- incfn(num.initial.cases)
+  inc_samples <- incubation_period(num.initial.cases)
 
   case_data <- data.table(exposure = rep(0, num.initial.cases), # Exposure time of 0 for all initial cases
                           asym = as.logical(rbinom(num.initial.cases, 1, prop.asym)),
@@ -45,7 +45,7 @@ outbreak_setup <- function(num.initial.cases, incfn, delayfn, k, prop.asym) {
                           new_cases = NA)
 
   # set isolation time for cluster to minimum time of onset of symptoms + draw from delay distribution
-  case_data <- case_data[, isolated_time := onset + delayfn(1)
+  case_data <- case_data[, isolated_time := onset + onset_to_isolation(1)
                          ][, isolated := FALSE]
 
   case_data$isolated_time[case_data$asym] <- Inf
diff --git a/R/outbreak_step.R b/R/outbreak_step.R
index e31adfe..7d00eef 100644
--- a/R/outbreak_step.R
+++ b/R/outbreak_step.R
@@ -5,10 +5,8 @@
 #' @param case_data a `data.table`: cases in outbreak so far; initially
 #'   generated by [outbreak_setup()]
 #' @inheritParams outbreak_model
-#' @param incfn a `function` that samples from incubation period Weibull
-#'   distribution; generated using [dist_setup()]
-#' @param delayfn a `function` that samples from the onset-to-hospitalisation
-#'   delay Weibull distribution; generated using [dist_setup()]
+#' @param onset_to_isolation a `function` that samples from the
+#'   onset-to-hospitalisation delay Weibull distribution
 #'
 #' @importFrom data.table data.table rbindlist
 #' @importFrom purrr map2 map2_dbl map_lgl
@@ -24,14 +22,14 @@
 #'
 #' @examples
 #' # incubation period sampling function
-#' incfn <- dist_setup(dist_shape = 2.32, dist_scale = 6.49)
+#' incubation_period <- \(x) stats::rweibull(n = x, shape = 2.32, scale = 6.49)
 #' # delay distribution sampling function
-#' delayfn <- dist_setup(dist_shape = 1.65, dist_scale = 4.28)
+#' onset_to_isolation <- \(x) stats::rweibull(n = x, shape = 1.65, scale = 4.28)
 #' # generate initial cases
 #' case_data <- outbreak_setup(
 #'   num.initial.cases = 5,
-#'   incfn = incfn,
-#'   delayfn = delayfn,
+#'   incubation_period = incubation_period,
+#'   onset_to_isolation = onset_to_isolation,
 #'   k = 1.95,
 #'   prop.asym = 0
 #' )
@@ -46,8 +44,8 @@
 #'   r0subclin = 1.25,
 #'   r0community = 2.5,
 #'   prop.asym = 0,
-#'   incfn = incfn,
-#'   delayfn = delayfn,
+#'   incubation_period = incubation_period,
+#'   onset_to_isolation = onset_to_isolation,
 #'   prop.ascertain = 0,
 #'   k = 1.95,
 #'   quarantine = FALSE
@@ -56,9 +54,10 @@
 #' case_data
 outbreak_step <- function(case_data = NULL, disp.iso = NULL, disp.com = NULL,
                           r0isolated = NULL, r0community = NULL,
-                          prop.asym = NULL, incfn = NULL, delayfn = NULL,
-                          prop.ascertain = NULL, k = NULL, quarantine = NULL,
-                          r0subclin = NULL, disp.subclin = NULL) {
+                          prop.asym = NULL, incubation_period = NULL,
+                          onset_to_isolation = NULL, prop.ascertain = NULL,
+                          k = NULL, quarantine = NULL, r0subclin = NULL,
+                          disp.subclin = NULL) {
 
   # For each case in case_data, draw new_cases from a negative binomial distribution
   # with an R0 and dispersion dependent on if isolated=TRUE
@@ -91,7 +90,7 @@ outbreak_step <- function(case_data = NULL, disp.iso = NULL, disp.com = NULL,
   }
 
   # Compile a data.table for all new cases, new_cases is the amount of people that each infector has infected
-  inc_samples <- incfn(total_new_cases)
+  inc_samples <- incubation_period(total_new_cases)
 
   prob_samples <- data.table(
     # time when new cases were exposed, a draw from serial interval based on infector's onset
@@ -137,11 +136,11 @@ outbreak_step <- function(case_data = NULL, disp.iso = NULL, disp.com = NULL,
   prob_samples[, isolated_time := ifelse(vect_isTRUE(asym), Inf,
                                         # If you are not asymptomatic, but you are missed,
                                         # you are isolated at your symptom onset
-                                        ifelse(vect_isTRUE(missed), onset + delayfn(1),
+                                        ifelse(vect_isTRUE(missed), onset + onset_to_isolation(1),
                                                # If you are not asymptomatic and you are traced,
                                                # you are isolated at max(onset,infector isolation time) # max(onset,infector_iso_time)
                                                ifelse(!vect_isTRUE(rep(quarantine, total_new_cases)),
-                                                      pmin(onset + delayfn(1), pmax(onset, infector_iso_time)),
+                                                      pmin(onset + onset_to_isolation(1), pmax(onset, infector_iso_time)),
                                                       infector_iso_time)))]
 
 
diff --git a/R/parameter_sweep.R b/R/parameter_sweep.R
index bd01410..05b192c 100644
--- a/R/parameter_sweep.R
+++ b/R/parameter_sweep.R
@@ -25,8 +25,10 @@
 #'   expand.grid(
 #'     delay_group = list(data.table(
 #'       delay = c("SARS", "Wuhan"),
-#'       delay_shape = c(1.651524, 2.305172),
-#'       delay_scale = c(4.287786, 9.483875)
+#'       onset_to_isolation = c(
+#'         \(x) stats::rweibull(n = x, shape = 1.651524, scale = 4.287786),
+#'         \(x) stats::rweibull(n = x, shape = 2.305172, scale = 9.483875)
+#'       )
 #'     )),
 #'     k_group = list(data.table(
 #'       theta = c("<1%", "15%", "30%"),
@@ -96,8 +98,7 @@ parameter_sweep <- function(scenarios = NULL, samples = 1,
       r0subclin = ifelse(
         "subclin_R0" %in% names(scenarios), x$subclin_R0, x$index_R0),
       k = x$k,
-      delay_shape = x$delay_shape,
-      delay_scale = x$delay_scale,
+      onset_to_isolation = x$onset_to_isolation,
       prop.ascertain = x$control_effectiveness,
       quarantine = x$quarantine,
       prop.asym = x$prop.asym
diff --git a/R/scenario_sim.R b/R/scenario_sim.R
index 74044c5..2870157 100644
--- a/R/scenario_sim.R
+++ b/R/scenario_sim.R
@@ -35,16 +35,17 @@
 #'   disp.iso = 1,
 #'   disp.com = 0.16,
 #'   k = 0.7,
-#'   delay_shape = 2.5,
-#'   delay_scale = 5,
+#'   onset_to_isolation = \(x) stats::rweibull(n = x, shape = 2.5, scale = 5),
+#'   incubation_period = \(x) stats::rweibull(n = x, shape = 2.32, scale = 6.49),
 #'   prop.asym = 0,
 #'   prop.ascertain = 0
 #' )
 #' res
 scenario_sim <- function(n.sim, prop.ascertain, cap_max_days, cap_cases,
                          r0isolated, r0community, disp.iso, disp.com, k,
-                         delay_shape, delay_scale, num.initial.cases, prop.asym,
-                         quarantine, r0subclin = NULL, disp.subclin = NULL) {
+                         onset_to_isolation, incubation_period,
+                         num.initial.cases, prop.asym, quarantine,
+                         r0subclin = NULL, disp.subclin = NULL) {
 
   # Set infectiousness of subclinical cases to be equal to clinical cases unless specified otherwise
   if(is.null(r0subclin)) {
@@ -65,8 +66,8 @@ scenario_sim <- function(n.sim, prop.ascertain, cap_max_days, cap_cases,
                                              disp.subclin = disp.subclin,
                                              disp.iso = disp.iso,
                                              disp.com = disp.com,
-                                             delay_shape = delay_shape,
-                                             delay_scale = delay_scale,
+                                             onset_to_isolation = onset_to_isolation,
+                                             incubation_period = incubation_period,
                                              k = k,
                                              prop.asym = prop.asym,
                                              quarantine = quarantine))
diff --git a/README.Rmd b/README.Rmd
index 791400d..2d6cea1 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -55,8 +55,8 @@ res <- scenario_sim(
   r0community = 2.5, ## non-isolated individuals have R0 of 2.5
   disp.com = 0.16, ## dispersion parameter in the community
   disp.iso = 1, ## dispersion  parameter of those isolated
-  delay_shape = 1.651524, ## shape parameter of time from onset to isolation
-  delay_scale = 4.287786, ## scale parameter of time from onset to isolation
+  onset_to_isolation = \(x) stats::rweibull(n = x, shape = 1.651524, scale = 4.287786), ## time from onset to isolation
+  incubation_period = \(x) stats::rweibull(n = x, shape = 2.322737, scale = 6.492272), ## incubation period
   k = 0, ## skew of generation interval to be beyond onset of symptoms
   quarantine = FALSE ## whether quarantine is in effect
 )
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 36d2d69..77f4317 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -7,7 +7,6 @@ reference:
   - title: Helper functions
   - subtitle: Probability distribution manipulation
     contents:
-    - dist_setup
     - inf_fn
   - subtitle: Loop wrappers for scenario modelling
     contents:
diff --git a/man/detect_extinct.Rd b/man/detect_extinct.Rd
index cffe477..cce7710 100644
--- a/man/detect_extinct.Rd
+++ b/man/detect_extinct.Rd
@@ -36,8 +36,8 @@ res <- scenario_sim(
   r0community = 2.5,
   disp.com = 0.16,
   disp.iso = 1,
-  delay_shape = 1.65,
-  delay_scale = 4.28,
+  onset_to_isolation = \(x) stats::rweibull(n = x, shape = 1.65, scale = 4.28),
+  incubation_period = \(x) stats::rweibull(n = x, shape = 2.322737, scale = 6.492272),
   k = 0,
   quarantine = FALSE
 )
diff --git a/man/dist_setup.Rd b/man/dist_setup.Rd
deleted file mode 100644
index 7a0f4dd..0000000
--- a/man/dist_setup.Rd
+++ /dev/null
@@ -1,28 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/aux_functions.R
-\name{dist_setup}
-\alias{dist_setup}
-\title{Create partial function to sample from gamma distributions}
-\usage{
-dist_setup(dist_shape = NULL, dist_scale = NULL)
-}
-\arguments{
-\item{dist_shape}{a positive \code{numeric} scalar: shape parameter of Weibull
-distribution}
-
-\item{dist_scale}{a positive \code{numeric} scalar: scale parameter of Weibull
-distribution}
-}
-\value{
-partial function that takes a numeric argument for number of samples
-}
-\description{
-Create partial function to sample from gamma distributions
-}
-\examples{
-incfn <- dist_setup(dist_shape = 2.32, dist_scale = 6.49)
-incfn(5)
-}
-\author{
-Joel Hellewell
-}
diff --git a/man/extinct_prob.Rd b/man/extinct_prob.Rd
index 3a6aa37..ab71c48 100644
--- a/man/extinct_prob.Rd
+++ b/man/extinct_prob.Rd
@@ -34,8 +34,8 @@ res <- scenario_sim(
   r0community = 2.5,
   disp.com = 0.16,
   disp.iso = 1,
-  delay_shape = 1.65,
-  delay_scale = 4.28,
+  onset_to_isolation = \(x) stats::rweibull(n = x, shape = 1.65, scale = 4.28),
+  incubation_period = \(x) stats::rweibull(n = x, shape = 2.322737, scale = 6.492272),
   k = 0,
   quarantine = FALSE
 )
diff --git a/man/outbreak_model.Rd b/man/outbreak_model.Rd
index 49b0ad3..65ede0d 100644
--- a/man/outbreak_model.Rd
+++ b/man/outbreak_model.Rd
@@ -16,8 +16,8 @@ outbreak_model(
   disp.com = NULL,
   disp.subclin = NULL,
   k,
-  delay_shape = NULL,
-  delay_scale = NULL,
+  onset_to_isolation,
+  incubation_period,
   prop.asym = NULL,
   quarantine = NULL
 )
@@ -56,11 +56,15 @@ sub-clincial non-isolated cases (must be >0)}
 \item{k}{a \code{numeric} scalar: skew parameter for sampling the serial
 interval from the incubation period}
 
-\item{delay_shape}{a positive \code{numeric} scalar: shape parameter of delay
-distribution}
+\item{onset_to_isolation}{a \code{function}: a random number generating
+\code{function} that accepts a single \code{integer} argument specifying the
+length of the \code{function} output.}
 
-\item{delay_scale}{a positive \code{numeric} scalar: scale parameter of delay
-distribution}
+\item{incubation_period}{a \code{function}: a random number generating
+\code{function} that samples from incubation period distribution, the
+\code{function} accepts a single \code{integer} argument specifying the number of
+times to sample the incubation period (i.e. length of the \code{function}
+output).}
 
 \item{prop.asym}{a nonnegative \code{numeric} scalar: proportion of cases that
 are completely asymptomatic (sublinical) (between 0 and 1)}
@@ -96,8 +100,8 @@ out <- outbreak_model(
   disp.com = 0.16,
   disp.subclin = 0.16,
   k = 0,
-  delay_shape = 1.65,
-  delay_scale = 4.28,
+  onset_to_isolation = \(x) stats::rweibull(n = x, shape = 1.65, scale = 4.28),
+  incubation_period = \(x) stats::rweibull(n = x, shape = 2.32, scale = 6.49),
   prop.asym = 0,
   quarantine = FALSE
 )
diff --git a/man/outbreak_setup.Rd b/man/outbreak_setup.Rd
index 2721261..c72ba3f 100644
--- a/man/outbreak_setup.Rd
+++ b/man/outbreak_setup.Rd
@@ -4,17 +4,27 @@
 \alias{outbreak_setup}
 \title{Set up initial cases for branching process}
 \usage{
-outbreak_setup(num.initial.cases, incfn, delayfn, k, prop.asym)
+outbreak_setup(
+  num.initial.cases,
+  incubation_period,
+  onset_to_isolation,
+  k,
+  prop.asym
+)
 }
 \arguments{
 \item{num.initial.cases}{a nonnegative \code{integer} scalar: number of initial
 or starting cases which are all assumed to be missed.}
 
-\item{incfn}{a \code{function} that samples from incubation period Weibull
-distribution; generated using \code{\link[=dist_setup]{dist_setup()}}}
+\item{incubation_period}{a \code{function}: a random number generating
+\code{function} that samples from incubation period distribution, the
+\code{function} accepts a single \code{integer} argument specifying the number of
+times to sample the incubation period (i.e. length of the \code{function}
+output).}
 
-\item{delayfn}{a \code{function} that samples from the onset-to-hospitalisation
-delay Weibull distribution; generated using \code{\link[=dist_setup]{dist_setup()}}}
+\item{onset_to_isolation}{a \code{function}: a random number generating
+\code{function} that accepts a single \code{integer} argument specifying the
+length of the \code{function} output.}
 
 \item{k}{a \code{numeric} scalar: skew parameter for sampling the serial
 interval from the incubation period}
@@ -41,13 +51,13 @@ Set up initial cases for branching process
 }
 \examples{
 # incubation period sampling function
-incfn <- dist_setup(dist_shape = 2.32, dist_scale = 6.49)
+incubation_period <- \(x) stats::rweibull(n = x, shape = 2.32, scale = 6.49)
 # delay distribution sampling function
-delayfn <- dist_setup(dist_shape = 1.65, dist_scale = 4.28)
+onset_to_isolation <- \(x) stats::rweibull(n = x, shape = 1.65, scale = 4.28)
 out <- outbreak_setup(
   num.initial.cases = 1,
-  incfn = incfn,
-  delayfn = delayfn,
+  incubation_period = incubation_period,
+  onset_to_isolation = onset_to_isolation,
   k = 1.95,
   prop.asym = 0
 )
diff --git a/man/outbreak_step.Rd b/man/outbreak_step.Rd
index 1ecb26d..dd0969d 100644
--- a/man/outbreak_step.Rd
+++ b/man/outbreak_step.Rd
@@ -11,8 +11,8 @@ outbreak_step(
   r0isolated = NULL,
   r0community = NULL,
   prop.asym = NULL,
-  incfn = NULL,
-  delayfn = NULL,
+  incubation_period = NULL,
+  onset_to_isolation = NULL,
   prop.ascertain = NULL,
   k = NULL,
   quarantine = NULL,
@@ -39,11 +39,14 @@ non-isolated cases (must be >0)}
 \item{prop.asym}{a nonnegative \code{numeric} scalar: proportion of cases that
 are completely asymptomatic (sublinical) (between 0 and 1)}
 
-\item{incfn}{a \code{function} that samples from incubation period Weibull
-distribution; generated using \code{\link[=dist_setup]{dist_setup()}}}
+\item{incubation_period}{a \code{function}: a random number generating
+\code{function} that samples from incubation period distribution, the
+\code{function} accepts a single \code{integer} argument specifying the number of
+times to sample the incubation period (i.e. length of the \code{function}
+output).}
 
-\item{delayfn}{a \code{function} that samples from the onset-to-hospitalisation
-delay Weibull distribution; generated using \code{\link[=dist_setup]{dist_setup()}}}
+\item{onset_to_isolation}{a \code{function} that samples from the
+onset-to-hospitalisation delay Weibull distribution}
 
 \item{prop.ascertain}{a nonnegative \code{numeric} scalar: proportion of
 infectious contacts ascertained by contact tracing (must be 0<=x<=1)}
@@ -74,14 +77,14 @@ Move forward one generation in the branching process
 }
 \examples{
 # incubation period sampling function
-incfn <- dist_setup(dist_shape = 2.32, dist_scale = 6.49)
+incubation_period <- \(x) stats::rweibull(n = x, shape = 2.32, scale = 6.49)
 # delay distribution sampling function
-delayfn <- dist_setup(dist_shape = 1.65, dist_scale = 4.28)
+onset_to_isolation <- \(x) stats::rweibull(n = x, shape = 1.65, scale = 4.28)
 # generate initial cases
 case_data <- outbreak_setup(
   num.initial.cases = 5,
-  incfn = incfn,
-  delayfn = delayfn,
+  incubation_period = incubation_period,
+  onset_to_isolation = onset_to_isolation,
   k = 1.95,
   prop.asym = 0
 )
@@ -96,8 +99,8 @@ out <- outbreak_step(
   r0subclin = 1.25,
   r0community = 2.5,
   prop.asym = 0,
-  incfn = incfn,
-  delayfn = delayfn,
+  incubation_period = incubation_period,
+  onset_to_isolation = onset_to_isolation,
   prop.ascertain = 0,
   k = 1.95,
   quarantine = FALSE
diff --git a/man/parameter_sweep.Rd b/man/parameter_sweep.Rd
index fcfa0b9..1691619 100644
--- a/man/parameter_sweep.Rd
+++ b/man/parameter_sweep.Rd
@@ -33,8 +33,10 @@ scenarios <- data.table(
   expand.grid(
     delay_group = list(data.table(
       delay = c("SARS", "Wuhan"),
-      delay_shape = c(1.651524, 2.305172),
-      delay_scale = c(4.287786, 9.483875)
+      onset_to_isolation = c(
+        \(x) stats::rweibull(n = x, shape = 1.651524, scale = 4.287786),
+        \(x) stats::rweibull(n = x, shape = 2.305172, scale = 9.483875)
+      )
     )),
     k_group = list(data.table(
       theta = c("<1\%", "15\%", "30\%"),
diff --git a/man/scenario_sim.Rd b/man/scenario_sim.Rd
index 0ce0142..9076b33 100644
--- a/man/scenario_sim.Rd
+++ b/man/scenario_sim.Rd
@@ -14,8 +14,8 @@ scenario_sim(
   disp.iso,
   disp.com,
   k,
-  delay_shape,
-  delay_scale,
+  onset_to_isolation,
+  incubation_period,
   num.initial.cases,
   prop.asym,
   quarantine,
@@ -50,11 +50,15 @@ non-isolated cases (must be >0)}
 \item{k}{a \code{numeric} scalar: skew parameter for sampling the serial
 interval from the incubation period}
 
-\item{delay_shape}{a positive \code{numeric} scalar: shape parameter of delay
-distribution}
+\item{onset_to_isolation}{a \code{function}: a random number generating
+\code{function} that accepts a single \code{integer} argument specifying the
+length of the \code{function} output.}
 
-\item{delay_scale}{a positive \code{numeric} scalar: scale parameter of delay
-distribution}
+\item{incubation_period}{a \code{function}: a random number generating
+\code{function} that samples from incubation period distribution, the
+\code{function} accepts a single \code{integer} argument specifying the number of
+times to sample the incubation period (i.e. length of the \code{function}
+output).}
 
 \item{num.initial.cases}{a nonnegative \code{integer} scalar: number of initial
 or starting cases which are all assumed to be missed.}
@@ -97,8 +101,8 @@ res <- scenario_sim(
   disp.iso = 1,
   disp.com = 0.16,
   k = 0.7,
-  delay_shape = 2.5,
-  delay_scale = 5,
+  onset_to_isolation = \(x) stats::rweibull(n = x, shape = 2.5, scale = 5),
+  incubation_period = \(x) stats::rweibull(n = x, shape = 2.32, scale = 6.49),
   prop.asym = 0,
   prop.ascertain = 0
 )
diff --git a/tests/testthat/test-aux_functions.R b/tests/testthat/test-aux_functions.R
index 43886ae..04313ce 100644
--- a/tests/testthat/test-aux_functions.R
+++ b/tests/testthat/test-aux_functions.R
@@ -2,32 +2,6 @@ context("Test all aux functions")
 
 set.seed(515)
 
-test_that("dist_setup returns a partialised function", {
-  r1 <- dist_setup(0.1, 2)
-  expect_error(r1(20), NA)
-  expect_true(inherits(r1, 'purrr_function_partial'))
-})
-
-test_that("dist_setup parameters behave as expected", {
-  f1 <- dist_setup(1, 0.01)
-  f2 <- dist_setup(1, 1e5)
-  r1 <- f1(1e5)
-  r2 <- f2(1e5)
-
-  expect_lt(mean(r1), mean(r2))
-  expect_true(inherits(f1, 'purrr_function_partial'))
-  expect_true(is.numeric(r1))
-  expect_length(r1, 1e5)
-
-  f3 <- dist_setup(0.01, 10)
-  r3 <- f3(1e5)
-  expect_gt(mean(r3), median(r3))
-
-  f4 <- dist_setup(50, 10)
-  r4 <- f4(1e5)
-  expect_lt(mean(r4), median(r4))
-})
-
 test_that("inf_fn parameters behave as expected", {
   r1 <- inf_fn(c(1, 4, 1), 2)
   expect_length(r1, 3)
@@ -73,8 +47,8 @@ test_that('extinct_prob works as expected', {
     disp.iso = 1,
     disp.com = 0.16,
     k = 0.7,
-    delay_shape = 2.5,
-    delay_scale = 5,
+    onset_to_isolation = \(x) stats::rweibull(n = x, shape = 2.5, scale = 5),
+    incubation_period = \(x) stats::rweibull(n = x, shape = 2.322737, scale = 6.492272),
     prop.asym = 0,
     prop.ascertain = 0
   )
@@ -117,8 +91,8 @@ test_that('extinct_prob works as expected', {
     disp.iso = 1,
     disp.com = 0.16,
     k = 0.7,
-    delay_shape = 2.5,
-    delay_scale = 5,
+    onset_to_isolation = \(x) stats::rweibull(n = x, shape = 2.5, scale = 5),
+    incubation_period = \(x) stats::rweibull(n = x, shape = 2.322737, scale = 6.492272),
     prop.asym = 0,
     prop.ascertain = 0
   )
@@ -137,8 +111,8 @@ test_that('extinct_prob works as expected', {
     disp.iso = 1,
     disp.com = 0.16,
     k = 0.7,
-    delay_shape = 2.5,
-    delay_scale = 5,
+    onset_to_isolation = \(x) stats::rweibull(n = x, shape = 2.5, scale = 5),
+    incubation_period = \(x) stats::rweibull(n = x, shape = 2.322737, scale = 6.492272),
     prop.asym = 0,
     prop.ascertain = 0
   )
@@ -160,8 +134,8 @@ test_that('extinct_prob week_range argument works', {
     disp.iso = 1,
     disp.com = 0.16,
     k = 0.7,
-    delay_shape = 2.5,
-    delay_scale = 5,
+    onset_to_isolation = \(x) stats::rweibull(n = x, shape = 2.5, scale = 5),
+    incubation_period = \(x) stats::rweibull(n = x, shape = 2.322737, scale = 6.492272),
     prop.asym = 0,
     prop.ascertain = 0
   )
@@ -209,7 +183,7 @@ test_that('extinct_prob week_range argument works', {
 
   r5b <- extinct_prob(res5, cap_cases = cap, week_range = 2)
   expect_equal(r5b, 1)
-  
+
   # Case of cases in week 2 but not 1 (by all sensible definitions is not an
   # extinction). Test here that week_range 1:2 says no extinction and neither
   # does week_range 2.
@@ -236,8 +210,8 @@ test_that('detect_extinct works', {
     disp.iso = 1,
     disp.com = 0.16,
     k = 0.7,
-    delay_shape = 2.5,
-    delay_scale = 5,
+    onset_to_isolation = \(x) stats::rweibull(n = x, shape = 2.5, scale = 5),
+    incubation_period = \(x) stats::rweibull(n = x, shape = 2.322737, scale = 6.492272),
     prop.asym = 0,
     prop.ascertain = 0
   )
@@ -296,7 +270,7 @@ test_that('detect_extinct works', {
 
   expect5 <- data.table(sim = c(1.0), extinct = c(0.0))
   expect_equal(r5, expect5)
-  
+
   r5b <- detect_extinct(res5, cap_cases = cap, week_range = 2)
   # The types in the output is a bit random. So just force all to doubles.
   r5b <- data.table(r5b)[, lapply(.SD, as.double)]
diff --git a/tests/testthat/test-basicusage.R b/tests/testthat/test-basicusage.R
index b701f4d..6ccda94 100644
--- a/tests/testthat/test-basicusage.R
+++ b/tests/testthat/test-basicusage.R
@@ -4,7 +4,9 @@ set.seed(20200410)
 res <- ringbp::scenario_sim(
   n.sim = 2, num.initial.cases = 10, prop.asym = 0, prop.ascertain = 0.2,
   cap_cases = 20, cap_max_days = 100, r0isolated = 0, r0community = 2.5,
-  disp.com = 0.16, disp.iso = 1, delay_shape = 1.651524, delay_scale = 4.287786,
+  disp.com = 0.16, disp.iso = 1,
+  onset_to_isolation = \(x) stats::rweibull(n = x, shape = 1.651524, scale = 4.287786),
+  incubation_period = \(x) stats::rweibull(n = x, shape = 2.322737, scale = 6.492272),
   k = 0, quarantine = FALSE)
 
 test_that("A basic sim returns the correct object", {
diff --git a/tests/testthat/test-outbreak_setup.R b/tests/testthat/test-outbreak_setup.R
index 5933437..ea0178b 100644
--- a/tests/testthat/test-outbreak_setup.R
+++ b/tests/testthat/test-outbreak_setup.R
@@ -3,14 +3,14 @@ context("Test basic usage")
 set.seed(20200410)
 
 test_that("A basic sim setup returns the correct object", {
-  incfn <- dist_setup(dist_shape = 2.322737, dist_scale = 6.492272)
+  incubation_period <- \(x) stats::rweibull(n = x, shape = 2.322737, scale = 6.492272)
   # delay distribution sampling function
-  delayfn <- dist_setup(2, 4)
+  onset_to_isolation <- \(x) stats::rweibull(n = x, shape = 2, scale = 4)
   # generate initial cases
   case_data <- outbreak_setup(
     num.initial.cases = 5,
-    incfn = incfn,
-    delayfn = delayfn,
+    incubation_period = incubation_period,
+    onset_to_isolation = onset_to_isolation,
     k = 1.95,
     prop.asym = 0
   )
@@ -21,15 +21,15 @@ test_that("A basic sim setup returns the correct object", {
 })
 
 test_that("asym arg works properly", {
-  incfn <- dist_setup(dist_shape = 2.322737, dist_scale = 6.492272)
+  incubation_period <- \(x) stats::rweibull(n = x, shape = 2.322737, scale = 6.492272)
   # delay distribution sampling function
-  delayfn <- dist_setup(2, 4)
+  onset_to_isolation <- \(x) stats::rweibull(n = x, shape = 2, scale = 4)
   # generate initial cases
   # All asymptomatics
   all_asym <- outbreak_setup(
     num.initial.cases = 5,
-    incfn = incfn,
-    delayfn = delayfn,
+    incubation_period = incubation_period,
+    onset_to_isolation = onset_to_isolation,
     k = 1.95,
     prop.asym = 1
   )
@@ -40,8 +40,8 @@ test_that("asym arg works properly", {
   # machine precision
   mix <- outbreak_setup(
     num.initial.cases = 10000,
-    incfn = incfn,
-    delayfn  =  delayfn,
+    incubation_period = incubation_period,
+    onset_to_isolation = onset_to_isolation,
     k = 1.95,
     prop.asym = 0.5
   )
diff --git a/tests/testthat/test-outbreak_step.R b/tests/testthat/test-outbreak_step.R
index db85910..638921a 100644
--- a/tests/testthat/test-outbreak_step.R
+++ b/tests/testthat/test-outbreak_step.R
@@ -3,14 +3,14 @@ context("Test basic usage")
 set.seed(123456)
 
 test_that("A basic sim returns the correct object", {
-  incfn <- dist_setup(dist_shape = 2.322737,dist_scale = 6.492272)
+  incubation_period <- \(x) stats::rweibull(n = x, shape = 2.322737, scale = 6.492272)
   # delay distribution sampling function
-  delayfn <- dist_setup(2, 4)
+  onset_to_isolation <- \(x) stats::rweibull(n = x, shape = 2, scale = 4)
   # generate initial cases
   case_data <- outbreak_setup(
     num.initial.cases = 1,
-    incfn = incfn,
-    delayfn  =  delayfn,
+    incubation_period = incubation_period,
+    onset_to_isolation  =  onset_to_isolation,
     k = 1.95,
     prop.asym = 0
   )
@@ -23,8 +23,8 @@ test_that("A basic sim returns the correct object", {
     r0isolated = 0,
     r0community = 500, # almost guarentees to get new cases
     prop.asym = 0,
-    incfn = incfn,
-    delayfn = delayfn,
+    incubation_period = incubation_period,
+    onset_to_isolation = onset_to_isolation,
     prop.ascertain = 0,
     k = 1.95,
     quarantine = FALSE
@@ -44,8 +44,8 @@ test_that("A basic sim returns the correct object", {
     r0isolated = 0,
     r0community = 0, # almost guarentees to get new cases
     prop.asym = 0,
-    incfn = incfn,
-    delayfn = delayfn,
+    incubation_period = incubation_period,
+    onset_to_isolation = onset_to_isolation,
     prop.ascertain = 0,
     k = 1.95,
     quarantine = FALSE
@@ -55,14 +55,14 @@ test_that("A basic sim returns the correct object", {
 })
 
 test_that("Sim with multiple infectors makes senes", {
-  incfn <- dist_setup(dist_shape = 2.322737, dist_scale = 6.492272)
+  incubation_period <- \(x) stats::rweibull(n = x, shape = 2.322737, scale = 6.492272)
   # delay distribution sampling function
-  delayfn <- dist_setup(2, 4)
+  onset_to_isolation <- \(x) stats::rweibull(n = x, shape = 2, scale = 4)
   # generate initial cases
   case_data <- outbreak_setup(
     num.initial.cases = 2,
-    incfn = incfn,
-    delayfn = delayfn,
+    incubation_period = incubation_period,
+    onset_to_isolation = onset_to_isolation,
     k = 1.95,
     prop.asym = 0
   )
@@ -75,8 +75,8 @@ test_that("Sim with multiple infectors makes senes", {
     r0isolated = 0,
     r0community = 10000, # almost guarentees both index cases create infections
     prop.asym = 0,
-    incfn = incfn,
-    delayfn = delayfn,
+    incubation_period = incubation_period,
+    onset_to_isolation = onset_to_isolation,
     prop.ascertain = 0,
     k = 1.95,
     quarantine = FALSE
@@ -89,14 +89,14 @@ test_that("Sim with multiple infectors makes senes", {
 
 
 test_that("R0isolated is working properly", {
-  incfn <- dist_setup(dist_shape = 2.322737,dist_scale = 6.492272)
+  incubation_period <- \(x) stats::rweibull(n = x, shape = 2.322737, scale = 6.492272)
   # delay distribution sampling function
-  delayfn <- dist_setup(2, 4)
+  onset_to_isolation <- \(x) stats::rweibull(n = x, shape = 2, scale = 4)
   # generate initial cases
   case_data <- outbreak_setup(
     num.initial.cases = 1,
-    incfn = incfn,
-    delayfn = delayfn,
+    incubation_period = incubation_period,
+    onset_to_isolation = onset_to_isolation,
     k = 1.95,
     prop.asym = 0
   )
@@ -110,8 +110,8 @@ test_that("R0isolated is working properly", {
     r0isolated = 0, # Shoiuld get zero cases
     r0community = 500, # Case is isolated so irrelevent
     prop.asym = 0,
-    incfn = incfn,
-    delayfn = delayfn,
+    incubation_period = incubation_period,
+    onset_to_isolation = onset_to_isolation,
     prop.ascertain = 0,
     k = 1.95,
     quarantine = FALSE
@@ -126,8 +126,8 @@ test_that("R0isolated is working properly", {
     r0isolated = 500, # Shoiuld get lots of cases
     r0community = 0, # Case is isolated so irrelevent
     prop.asym = 0,
-    incfn = incfn,
-    delayfn = delayfn,
+    incubation_period = incubation_period,
+    onset_to_isolation = onset_to_isolation,
     prop.ascertain = 0,
     k = 1.95,
     quarantine = FALSE
@@ -136,14 +136,14 @@ test_that("R0isolated is working properly", {
 })
 
 test_that('Test a bunch of args',{
-  incfn <- dist_setup(dist_shape = 2.322737,dist_scale = 6.492272)
+  incubation_period <- \(x) stats::rweibull(n = x, shape = 2.322737, scale = 6.492272)
   # delay distribution sampling function
-  delayfn <- dist_setup(2, 4)
+  onset_to_isolation <- \(x) stats::rweibull(n = x, shape = 2, scale = 4)
   # generate initial cases
   case_data <- outbreak_setup(
     num.initial.cases = 1,
-    incfn = incfn,
-    delayfn = delayfn,
+    incubation_period = incubation_period,
+    onset_to_isolation = onset_to_isolation,
     k = 1.95,
     prop.asym = 0
   )
@@ -156,8 +156,8 @@ test_that('Test a bunch of args',{
     r0isolated = 0,
     r0community = 10000, # almost guarentees both index cases create infections
     prop.asym = 0,
-    incfn = incfn,
-    delayfn = delayfn,
+    incubation_period = incubation_period,
+    onset_to_isolation = onset_to_isolation,
     prop.ascertain = 0,
     k = 1.95,
     quarantine = FALSE
@@ -171,8 +171,8 @@ test_that('Test a bunch of args',{
     r0isolated = 0,
     r0community = 10000, # almost guarentees both index cases create infections
     prop.asym = 0,
-    incfn = incfn,
-    delayfn = delayfn,
+    incubation_period = incubation_period,
+    onset_to_isolation = onset_to_isolation,
     prop.ascertain = 1,
     k = 1.95,
     quarantine = FALSE
@@ -187,8 +187,8 @@ test_that('Test a bunch of args',{
     r0isolated = 0,
     r0community = 100000, # To test a mix make sure there's loads of cases.
     prop.asym = 0,
-    incfn = incfn,
-    delayfn = delayfn,
+    incubation_period = incubation_period,
+    onset_to_isolation = onset_to_isolation,
     prop.ascertain = 0.5,
     k = 1.95,
     quarantine = FALSE