diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R index 419b33d..76bfed7 100644 --- a/R/primary_censored_dist.R +++ b/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) +} diff --git a/tests/testthat/test-primary_censored_dist.R b/tests/testthat/test-primary_censored_dist.R index fd6094d..4df96b0 100644 --- a/tests/testthat/test-primary_censored_dist.R +++ b/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 + ) + ) + } + } + } + } +)