diff --git a/DESCRIPTION b/DESCRIPTION index 45d77da..16eb8a2 100644 --- a/DESCRIPTION +++ b/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, diff --git a/NAMESPACE b/NAMESPACE index c19ecc3..96309f2 100644 --- a/NAMESPACE +++ b/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) diff --git a/NEWS.md b/NEWS.md index 49c6e7a..08e6fe0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,11 +4,11 @@ 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 @@ -16,6 +16,7 @@ This is the development version of `primarycensoreddist` and is not yet ready fo * 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 diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R index 419b33d..b088895 100644 --- a/R/primary_censored_dist.R +++ b/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) +} diff --git a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan b/inst/stan/functions/primary_censored_dist_analytical_cdf.stan index 44a896a..ae96298 100644 --- a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan +++ b/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(); diff --git a/man/new_primary_censored_dist.Rd b/man/new_primary_censored_dist.Rd index b207dab..8413302 100644 --- a/man/new_primary_censored_dist.Rd +++ b/man/new_primary_censored_dist.Rd @@ -54,6 +54,7 @@ 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_pgamma_dunif}()}, -\code{\link{primary_censored_cdf.pcens_plnorm_dunif}()} +\code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}, +\code{\link{primary_censored_cdf.pcens_pweibull_dunif}()} } \concept{primary_censored_dist} diff --git a/man/primary_censored_cdf.Rd b/man/primary_censored_cdf.Rd index 002f81a..d425fd5 100644 --- a/man/primary_censored_cdf.Rd +++ b/man/primary_censored_cdf.Rd @@ -30,6 +30,7 @@ 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_pgamma_dunif}()}, -\code{\link{primary_censored_cdf.pcens_plnorm_dunif}()} +\code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}, +\code{\link{primary_censored_cdf.pcens_pweibull_dunif}()} } \concept{primary_censored_dist} diff --git a/man/primary_censored_cdf.default.Rd b/man/primary_censored_cdf.default.Rd index 648c05e..6efe5ff 100644 --- a/man/primary_censored_cdf.default.Rd +++ b/man/primary_censored_cdf.default.Rd @@ -38,6 +38,7 @@ 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_pgamma_dunif}()}, -\code{\link{primary_censored_cdf.pcens_plnorm_dunif}()} +\code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}, +\code{\link{primary_censored_cdf.pcens_pweibull_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 c333bd1..b88845a 100644 --- a/man/primary_censored_cdf.pcens_pgamma_dunif.Rd +++ b/man/primary_censored_cdf.pcens_pgamma_dunif.Rd @@ -27,6 +27,7 @@ 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_plnorm_dunif}()} +\code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}, +\code{\link{primary_censored_cdf.pcens_pweibull_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 e041be4..2059694 100644 --- a/man/primary_censored_cdf.pcens_plnorm_dunif.Rd +++ b/man/primary_censored_cdf.pcens_plnorm_dunif.Rd @@ -27,6 +27,7 @@ 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_pgamma_dunif}()}, +\code{\link{primary_censored_cdf.pcens_pweibull_dunif}()} } \concept{primary_censored_dist} diff --git a/man/primary_censored_cdf.pcens_pweibull_dunif.Rd b/man/primary_censored_cdf.pcens_pweibull_dunif.Rd new file mode 100644 index 0000000..3e1839e --- /dev/null +++ b/man/primary_censored_cdf.pcens_pweibull_dunif.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/primary_censored_dist.R +\name{primary_censored_cdf.pcens_pweibull_dunif} +\alias{primary_censored_cdf.pcens_pweibull_dunif} +\title{Method for Weibull delay with uniform primary} +\usage{ +\method{primary_censored_cdf}{pcens_pweibull_dunif}(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{ +Method for Weibull delay with uniform primary +} +\seealso{ +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/tests/testthat/test-fitdistdoublecens.R b/tests/testthat/test-fitdistdoublecens.R index 52762f3..d2a3abb 100644 --- a/tests/testthat/test-fitdistdoublecens.R +++ b/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" + ) + } +) diff --git a/tests/testthat/test-primary_censored_dist.R b/tests/testthat/test-primary_censored_dist.R index fd6094d..5a0bdf7 100644 --- a/tests/testthat/test-primary_censored_dist.R +++ b/tests/testthat/test-primary_censored_dist.R @@ -51,11 +51,23 @@ 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 + pwindow <- 2 expect_no_error( primary_censored_cdf(obj_gamma, q = q_values, pwindow = pwindow) @@ -63,6 +75,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) + ) } ) @@ -121,6 +136,31 @@ test_that( primary_censored_cdf(obj_lnorm_no_sdlog, q = 1, pwindow = 1), "sdlog parameter is required for Log-Normal distribution" ) + + pdist_name <- "pweibull" + pdist <- pweibull + + obj_weibull_no_shape <- new_primary_censored_dist( + pdist, dprimary, list(), + pdist_name, dprimary_name, + scale = 1 + ) + + expect_error( + primary_censored_cdf(obj_weibull_no_shape, q = 1, pwindow = 1), + "shape parameter is required for Weibull distribution" + ) + + obj_weibull_no_scale <- new_primary_censored_dist( + pdist, dprimary, list(), + pdist_name, dprimary_name, + shape = 2 + ) + + expect_error( + primary_censored_cdf(obj_weibull_no_scale, q = 1, pwindow = 1), + "scale parameter is required for Weibull distribution" + ) } ) @@ -232,3 +272,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) + scales <- c(0.5, 1, 2) + pwindows <- c(1, 2, 3, 4, 5) + + 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 + ) + ) + } + } + } + } +) diff --git a/tests/testthat/test-rpd-primarycensoreddist.R b/tests/testthat/test-rpd-primarycensoreddist.R index 4d3154a..b85c97d 100644 --- a/tests/testthat/test-rpd-primarycensoreddist.R +++ b/tests/testthat/test-rpd-primarycensoreddist.R @@ -104,6 +104,8 @@ test_that( samples <- rpcens( n, rweibull, pwindow, swindow, + rprimary = rexpgrowth, + rprimary_args = list(r = 0.5), D = D, shape = shape, scale = scale ) @@ -115,6 +117,8 @@ test_that( x_values <- seq(0, D - swindow, by = swindow) pmf <- dpcens( x_values, pweibull, pwindow, swindow, + dprimary = dexpgrowth, + dprimary_args = list(r = 0.5), D = D, shape = shape, scale = scale ) theoretical_mean <- sum(x_values * pmf) @@ -132,6 +136,8 @@ test_that( empirical_cdf <- ecdf(samples)(x_values) theoretical_cdf <- ppcens( c(x_values[-1], D), pweibull, pwindow, D, + dprimary = dexpgrowth, + dprimary_args = list(r = 0.5), shape = shape, scale = scale ) expect_equal(cumsum(pmf), theoretical_cdf, tolerance = 0.01) diff --git a/tests/testthat/test-stan-primary_censored_dist_analytical_cdf.R b/tests/testthat/test-stan-primary_censored_dist_analytical_cdf.R new file mode 100644 index 0000000..7698746 --- /dev/null +++ b/tests/testthat/test-stan-primary_censored_dist_analytical_cdf.R @@ -0,0 +1,140 @@ +skip_on_cran() +if (on_ci()) { + skip_on_os("windows") + skip_on_os("mac") +} + +test_that( + "Stan primary_censored_dist_analytical_lcdf matches R implementation for + Gamma", + { + 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( + pgamma, + dunif, list(), + "pgamma", "dunif", + shape = shape, rate = rate + ) + + q_values <- seq(0, 30, by = 1) + r_result <- primary_censored_cdf( + obj, + q = q_values, pwindow = pwindow, use_numeric = FALSE + ) + + stan_result <- vapply(q_values, function(q) { + primary_censored_dist_analytical_cdf( + q, 2, c(shape, rate), pwindow, Inf, 1, numeric(0) + ) + }, numeric(1)) + + stan_result <- ifelse(is.nan(stan_result), 1, stan_result) + expect_equal( + r_result, stan_result, + tolerance = 1e-6, + info = sprintf( + "Mismatch for shape = %s, rate = %s, pwindow = %s", + shape, rate, pwindow + ) + ) + } + } + } + } +) + +test_that( + "Stan primary_censored_dist_analytical_lcdf matches R implementation for + Lognormal", + { + 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( + plnorm, + dunif, list(), + "plnorm", "dunif", + meanlog = meanlog, sdlog = sdlog + ) + + q_values <- seq(0, 30, by = 1) + r_result <- primary_censored_cdf( + obj, + q = q_values, pwindow = pwindow, use_numeric = FALSE + ) + + stan_result <- vapply(q_values, function(q) { + primary_censored_dist_analytical_cdf( + q, 1, c(meanlog, sdlog), pwindow, Inf, 1, numeric(0) + ) + }, numeric(1)) + + stan_result <- ifelse(is.nan(stan_result), 1, stan_result) + expect_equal( + r_result, stan_result, + tolerance = 1e-6, + info = sprintf( + "Mismatch for meanlog = %s, sdlog = %s, pwindow = %s", + meanlog, sdlog, pwindow + ) + ) + } + } + } + } +) + +test_that( + "Stan primary_censored_dist_analytical_lcdf matches R implementation for + Weibull", + { + 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( + pweibull, + dunif, list(), + "pweibull", "dunif", + shape = shape, scale = scale + ) + + q_values <- seq(0, 30, by = 1) + r_result <- primary_censored_cdf( + obj, + q = q_values, pwindow = pwindow, use_numeric = FALSE + ) + + stan_result <- vapply(q_values, function(q) { + primary_censored_dist_analytical_cdf( + q, 3, c(shape, scale), pwindow, Inf, 1, numeric(0) + ) + }, numeric(1)) + + stan_result <- ifelse(is.nan(stan_result), 1, stan_result) + expect_equal( + r_result, stan_result, + tolerance = 1e-6, + info = sprintf( + "Mismatch for shape = %s, scale = %s, pwindow = %s", + shape, scale, pwindow + ) + ) + } + } + } + } +)