diff --git a/R/dist.R b/R/dist.R
index 5594ea20f..80c334ddc 100644
--- a/R/dist.R
+++ b/R/dist.R
@@ -84,50 +84,30 @@
 #' )
 dist_skel <- function(n, dist = FALSE, cum = TRUE, model,
                       discrete = FALSE, params, max_value = 120) {
+  ## define unnormalised support function
   if (model == "exp") {
     # define support functions for exponential dist
-    rdist <- function(n) {
-      rexp(n, params$rate)
-    }
-    pdist <- function(n) {
-      pexp(n, params$rate) / pexp(max_value, params$rate)
-    }
-    ddist <- function(n) {
-      (pexp(n + 1, params$rate) -
-        pexp(n, params$rate)) /
-        pexp(max_value + 1, params$rate)
+    updist <- function(n) {
+      pexp(n, params$rate)
     }
   } else if (model == "gamma") {
-    rdist <- function(n) {
-      rgamma(n, params$shape, params$scale)
-    }
-    pdist <- function(n) {
-      pgamma(n, params$shape, params$scale) /
-        pgamma(max_value + 1, params$shape, params$scale)
-    }
-    ddist <- function(n) {
-      (pgamma(n + 1, params$shape, params$scale) -
-        pgamma(n, params$shape, params$scale)) /
-        pgamma(max_value + 1, params$shape, params$scale)
+    # define support functions for gamma dist
+    updist <- function(n) {
+      pgamma(n, params$shape, params$scale)
     }
   } else if (model == "lognormal") {
-    rdist <- function(n) {
-      rlnorm(n, params$mean, params$sd)
-    }
-    pdist <- function(n) {
-      plnorm(n, params$mean, params$sd) /
-        plnorm(max_value + 1, params$mean, params$sd)
-    }
-    ddist <- function(n) {
-      (plnorm(n + 1, params$mean, params$sd) -
-        plnorm(n, params$mean, params$sd)) /
-        plnorm(max_value + 1, params$mean, params$sd)
+    updist <- function(n) {
+      plnorm(n, params$mean, params$sd)
     }
   }
 
+  ## define normalised support functions
   if (discrete) {
-    cmf <- c(0, pdist(seq_len(max_value + 2)))
-    pmf <- (cmf[-1] - c(0, cmf[seq(1, length(cmf) - 2)])) / sum(tail(cmf, 2))
+    cmf <- c(0, updist(1),
+      updist(seq_len(max_value)) + updist(seq_len(max_value) + 1)
+    ) /
+      (updist(max_value) + updist(max_value + 1))
+    pmf <- diff(cmf)
     rdist <- function(n) {
       sample(x = seq_len(max_value + 1) - 1, size = n, prob = pmf)
     }
@@ -137,6 +117,26 @@ dist_skel <- function(n, dist = FALSE, cum = TRUE, model,
     ddist <- function(n) {
       pmf[n + 1]
     }
+  } else {
+    pdist <- function(n) {
+      updist(n) / updist(max_value + 1)
+    }
+    ddist <- function(n) {
+      pdist(n + 1) - pdist(n)
+    }
+    if (model == "exp") {
+      rdist <- function(n) {
+        rexp(n, params$rate)
+      }
+    } else if (model == "gamma") {
+      rdist <- function(n) {
+        rgamma(n, params$shape, params$scale)
+      }
+    } else if (model == "lognormal") {
+      rdist <- function(n) {
+        rlnorm(n, params$mean, params$sd)
+      }
+    }
   }
 
   # define internal sampling function
@@ -161,7 +161,6 @@ dist_skel <- function(n, dist = FALSE, cum = TRUE, model,
       while (any(!is.na(n) & n >= max_value)) {
         n <- ifelse(n >= max_value, inner_skel(n), n)
       }
-
       n <- as.integer(n)
     }
     return(n)