diff --git a/NAMESPACE b/NAMESPACE
index 989308e..93b5927 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,7 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(primary_censored_cdf,default)
-S3method(primary_censored_cdf,pcens_numeric)
 S3method(primary_censored_cdf,pcens_pgamma_dunif)
 S3method(primary_censored_cdf,pcens_plnorm_dunif)
 export(check_dprimary)
diff --git a/NEWS.md b/NEWS.md
index 626bbc2..553aa86 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -7,6 +7,7 @@ This is the current development version.
 * Add `{touchstone}` based benchmarks for benchmarking R utility functions, and fitting the `stan` and `fitdistplus` models.
 * Added a "How it works" vignette.
 * Added R infrastructure for analytical solutions via the `primary_censored_dist` S3 class.
+* Added analytical solutions for the gamma and lognormal distributions with uniform primary censoring.
 
 # primarycensoreddist 0.4.0
 
diff --git a/R/fitdistdoublecens.R b/R/fitdistdoublecens.R
index 72631d7..16d1f8d 100644
--- a/R/fitdistdoublecens.R
+++ b/R/fitdistdoublecens.R
@@ -55,6 +55,7 @@
 fitdistdoublecens <- function(censdata, distr,
                               pwindow = 1, D = Inf,
                               dprimary = stats::dunif,
+                              dprimary_name = NULL,
                               dprimary_args = list(), ...) {
   # Check if fitdistrplus is available
   if (!requireNamespace("fitdistrplus", quietly = TRUE)) {
@@ -67,9 +68,13 @@ fitdistdoublecens <- function(censdata, distr,
     stop("censdata must contain 'left' and 'right' columns")
   }
 
-  # Get the distribution functions
-  pdist <- get(paste0("p", distr))
+  if (is.null(dprimary_name)) {
+    dprimary_name <- .extract_function_name(dprimary)
+  }
 
+  # Get the distribution functions
+  pdist_name <- paste0("p", distr)
+  pdist <- get(pdist_name)
   swindows <- censdata$right - censdata$left
 
   # Create the function definition with named arguments for dpcens
@@ -83,7 +88,9 @@ fitdistdoublecens <- function(censdata, distr,
         pwindow = pwindow,
         D = D,
         dprimary = dprimary,
-        dprimary_args = dprimary_args
+        dprimary_args = dprimary_args,
+        pdist_name = pdist_name,
+        dprimary_name = dprimary_name
       )
     ))
   }
@@ -99,7 +106,9 @@ fitdistdoublecens <- function(censdata, distr,
         pwindow = pwindow,
         D = D,
         dprimary = dprimary,
-        dprimary_args = dprimary_args
+        dprimary_args = dprimary_args,
+        pdist_name = pdist_name,
+        dprimary_name = dprimary_name
       )
     ))
   }
@@ -125,14 +134,15 @@ fitdistdoublecens <- function(censdata, distr,
 #' each element in x
 #' @keywords internal
 .dpcens <- function(x, swindows, pdist, pwindow, D, dprimary,
-                    dprimary_args, ...) {
+                    dprimary_args, pdist_name, dprimary_name, ...) {
   tryCatch(
     {
       if (length(unique(swindows)) == 1) {
         dprimarycensoreddist(
           x, pdist,
           pwindow = pwindow, swindow = swindows[1], D = D, dprimary = dprimary,
-          dprimary_args = dprimary_args, ...
+          dprimary_args = dprimary_args, pdist_name = pdist_name,
+          dprimary_name = dprimary_name, ...
         )
       } else {
         # Group x and swindows by unique swindow values
@@ -144,7 +154,8 @@ fitdistdoublecens <- function(censdata, distr,
           result[mask] <- dprimarycensoreddist(
             x[mask], pdist,
             pwindow = pwindow, swindow = sw, D = D,
-            dprimary = dprimary, dprimary_args = dprimary_args, ...
+            dprimary = dprimary, dprimary_args = dprimary_args,
+            pdist_name = pdist_name, dprimary_name = dprimary_name, ...
           )
         }
 
@@ -160,13 +171,15 @@ fitdistdoublecens <- function(censdata, distr,
 #' Define a fitdistrplus compatible wrapper around pprimarycensoreddist
 #' @inheritParams pprimarycensoreddist
 #' @keywords internal
-.ppcens <- function(q, pdist, pwindow, D, dprimary, dprimary_args, ...) {
+.ppcens <- function(q, pdist, pwindow, D, dprimary, dprimary_args,
+                    pdist_name, dprimary_name, ...) {
   tryCatch(
     {
       pprimarycensoreddist(
         q, pdist,
         pwindow = pwindow,
-        D = D, dprimary = dprimary, dprimary_args = dprimary_args, ...
+        D = D, dprimary = dprimary, dprimary_args = dprimary_args,
+        pdist_name = pdist_name, dprimary_name = dprimary_name, ...
       )
     },
     error = function(e) {
diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R
index 303cfb9..c1df072 100644
--- a/R/primary_censored_dist.R
+++ b/R/primary_censored_dist.R
@@ -64,21 +64,6 @@ primary_censored_cdf <- function(
 #' the numeric integration method.
 #'
 #' @inheritParams primary_censored_cdf
-#'
-#' @family primary_censored_dist
-#'
-#' @export
-primary_censored_cdf.default <- function(
-    object, q, pwindow, use_numeric = FALSE) {
-  primary_censored_cdf.pcens_numeric(object, q, pwindow, use_numeric)
-}
-
-#' Numeric method for computing primary event censored CDF
-#'
-#' This method uses numerical integration to compute the primary event censored
-#' CDF for any combination of delay distribution and primary event distribution.
-#'
-#' @inheritParams primary_censored_cdf
 #' @inheritParams pprimarycensoreddist
 #'
 #' @details
@@ -93,7 +78,7 @@ primary_censored_cdf.default <- function(
 #' @family primary_censored_dist
 #'
 #' @export
-primary_censored_cdf.pcens_numeric <- function(
+primary_censored_cdf.default <- function(
     object, q, pwindow, use_numeric = FALSE) {
   result <- vapply(q, function(d) {
     if (d <= 0) {
@@ -124,16 +109,63 @@ primary_censored_cdf.pcens_numeric <- function(
 #' @export
 primary_censored_cdf.pcens_pgamma_dunif <- function(
     object, q, pwindow, use_numeric = FALSE) {
-  use_numeric <- TRUE
   if (isTRUE(use_numeric)) {
     return(
-      primary_censored_cdf.pcens_numeric(object, q, pwindow, use_numeric)
+      primary_censored_cdf.default(object, q, pwindow, use_numeric)
     )
   }
+  # Extract Gamma distribution parameters
+  shape <- object$args$shape
+  scale <- object$args$scale
+  rate <- object$args$rate
+  # if we don't have scale get fromm rate
+  if (is.null(scale) && !is.null(rate)) {
+    scale <- 1 / rate
+  }
+  if (is.null(shape)) {
+    stop("shape parameter is required for Gamma distribution")
+  }
+  if (is.null(scale)) {
+    stop("scale or rate parameter is required for Gamma distribution")
+  }
 
-  result <- vapply(q, function(n) {
-    # Implement analytical solution here
-  }, numeric(1))
+  partial_pgamma <- function(q) {
+    pgamma(q, shape = shape, scale = scale)
+  }
+  partial_pgamm_k_1 <- function(q) {
+    pgamma(q, shape = shape + 1, scale = scale)
+  }
+  # Adjust q so that we have [q-pwindow, q]
+  q <- q - pwindow
+  # Handle cases where q + pwindow <= 0
+  zero_cases <- q + pwindow <= 0
+  result <- ifelse(zero_cases, 0, NA)
+
+  # Process non-zero cases only if there are any
+  if (!all(zero_cases)) {
+    non_zero_q <- q[!zero_cases]
+
+    # Compute necessary survival and distribution functions
+    pgamma_q <- partial_pgamma(non_zero_q)
+    pgamma_q_pwindow <- partial_pgamma(non_zero_q + pwindow)
+    pgamma_q_1 <- partial_pgamm_k_1(non_zero_q)
+    pgamma_q_pwindow_1 <- partial_pgamm_k_1(non_zero_q + pwindow)
+
+    Q_T <- 1 - pgamma_q_pwindow
+    Delta_F_T_kp1 <- pgamma_q_pwindow_1 - pgamma_q_1
+    Delta_F_T_k <- pgamma_q_pwindow - pgamma_q
+
+    # Calculate Q_Splus using the analytical formula
+    Q_Splus <- Q_T +
+      (shape * scale / pwindow) * Delta_F_T_kp1 -
+      (non_zero_q / pwindow) * Delta_F_T_k
+
+    # Compute the CDF as 1 - Q_Splus
+    non_zero_result <- 1 - Q_Splus
+
+    # Assign non-zero results back to the main result vector
+    result[!zero_cases] <- non_zero_result
+  }
 
   return(result)
 }
@@ -147,16 +179,62 @@ primary_censored_cdf.pcens_pgamma_dunif <- function(
 #' @export
 primary_censored_cdf.pcens_plnorm_dunif <- function(
     object, q, pwindow, use_numeric = FALSE) {
-  use_numeric <- TRUE
   if (isTRUE(use_numeric)) {
     return(
-      primary_censored_cdf.pcens_numeric(object, q, pwindow, use_numeric)
+      primary_censored_cdf.default(object, q, pwindow, use_numeric)
     )
   }
 
-  result <- vapply(q, function(n) {
-    # Implement analytical solution here
-  }, numeric(1))
+  # Extract Log-Normal distribution parameters
+  mu <- object$args$meanlog
+  sigma <- object$args$sdlog
+  if (is.null(mu)) {
+    stop("meanlog parameter is required for Log-Normal distribution")
+  }
+  if (is.null(sigma)) {
+    stop("sdlog parameter is required for Log-Normal distribution")
+  }
+
+  partial_plnorm <- function(q) {
+    plnorm(q, meanlog = mu, sdlog = sigma)
+  }
+  partial_plnorm_sigma2 <- function(q) {
+    plnorm(q, meanlog = mu + sigma^2, sdlog = sigma)
+  }
+  # Adjust q so that we have [q-pwindow, q]
+  q <- q - pwindow
+
+  # Handle cases where q + pwindow <= 0
+  zero_cases <- q + pwindow <= 0
+  result <- ifelse(zero_cases, 0, NA)
+
+  # Process non-zero cases only if there are any
+  if (!all(zero_cases)) {
+    non_zero_q <- q[!zero_cases]
+
+    # Compute necessary survival and distribution functions
+    plnorm_q <- partial_plnorm(non_zero_q)
+    plnorm_q_pwindow <- partial_plnorm(non_zero_q + pwindow)
+    plnorm_q_sigma2 <- partial_plnorm_sigma2(non_zero_q)
+    plnorm_q_pwindow_sigma2 <- partial_plnorm_sigma2(
+      non_zero_q + pwindow
+    )
+
+    Q_T <- 1 - plnorm_q_pwindow
+    Delta_F_T_mu_sigma <- plnorm_q_pwindow_sigma2 - plnorm_q_sigma2
+    Delta_F_T <- plnorm_q_pwindow - plnorm_q
+
+    # Calculate Q_Splus using the analytical formula
+    Q_Splus <- Q_T +
+      (exp(mu + 0.5 * sigma^2) / pwindow) * Delta_F_T_mu_sigma -
+      (non_zero_q / pwindow) * Delta_F_T
+
+    # Compute the CDF as 1 - Q_Splus
+    non_zero_result <- 1 - Q_Splus
+
+    # Assign non-zero results back to the main result vector
+    result[!zero_cases] <- non_zero_result
+  }
 
   return(result)
 }
diff --git a/man/dot-dpcens.Rd b/man/dot-dpcens.Rd
index 7d4fc60..135baaf 100644
--- a/man/dot-dpcens.Rd
+++ b/man/dot-dpcens.Rd
@@ -4,7 +4,18 @@
 \alias{.dpcens}
 \title{Define a fitdistrplus compatible wrapper around dprimarycensoreddist}
 \usage{
-.dpcens(x, swindows, pdist, pwindow, D, dprimary, dprimary_args, ...)
+.dpcens(
+  x,
+  swindows,
+  pdist,
+  pwindow,
+  D,
+  dprimary,
+  dprimary_args,
+  pdist_name,
+  dprimary_name,
+  ...
+)
 }
 \arguments{
 \item{x}{Vector of quantiles}
@@ -32,6 +43,18 @@ dprimary. For example, when using \code{dexpgrowth}, you would
 pass \code{list(min = 0, max = pwindow, r = 0.2)} to set the minimum, maximum,
 and rate parameters}
 
+\item{pdist_name}{A string specifying the name of the delay distribution
+function. If NULL, the function name is extracted using
+\code{\link[=.extract_function_name]{.extract_function_name()}}. Used to determine if a analytical solution
+exists for the primary censored distribution. Must be set if \code{pdist} is
+passed a pre-assigned variable rather than a function name.}
+
+\item{dprimary_name}{A string specifying the name of the primary event
+distribution function. If NULL, the function name is extracted using
+\code{\link[=.extract_function_name]{.extract_function_name()}}. Used to determine if a analytical solution
+exists for the primary censored distribution. Must be set if \code{dprimary} is
+passed a pre-assigned variable rather than a function name.}
+
 \item{...}{Additional arguments to be passed to the distribution function}
 }
 \description{
diff --git a/man/dot-ppcens.Rd b/man/dot-ppcens.Rd
index 3988afc..65f5428 100644
--- a/man/dot-ppcens.Rd
+++ b/man/dot-ppcens.Rd
@@ -4,7 +4,17 @@
 \alias{.ppcens}
 \title{Define a fitdistrplus compatible wrapper around pprimarycensoreddist}
 \usage{
-.ppcens(q, pdist, pwindow, D, dprimary, dprimary_args, ...)
+.ppcens(
+  q,
+  pdist,
+  pwindow,
+  D,
+  dprimary,
+  dprimary_args,
+  pdist_name,
+  dprimary_name,
+  ...
+)
 }
 \arguments{
 \item{q}{Vector of quantiles}
@@ -29,6 +39,18 @@ dprimary. For example, when using \code{dexpgrowth}, you would
 pass \code{list(min = 0, max = pwindow, r = 0.2)} to set the minimum, maximum,
 and rate parameters}
 
+\item{pdist_name}{A string specifying the name of the delay distribution
+function. If NULL, the function name is extracted using
+\code{\link[=.extract_function_name]{.extract_function_name()}}. Used to determine if a analytical solution
+exists for the primary censored distribution. Must be set if \code{pdist} is
+passed a pre-assigned variable rather than a function name.}
+
+\item{dprimary_name}{A string specifying the name of the primary event
+distribution function. If NULL, the function name is extracted using
+\code{\link[=.extract_function_name]{.extract_function_name()}}. Used to determine if a analytical solution
+exists for the primary censored distribution. Must be set if \code{dprimary} is
+passed a pre-assigned variable rather than a function name.}
+
 \item{...}{Additional arguments to be passed to pdist}
 }
 \description{
diff --git a/man/fitdistdoublecens.Rd b/man/fitdistdoublecens.Rd
index 7e3f63a..48dad8a 100644
--- a/man/fitdistdoublecens.Rd
+++ b/man/fitdistdoublecens.Rd
@@ -10,6 +10,7 @@ fitdistdoublecens(
   pwindow = 1,
   D = Inf,
   dprimary = stats::dunif,
+  dprimary_name = NULL,
   dprimary_args = list(),
   ...
 )
@@ -35,6 +36,12 @@ distribution over [0, pwindow]. Users can provide custom functions or use
 helper functions like \code{dexpgrowth} for an exponential growth distribution.
 See \code{primary_dists.R} for examples.}
 
+\item{dprimary_name}{A string specifying the name of the primary event
+distribution function. If NULL, the function name is extracted using
+\code{\link[=.extract_function_name]{.extract_function_name()}}. Used to determine if a analytical solution
+exists for the primary censored distribution. Must be set if \code{dprimary} is
+passed a pre-assigned variable rather than a function name.}
+
 \item{dprimary_args}{List of additional arguments to be passed to
 dprimary. For example, when using \code{dexpgrowth}, you would
 pass \code{list(min = 0, max = pwindow, r = 0.2)} to set the minimum, maximum,
diff --git a/man/new_primary_censored_dist.Rd b/man/new_primary_censored_dist.Rd
index 7c68b33..b207dab 100644
--- a/man/new_primary_censored_dist.Rd
+++ b/man/new_primary_censored_dist.Rd
@@ -53,7 +53,6 @@ S3 class for primary event censored distribution computation
 Low level primary event censored distribution objects and methods
 \code{\link{primary_censored_cdf}()},
 \code{\link{primary_censored_cdf.default}()},
-\code{\link{primary_censored_cdf.pcens_numeric}()},
 \code{\link{primary_censored_cdf.pcens_pgamma_dunif}()},
 \code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}
 }
diff --git a/man/primary_censored_cdf.Rd b/man/primary_censored_cdf.Rd
index 4dd2d7d..002f81a 100644
--- a/man/primary_censored_cdf.Rd
+++ b/man/primary_censored_cdf.Rd
@@ -29,7 +29,6 @@ Compute primary event censored CDF
 Low level primary event censored distribution objects and methods
 \code{\link{new_primary_censored_dist}()},
 \code{\link{primary_censored_cdf.default}()},
-\code{\link{primary_censored_cdf.pcens_numeric}()},
 \code{\link{primary_censored_cdf.pcens_pgamma_dunif}()},
 \code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}
 }
diff --git a/man/primary_censored_cdf.default.Rd b/man/primary_censored_cdf.default.Rd
index 2359b20..648c05e 100644
--- a/man/primary_censored_cdf.default.Rd
+++ b/man/primary_censored_cdf.default.Rd
@@ -24,11 +24,19 @@ This method serves as a fallback for combinations of delay and primary
 event distributions that don't have specific implementations. It uses
 the numeric integration method.
 }
+\details{
+This method implements the numerical integration approach for computing
+the primary event censored CDF. It uses the same mathematical formulation
+as described in the details section of \code{\link[=pprimarycensoreddist]{pprimarycensoreddist()}}, but
+applies numerical integration instead of analytical solutions.
+}
 \seealso{
+\code{\link[=pprimarycensoreddist]{pprimarycensoreddist()}} for the mathematical details of the
+primary event censored CDF computation.
+
 Low level primary event censored distribution objects and methods
 \code{\link{new_primary_censored_dist}()},
 \code{\link{primary_censored_cdf}()},
-\code{\link{primary_censored_cdf.pcens_numeric}()},
 \code{\link{primary_censored_cdf.pcens_pgamma_dunif}()},
 \code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}
 }
diff --git a/man/primary_censored_cdf.pcens_numeric.Rd b/man/primary_censored_cdf.pcens_numeric.Rd
deleted file mode 100644
index b43974f..0000000
--- a/man/primary_censored_cdf.pcens_numeric.Rd
+++ /dev/null
@@ -1,43 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/primary_censored_dist.R
-\name{primary_censored_cdf.pcens_numeric}
-\alias{primary_censored_cdf.pcens_numeric}
-\title{Numeric method for computing primary event censored CDF}
-\usage{
-\method{primary_censored_cdf}{pcens_numeric}(object, q, pwindow, use_numeric = FALSE)
-}
-\arguments{
-\item{object}{A \code{primary_censored_dist} object as created by
-\code{\link[=new_primary_censored_dist]{new_primary_censored_dist()}}.}
-
-\item{q}{Vector of quantiles}
-
-\item{pwindow}{Primary event window}
-
-\item{use_numeric}{Logical, if TRUE forces use of numeric integration
-even for distributions with analytical solutions. This is primarily
-useful for testing purposes or for settings where the analytical solution
-breaks down.}
-}
-\description{
-This method uses numerical integration to compute the primary event censored
-CDF for any combination of delay distribution and primary event distribution.
-}
-\details{
-This method implements the numerical integration approach for computing
-the primary event censored CDF. It uses the same mathematical formulation
-as described in the details section of \code{\link[=pprimarycensoreddist]{pprimarycensoreddist()}}, but
-applies numerical integration instead of analytical solutions.
-}
-\seealso{
-\code{\link[=pprimarycensoreddist]{pprimarycensoreddist()}} for the mathematical details of the
-primary event censored CDF computation.
-
-Low level primary event censored distribution objects and methods
-\code{\link{new_primary_censored_dist}()},
-\code{\link{primary_censored_cdf}()},
-\code{\link{primary_censored_cdf.default}()},
-\code{\link{primary_censored_cdf.pcens_pgamma_dunif}()},
-\code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}
-}
-\concept{primary_censored_dist}
diff --git a/man/primary_censored_cdf.pcens_pgamma_dunif.Rd b/man/primary_censored_cdf.pcens_pgamma_dunif.Rd
index 1ee894f..c333bd1 100644
--- a/man/primary_censored_cdf.pcens_pgamma_dunif.Rd
+++ b/man/primary_censored_cdf.pcens_pgamma_dunif.Rd
@@ -27,7 +27,6 @@ Low level primary event censored distribution objects and methods
 \code{\link{new_primary_censored_dist}()},
 \code{\link{primary_censored_cdf}()},
 \code{\link{primary_censored_cdf.default}()},
-\code{\link{primary_censored_cdf.pcens_numeric}()},
 \code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}
 }
 \concept{primary_censored_dist}
diff --git a/man/primary_censored_cdf.pcens_plnorm_dunif.Rd b/man/primary_censored_cdf.pcens_plnorm_dunif.Rd
index 387070f..e041be4 100644
--- a/man/primary_censored_cdf.pcens_plnorm_dunif.Rd
+++ b/man/primary_censored_cdf.pcens_plnorm_dunif.Rd
@@ -27,7 +27,6 @@ Low level primary event censored distribution objects and methods
 \code{\link{new_primary_censored_dist}()},
 \code{\link{primary_censored_cdf}()},
 \code{\link{primary_censored_cdf.default}()},
-\code{\link{primary_censored_cdf.pcens_numeric}()},
 \code{\link{primary_censored_cdf.pcens_pgamma_dunif}()}
 }
 \concept{primary_censored_dist}
diff --git a/tests/testthat/test-primary_censored_dist.R b/tests/testthat/test-primary_censored_dist.R
index 921a14b..fd6094d 100644
--- a/tests/testthat/test-primary_censored_dist.R
+++ b/tests/testthat/test-primary_censored_dist.R
@@ -25,120 +25,210 @@ test_that("new_primary_censored_dist creates object with correct structure", {
   expect_identical(obj, new_obj)
 })
 
-test_that("primary_censored_cdf.pcens_numeric computes correct values", {
-  pdist_name <- "pgamma"
-  pdist <- pgamma
-  dprimary_name <- "dunif"
-  dprimary <- dunif
-  shape <- 2
-  rate <- 1
-
-  obj <- new_primary_censored_dist(
-    pdist,
-    dprimary, list(),
-    pdist_name, dprimary_name,
-    shape = shape, rate = rate
-  )
-
-  q_values <- c(0, 5, 10, 15)
-  pwindow <- 10
-
-  result <- primary_censored_cdf(obj, q = q_values, pwindow = pwindow)
-
-  expect_true(all(result >= 0))
-  expect_true(all(result <= 1))
-  expect_type(result, "double")
-  expect_length(result, length(q_values))
-  expect_true(all(diff(result) >= 0)) # Ensure CDF is non-decreasing
-})
-
-test_that("primary_censored_cdf methods dispatch correctly", {
-  pdist_name <- "pgamma"
-  pdist <- pgamma
-  dprimary_name <- "dunif"
-  dprimary <- dunif
+test_that(
+  "primary_censored_cdf methods dispatch correctly to existing
+   analytical solutions",
+  {
+    pdist_name <- "pgamma"
+    pdist <- pgamma
+    dprimary_name <- "dunif"
+    dprimary <- dunif
+
+    obj_gamma <- new_primary_censored_dist(
+      pdist, dprimary, list(),
+      pdist_name, dprimary_name,
+      shape = 2, rate = 1
+    )
 
-  obj_gamma <- new_primary_censored_dist(
-    pdist, dprimary, list(),
-    pdist_name, dprimary_name,
-    shape = 2, rate = 1
-  )
+    pdist_name <- "plnorm"
+    pdist <- plnorm
+    dprimary_name <- "dunif"
+    dprimary <- dunif
 
-  pdist_name <- "plnorm"
-  pdist <- plnorm
-  dprimary_name <- "dunif"
-  dprimary <- dunif
+    obj_lnorm <- new_primary_censored_dist(
+      pdist, dprimary, list(),
+      pdist_name, dprimary_name,
+      meanlog = 0, sdlog = 1
+    )
 
-  obj_lnorm <- new_primary_censored_dist(
-    pdist, dprimary, list(),
-    pdist_name, dprimary_name,
-    meanlog = 0, sdlog = 1
-  )
+    expect_s3_class(obj_gamma, "pcens_pgamma_dunif")
+    expect_s3_class(obj_lnorm, "pcens_plnorm_dunif")
 
-  expect_s3_class(obj_gamma, "pcens_pgamma_dunif")
-  expect_s3_class(obj_lnorm, "pcens_plnorm_dunif")
+    q_values <- c(5, 10)
+    pwindow <- 10
 
-  q_values <- c(5, 10)
-  pwindow <- 10
+    expect_no_error(
+      primary_censored_cdf(obj_gamma, q = q_values, pwindow = pwindow)
+    )
+    expect_no_error(
+      primary_censored_cdf(obj_lnorm, q = q_values, pwindow = pwindow)
+    )
+  }
+)
+
+test_that(
+  "primary_censored_cdf errors as expected when the wrong distributional
+   parameters are supplied",
+  {
+    pdist_name <- "pgamma"
+    pdist <- pgamma
+    dprimary_name <- "dunif"
+    dprimary <- dunif
+
+    obj_gamma <- new_primary_censored_dist(
+      pdist, dprimary, list(),
+      pdist_name, dprimary_name,
+      rate = 1
+    )
 
-  expect_no_error(
-    primary_censored_cdf(obj_gamma, q = q_values, pwindow = pwindow)
-  )
-  expect_no_error(
-    primary_censored_cdf(obj_lnorm, q = q_values, pwindow = pwindow)
-  )
-})
+    expect_error(
+      primary_censored_cdf(obj_gamma, q = 1, pwindow = 1),
+      "shape parameter is required for Gamma distribution"
+    )
 
-test_that("primary_censored_cdf.pcens_pgamma_dunif uses numeric method", {
-  pdist_name <- "pgamma"
-  pdist <- pgamma
-  dprimary_name <- "dunif"
-  dprimary <- dunif
-  shape <- 2
-  rate <- 1
+    obj_gamma_no_rate <- new_primary_censored_dist(
+      pdist, dprimary, list(),
+      pdist_name, dprimary_name,
+      shape = 2
+    )
 
-  obj <- new_primary_censored_dist(
-    pdist,
-    dprimary, list(),
-    pdist_name, dprimary_name,
-    shape = shape, rate = rate
-  )
+    expect_error(
+      primary_censored_cdf(obj_gamma_no_rate, q = 1, pwindow = 1),
+      "scale or rate parameter is required for Gamma distribution"
+    )
 
-  q_values <- c(5, 10)
-  pwindow <- 10
+    pdist_name <- "plnorm"
+    pdist <- plnorm
 
-  expect_identical(
-    primary_censored_cdf(obj, q = q_values, pwindow = pwindow),
-    primary_censored_cdf(
-      obj,
-      q = q_values, pwindow = pwindow, use_numeric = TRUE
+    obj_lnorm_no_meanlog <- new_primary_censored_dist(
+      pdist, dprimary, list(),
+      pdist_name, dprimary_name,
+      sdlog = 1
     )
-  )
-})
 
-test_that("primary_censored_cdf.pcens_plnorm_dunif uses numeric method", {
-  pdist_name <- "plnorm"
-  pdist <- plnorm
-  dprimary_name <- "dunif"
-  dprimary <- dunif
-  meanlog <- 0
-  sdlog <- 1
-
-  obj <- new_primary_censored_dist(
-    pdist,
-    dprimary, list(),
-    pdist_name, dprimary_name,
-    meanlog = meanlog, sdlog = sdlog
-  )
+    expect_error(
+      primary_censored_cdf(obj_lnorm_no_meanlog, q = 1, pwindow = 1),
+      "meanlog parameter is required for Log-Normal distribution"
+    )
 
-  q_values <- c(5, 10)
-  pwindow <- 10
+    obj_lnorm_no_sdlog <- new_primary_censored_dist(
+      pdist, dprimary, list(),
+      pdist_name, dprimary_name,
+      meanlog = 0
+    )
 
-  expect_identical(
-    primary_censored_cdf(obj, q = q_values, pwindow = pwindow),
-    primary_censored_cdf(
-      obj,
-      q = q_values, pwindow = pwindow, use_numeric = TRUE
+    expect_error(
+      primary_censored_cdf(obj_lnorm_no_sdlog, q = 1, pwindow = 1),
+      "sdlog parameter is required for Log-Normal distribution"
     )
-  )
-})
+  }
+)
+
+test_that(
+  "primary_censored_cdf.default computes the same values as
+   primary_censored_cdf.pcens_pgamma_dunif",
+  {
+    pdist_name <- "pgamma"
+    pdist <- pgamma
+    dprimary_name <- "dunif"
+    dprimary <- dunif
+
+    shapes <- c(0.5, 1, 2, 5)
+    rates <- c(0.1, 0.5, 1, 2)
+    pwindows <- c(1, 2, 5, 10)
+
+    for (shape in shapes) {
+      for (rate in rates) {
+        for (pwindow in pwindows) {
+          obj <- new_primary_censored_dist(
+            pdist,
+            dprimary, list(),
+            pdist_name, dprimary_name,
+            shape = shape, rate = rate
+          )
+
+          q_values <- seq(0, 30, by = 0.1)
+          result_numeric <- primary_censored_cdf(
+            obj,
+            q = q_values, pwindow = pwindow, use_numeric = TRUE
+          )
+          result_analytical <- primary_censored_cdf(
+            obj,
+            q = q_values, pwindow = pwindow, use_numeric = FALSE
+          )
+
+          # Check properties of numeric result
+          expect_type(result_numeric, "double")
+          expect_length(result_numeric, length(q_values))
+          expect_true(
+            all(diff(result_numeric) >= 0)
+          ) # Ensure CDF is non-decreasing
+
+          # Check that analytical and numeric results are the same
+          expect_equal(
+            result_numeric, result_analytical,
+            tolerance = 1e-5,
+            info = sprintf(
+              "Mismatch for shape = %s, rate = %s, pwindow = %s",
+              shape, rate, pwindow
+            )
+          )
+        }
+      }
+    }
+  }
+)
+
+test_that(
+  "primary_censored_cdf.default computes the same values as
+   primary_censored_cdf.pcens_plnorm_dunif",
+  {
+    pdist_name <- "plnorm"
+    pdist <- plnorm
+    dprimary_name <- "dunif"
+    dprimary <- dunif
+
+    meanlogs <- c(-1, 0, 1, 2)
+    sdlogs <- c(0.5, 1, 1.5)
+    pwindows <- c(1, 2, 5, 8)
+
+    for (meanlog in meanlogs) {
+      for (sdlog in sdlogs) {
+        for (pwindow in pwindows) {
+          obj <- new_primary_censored_dist(
+            pdist,
+            dprimary, list(),
+            pdist_name, dprimary_name,
+            meanlog = meanlog, sdlog = sdlog
+          )
+
+          q_values <- seq(0, 30, by = 0.1)
+          result_numeric <- primary_censored_cdf(
+            obj,
+            q = q_values, pwindow = pwindow, use_numeric = TRUE
+          )
+          result_analytical <- primary_censored_cdf(
+            obj,
+            q = q_values, pwindow = pwindow, use_numeric = FALSE
+          )
+          # Check properties of numeric result
+          expect_type(result_numeric, "double")
+          expect_length(result_numeric, length(q_values))
+          expect_true(
+            all(diff(result_numeric) >= 0)
+          ) # Ensure CDF is non-decreasing
+
+          # Check that analytical and numeric results are the same
+          expect_equal(
+            result_numeric, result_analytical,
+            tolerance = 1e-5,
+            info = sprintf(
+              "Mismatch for meanlog = %s, sdlog = %s, pwindow = %s",
+              meanlog, sdlog, pwindow
+            )
+          )
+        }
+      }
+    }
+  }
+)
diff --git a/touchstone/script.R b/touchstone/script.R
index 1b7b637..74d8804 100644
--- a/touchstone/script.R
+++ b/touchstone/script.R
@@ -22,7 +22,7 @@ touchstone::benchmark_run(
     library(primarycensoreddist)
     q <- seq(0, 10, by = 0.01)
   },
-  pprimarycensoreddist_expgrowth = {
+  pprimarycensoreddist_plnrom_expgrowth = {
     pprimarycensoreddist(
       q, plnorm,
       dprimary = dexpgrowth,
@@ -51,7 +51,7 @@ touchstone::benchmark_run(
     library(primarycensoreddist)
     x <- seq(0, 10, by = 1)
   },
-  dprimarycensoreddist_expgrowth = {
+  dprimarycensoreddist_pweibull_expgrowth = {
     dprimarycensoreddist(
       x, pweibull,
       dprimary = dexpgrowth,
@@ -118,7 +118,7 @@ touchstone::benchmark_run(
       right = samples + swindow
     )
   },
-  fitdistdoublecens_gamma = {
+  fitdistdoublecens_gamma_expgrowth = {
     fitdistdoublecens(
       delay_data,
       distr = "gamma",