add R weibull solution

R/primary_censored_dist.R

@@ -247,3 +247,80 @@ primary_censored_cdf.pcens_plnorm_dunif <- function(
 
   return(result)
 }
+
+#' Method for Weibull delay with uniform primary
+#'
+#' @inheritParams primary_censored_cdf
+#'
+#' @family primary_censored_dist
+#'
+#' @export
+primary_censored_cdf.pcens_pweibull_dunif <- function(
+    object, q, pwindow, use_numeric = FALSE) {
+  if (isTRUE(use_numeric)) {
+    return(
+      primary_censored_cdf.default(object, q, pwindow, use_numeric)
+    )
+  }
+
+  # Extract Weibull distribution parameters
+  shape <- object$args$shape
+  scale <- object$args$scale
+  if (is.null(shape)) {
+    stop("shape parameter is required for Weibull distribution")
+  }
+  if (is.null(scale)) {
+    stop("scale parameter is required for Weibull distribution")
+  }
+
+  partial_pweibull <- function(q) {
+    stats::pweibull(q, shape = shape, scale = scale)
+  }
+
+  g <- function(t) {
+    gamma_1k <- pgamma(
+      (t / scale)^shape,
+      shape = 1 / shape,
+      lower.tail = TRUE
+    ) / shape
+    gamma_1k - (t / scale) * exp(-(t / scale)^shape)
+  }
+
+  # 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
+    pweibull_q <- partial_pweibull(non_zero_q)
+    pweibull_q_pwindow <- partial_pweibull(non_zero_q + pwindow)
+    g_q <- g(non_zero_q)
+    g_q_pwindow <- g(non_zero_q + pwindow)
+
+    Q_T <- 1 - pweibull_q_pwindow
+    Delta_g <- g_q_pwindow - g_q
+    Delta_F_T <- pweibull_q_pwindow - pweibull_q
+
+    # Calculate Q_Splus using the analytical formula
+    Q_Splus <- Q_T +
+      (scale / pwindow) * Delta_g -
+      (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
+  }
+
+  # Ensure the result is not greater than 1 (accounts for numerical errors)
+  result <- pmin(1, result)
+
+  return(result)
+}

tests/testthat/test-primary_censored_dist.R

@@ -24,7 +24,6 @@ test_that("new_primary_censored_dist creates object with correct structure", {
   )
   expect_identical(obj, new_obj)
 })
-
 test_that(
   "primary_censored_cdf methods dispatch correctly to existing
    analytical solutions",
@@ -51,8 +50,20 @@ test_that(
       meanlog = 0, sdlog = 1
     )
 
+    pdist_name <- "pweibull"
+    pdist <- pweibull
+    dprimary_name <- "dunif"
+    dprimary <- dunif
+
+    obj_weibull <- new_primary_censored_dist(
+      pdist, dprimary, list(),
+      pdist_name, dprimary_name,
+      shape = 2, scale = 1
+    )
+
     expect_s3_class(obj_gamma, "pcens_pgamma_dunif")
     expect_s3_class(obj_lnorm, "pcens_plnorm_dunif")
+    expect_s3_class(obj_weibull, "pcens_pweibull_dunif")
 
     q_values <- c(5, 10)
     pwindow <- 10
@@ -63,6 +74,9 @@ test_that(
     expect_no_error(
       primary_censored_cdf(obj_lnorm, q = q_values, pwindow = pwindow)
     )
+    expect_no_error(
+      primary_censored_cdf(obj_weibull, q = q_values, pwindow = pwindow)
+    )
   }
 )
 
@@ -232,3 +246,58 @@ test_that(
     }
   }
 )
+
+test_that(
+  "primary_censored_cdf.default computes the same values as
+   primary_censored_cdf.pcens_pweibull_dunif",
+  {
+    pdist_name <- "pweibull"
+    pdist <- pweibull
+    dprimary_name <- "dunif"
+    dprimary <- dunif
+
+    shapes <- c(0.5, 1, 2, 3)
+    scales <- c(0.5, 1, 2, 5)
+    pwindows <- c(1, 2, 5, 10)
+
+    for (shape in shapes) {
+      for (scale in scales) {
+        for (pwindow in pwindows) {
+          obj <- new_primary_censored_dist(
+            pdist,
+            dprimary, list(),
+            pdist_name, dprimary_name,
+            shape = shape, scale = scale
+          )
+
+          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, scale = %s, pwindow = %s",
+              shape, scale, pwindow
+            )
+          )
+        }
+      }
+    }
+  }
+)