Issue 94: Add R and Stan for analytical Weibull

DESCRIPTION

@@ -35,6 +35,8 @@ URL: https://primarycensoreddist.epinowcast.org,
 BugReports: https://github.com/epinowcast/primarycensoreddist/issues/
 Depends:
     R (>= 4.0.0)
+Imports:
+    pracma
 Suggests:
     bookdown,
     cmdstanr,

NAMESPACE

@@ -3,6 +3,7 @@
 S3method(primary_censored_cdf,default)
 S3method(primary_censored_cdf,pcens_pgamma_dunif)
 S3method(primary_censored_cdf,pcens_plnorm_dunif)
+S3method(primary_censored_cdf,pcens_pweibull_dunif)
 export(check_dprimary)
 export(check_pdist)
 export(check_truncation)

NEWS.md

@@ -4,18 +4,19 @@ This is the development version of `primarycensoreddist` and is not yet ready fo
 
 ## Package
 
-* Split "Why it works" vignette into two separate vignettes, "Why it works" and "Analytic solutions for censored delay distributions".
 * Removed the need to assign functions to the global environment for `fitdistdoublecens()` by using `withr`.
 * Added a `check_truncation()` function to check if the truncation time is larger than the maximum observed delay. This is used in `fitdistdoublecens()` and `pcd_as_stan_data()` to ensure that the truncation time is appropriate to maximise computational efficiency.
 * `pcd_as_cmdstan_data()` has been renamed to `pcd_as_stan_data()` to better reflect that it is used for `Stan` models in general rather than just the `CmdStan` models.
 * The stan code has been refactored into a folder of functions within the current `stan` folder and the `stan` model has been moved into the `stan` folder. All paths to the stan code have been updated to reflect this.
+* Added R and stan implementations of the primary censored cdf for the weibull distribution with uniform primary censoring.
 
 ## Documentation
 
 * Simplified the "Analytic solutions" vignette by removing verbose derivation details.
 * Added links between vignettes to make it easier to navigate the documentation.
 * Added explicit usage of `pdist`, `dprimary`, `rdist`, and `rprimary` arguments in the getting started vignette to make it easier to link to mathematical details.
 * Fixed error in "Analytic solutions" vignette where the Weibull density was not being treated as zero for negative delays.
+* Split "Why it works" vignette into two separate vignettes, "Why it works" and "Analytic solutions for censored delay distributions".
 
 # primarycensoreddist 0.5.0
 

R/primary_censored_dist.R

@@ -247,3 +247,98 @@ 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)
+    )
+  }
+
+  if (pwindow > 3) {
+    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)
+  }
+
+  # Precompute constants
+  inv_shape <- 1 / shape
+  inv_scale <- 1 / scale
+
+  g <- function(t) {
+    # Use the lower incomplete gamma function
+    scaled_t <- (t * inv_scale)^shape
+    g_out <- vapply(scaled_t, function(x) {
+      a <- 1 + inv_shape
+      if (abs(-x + a * log(x)) > 700 || abs(a) > 170) {
+        return(0)
+      } else {
+        result <- pracma::gammainc(x, a)["lowinc"]
+      }
+      return(result)
+    }, numeric(1))
+    return(g_out)
+  }
+
+  # 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]
+    q_pwindow <- pmax(non_zero_q + pwindow, 0)
+    non_zero_q_pos <- pmax(non_zero_q, 0)
+
+    # Compute necessary survival and distribution functions
+    pweibull_q <- partial_pweibull(non_zero_q_pos)
+    pweibull_q_pwindow <- partial_pweibull(q_pwindow)
+    g_q <- g(non_zero_q_pos)
+    g_q_pwindow <- g(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)
+}

inst/stan/functions/primary_censored_dist_analytical_cdf.stan

@@ -10,6 +10,7 @@
 int check_for_analytical(int dist_id, int primary_dist_id) {
   if (dist_id == 2 && primary_dist_id == 1) return 1; // Gamma delay with Uniform primary
   if (dist_id == 1 && primary_dist_id == 1) return 1; // Lognormal delay with Uniform primary
+  if (dist_id == 3 && primary_dist_id == 1) return 1; // Weibull delay with Uniform primary
   return 0; // No analytical solution for other combinations
 }
 
@@ -125,6 +126,79 @@ real primary_censored_lognormal_uniform_lcdf(data real d, real q, array[] real p
   return log_F_Splus;
 }
 
+/**
+  * Compute the log of the lower incomplete gamma function
+  *
+  * This function is used in the analytical solution for the primary censored
+  * Weibull distribution with uniform primary censoring. It corresponds to the
+  * g(t; λ, k) function described in the analytic solutions document.
+  *
+  * @param t Upper bound of integration
+  * @param shape Shape parameter (k) of the Weibull distribution
+  * @param scale Scale parameter (λ) of the Weibull distribution
+  *
+  * @return Log of g(t; λ, k) = γ(1 + 1/k, (t/λ)^k)
+  */
+real log_weibull_g(real t, real shape, real scale) {
+  real x = pow(t * inv(scale), shape);
+  real a = 1 + inv(shape);
+  return log(gamma_p(a, x)) + lgamma(a);
+}
+
+/**
+  * Compute the primary event censored log CDF analytically for Weibull delay with Uniform primary
+  *
+  * @param d Delay time
+  * @param q Lower bound of integration (max(d - pwindow, 0))
+  * @param params Array of Weibull distribution parameters [shape, scale]
+  * @param pwindow Primary event window
+  *
+  * @return Log of the primary event censored CDF for Weibull delay with
+  * Uniform primary
+  */
+real primary_censored_weibull_uniform_lcdf(data real d, real q, array[] real params, data real pwindow) {
+  real shape = params[1];
+  real scale = params[2];
+  real log_window = log(pwindow);
+
+  real log_F_T = weibull_lcdf(d | shape, scale);
+
+  real log_delta_g;
+  real log_delta_F_T;
+  real log_F_Splus;
+
+  if (q != 0) {
+    real log_F_T_q = weibull_lcdf(q | shape, scale);
+
+    log_delta_g = log_diff_exp(
+      log_weibull_g(d, shape, scale),
+      log_weibull_g(q, shape, scale)
+    );
+    log_delta_F_T = log_diff_exp(log_F_T, log_F_T_q);
+
+    log_F_Splus = log_diff_exp(
+      log_F_T,
+      log_diff_exp(
+        log(scale) + log_delta_g,
+        log(d - pwindow) + log_delta_F_T
+      ) - log_window
+    );
+  } else {
+    log_delta_g = log_weibull_g(d, shape, scale);
+    log_delta_F_T = log_F_T;
+
+    log_F_Splus = log_diff_exp(
+      log_F_T,
+      log_sum_exp(
+        log(scale) + log_delta_g,
+        log(pwindow - d) + log_delta_F_T
+      ) - log_window
+    );
+  }
+
+  return log_F_Splus;
+}
+
 /**
   * Compute the primary event censored log CDF analytically for a single delay
   *
@@ -157,6 +231,9 @@ real primary_censored_dist_analytical_lcdf(data real d, int dist_id,
   } else if (dist_id == 1 && primary_dist_id == 1) {
     // Lognormal delay with Uniform primary
     result = primary_censored_lognormal_uniform_lcdf(d | q, params, pwindow);
+  } else if (dist_id == 3 && primary_dist_id == 1) {
+    // Weibull delay with Uniform primary
+    result = primary_censored_weibull_uniform_lcdf(d | q, params, pwindow);
   } else {
     // No analytical solution available
     return negative_infinity();

tests/testthat/test-fitdistdoublecens.R

@@ -139,19 +139,59 @@ test_that("fitdistdoublecens works with mixed secondary windows", {
   expect_equal(unname(fit_gamma$estimate["rate"]), true_rate, tolerance = 0.2)
 })
 
-test_that("fitdistdoublecens throws error when fitdistrplus is not installed", {
-  with_mocked_bindings(
-    {
-      # Create dummy data
-      dummy_data <- data.frame(left = 1:5, right = 2:6)
-
-      # Expect an error when trying to use fitdistdoublecens
-      expect_error(
-        fitdistdoublecens(dummy_data, "norm"),
-        "Package 'fitdistrplus' is required but not installed for this"
-      )
-    },
-    requireNamespace = function(...) FALSE,
-    .package = "base"
-  )
-})
+test_that(
+  "fitdistdoublecens throws error when required packages are not installed",
+  {
+    # Create dummy data
+    dummy_data <- data.frame(left = 1:5, right = 2:6)
+
+    # Test for fitdistrplus
+    with_mocked_bindings(
+      {
+        expect_error(
+          fitdistdoublecens(dummy_data, "norm"),
+          "Package 'fitdistrplus' is required but not installed for this",
+          fixed = TRUE
+        )
+      },
+      requireNamespace = function(pkg, ...) {
+        if (pkg == "fitdistrplus") {
+          return(FALSE)
+        }
+        TRUE
+      },
+      .package = "base"
+    )
+
+    # Test for withr
+    with_mocked_bindings(
+      {
+        expect_error(
+          fitdistdoublecens(dummy_data, "norm"),
+          "Package 'withr' is required but not installed for this function.",
+          fixed = TRUE
+        )
+      },
+      requireNamespace = function(pkg, ...) {
+        if (pkg == "withr") {
+          return(FALSE)
+        }
+        TRUE
+      },
+      .package = "base"
+    )
+
+    # Test when both packages are missing
+    with_mocked_bindings(
+      {
+        expect_error(
+          fitdistdoublecens(dummy_data, "norm"),
+          "Package 'fitdistrplus' is required but not installed",
+          fixed = TRUE
+        )
+      },
+      requireNamespace = function(...) FALSE,
+      .package = "base"
+    )
+  }
+)