From 1a456b60907c1776f375b5152b18ec0e8896e637 Mon Sep 17 00:00:00 2001 From: Sam Date: Fri, 27 Sep 2024 18:48:46 +0100 Subject: [PATCH 01/24] add R weibull solution --- R/primary_censored_dist.R | 77 +++++++++++++++++++++ tests/testthat/test-primary_censored_dist.R | 71 ++++++++++++++++++- 2 files changed, 147 insertions(+), 1 deletion(-) 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 + ) + ) + } + } + } + } +) From e91da869288d1093fd01c103c7e6a9d30dc6c49a Mon Sep 17 00:00:00 2001 From: Sam Date: Fri, 27 Sep 2024 19:58:30 +0100 Subject: [PATCH 02/24] update docs and news --- NAMESPACE | 1 + NEWS.md | 4 +++ man/new_primary_censored_dist.Rd | 3 +- man/primary_censored_cdf.Rd | 3 +- man/primary_censored_cdf.default.Rd | 3 +- ...primary_censored_cdf.pcens_pgamma_dunif.Rd | 3 +- ...primary_censored_cdf.pcens_plnorm_dunif.Rd | 3 +- ...imary_censored_cdf.pcens_pweibull_dunif.Rd | 33 +++++++++++++++++++ 8 files changed, 48 insertions(+), 5 deletions(-) create mode 100644 man/primary_censored_cdf.pcens_pweibull_dunif.Rd 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..c8aa81d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,10 @@ This is the development version of `primarycensoreddist` and is not yet ready fo ## Package +* Added R and stan implementations of the primary censored cdf for the weibull distribution with uniform primary censoring. + +## Documentation + * 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. 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} From 196c45bab12eb129af8d4beadffb77787f9ca02a Mon Sep 17 00:00:00 2001 From: Sam Date: Fri, 27 Sep 2024 20:22:46 +0100 Subject: [PATCH 03/24] add stan implementation --- .../primary_censored_dist_analytical_cdf.stan | 72 +++++++++++++++++++ 1 file changed, 72 insertions(+) diff --git a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan b/inst/stan/functions/primary_censored_dist_analytical_cdf.stan index 44a896a..4a1ee5c 100644 --- a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan +++ b/inst/stan/functions/primary_censored_dist_analytical_cdf.stan @@ -125,6 +125,78 @@ real primary_censored_lognormal_uniform_lcdf(data real d, real q, array[] real p return log_F_Splus; } +/** + * Compute the log of the integral of the Weibull PDF from 0 to t + * + * 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) / k - (t/λ) * exp(-(t/λ)^k) + */ +real log_weibull_g(real t, real shape, real scale) { + real x = (t / scale)^shape; + return log(gamma_p(1.0 / shape, x) / shape) - log(scale) - x; +} + +/** + * 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 * From 7276585efaaddfb9cdf04180771fc02fbe90b695 Mon Sep 17 00:00:00 2001 From: Sam Date: Fri, 27 Sep 2024 21:09:59 +0100 Subject: [PATCH 04/24] do some target stan tests --- ...tan-primary_censored_dist_analytical_cdf.R | 108 ++++++++++++++++++ 1 file changed, 108 insertions(+) create mode 100644 tests/testthat/test-stan-primary_censored_dist_analytical_cdf.R 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..8993d86 --- /dev/null +++ b/tests/testthat/test-stan-primary_censored_dist_analytical_cdf.R @@ -0,0 +1,108 @@ +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)) + expect_equal(r_result, stan_result, tolerance = 1e-6) + } + } + } + } +) + +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)) + expect_equal(r_result, stan_result, tolerance = 1e-6) + } + } + } + } +) + +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)) + + expect_equal(r_result, stan_result, tolerance = 1e-6) + } + } + } + } +) From 45cb93117adc8a6ed07d3511798ee61089d5963f Mon Sep 17 00:00:00 2001 From: Sam Date: Mon, 7 Oct 2024 14:22:37 +0100 Subject: [PATCH 05/24] revert skipped linebreak --- tests/testthat/test-primary_censored_dist.R | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/testthat/test-primary_censored_dist.R b/tests/testthat/test-primary_censored_dist.R index 4df96b0..d97d258 100644 --- a/tests/testthat/test-primary_censored_dist.R +++ b/tests/testthat/test-primary_censored_dist.R @@ -24,6 +24,7 @@ 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", From 691235defcd86a0f967ee229b608ad86d7de0ded Mon Sep 17 00:00:00 2001 From: Sam Date: Tue, 8 Oct 2024 10:39:57 +0100 Subject: [PATCH 06/24] push refactored g --- DESCRIPTION | 2 ++ R/primary_censored_dist.R | 18 +++++++++++------- .../primary_censored_dist_analytical_cdf.stan | 7 +++++-- 3 files changed, 18 insertions(+), 9 deletions(-) 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/R/primary_censored_dist.R b/R/primary_censored_dist.R index 76bfed7..acb2c91 100644 --- a/R/primary_censored_dist.R +++ b/R/primary_censored_dist.R @@ -277,13 +277,17 @@ primary_censored_cdf.pcens_pweibull_dunif <- function( stats::pweibull(q, shape = shape, scale = scale) } + # Precompute constants + inv_shape <- 1 / shape + inv_scale <- 1 / 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) + # Use the lower incomplete gamma function + scaled_t <- (t * inv_scale)^shape + gamma_1k <- vapply(scaled_t, function(x) { + pracma::gammainc(inv_shape, x)["lowinc"] + }, numeric(1)) * inv_shape + scale * gamma_1k - t * exp(-scaled_t) } # Adjust q so that we have [q-pwindow, q] @@ -309,7 +313,7 @@ primary_censored_cdf.pcens_pweibull_dunif <- function( # Calculate Q_Splus using the analytical formula Q_Splus <- Q_T + - (scale / pwindow) * Delta_g - + (1 / pwindow) * Delta_g - (non_zero_q / pwindow) * Delta_F_T # Compute the CDF as 1 - Q_Splus diff --git a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan b/inst/stan/functions/primary_censored_dist_analytical_cdf.stan index 4a1ee5c..20ff130 100644 --- a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan +++ b/inst/stan/functions/primary_censored_dist_analytical_cdf.stan @@ -139,8 +139,11 @@ real primary_censored_lognormal_uniform_lcdf(data real d, real q, array[] real p * @return Log of g(t; λ, k) = γ(1 + 1/k, (t/λ)^k) / k - (t/λ) * exp(-(t/λ)^k) */ real log_weibull_g(real t, real shape, real scale) { - real x = (t / scale)^shape; - return log(gamma_p(1.0 / shape, x) / shape) - log(scale) - x; + real x = pow(t / scale, shape); + real log_gamma_term = log(gamma_p(1 + inv(shape), x)) - log(shape); + real log_exp_term = log(t) - log(scale) - x; + + return log_diff_exp(log_gamma_term, log_exp_term); } /** From ed6fe892079905f22bc95320acf97abf31e75ada Mon Sep 17 00:00:00 2001 From: Sam Date: Tue, 8 Oct 2024 10:44:46 +0100 Subject: [PATCH 07/24] correct reversion --- R/primary_censored_dist.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R index acb2c91..ea24493 100644 --- a/R/primary_censored_dist.R +++ b/R/primary_censored_dist.R @@ -285,7 +285,7 @@ primary_censored_cdf.pcens_pweibull_dunif <- function( # Use the lower incomplete gamma function scaled_t <- (t * inv_scale)^shape gamma_1k <- vapply(scaled_t, function(x) { - pracma::gammainc(inv_shape, x)["lowinc"] + pracma::gammainc(1 + inv_shape, x)["lowinc"] }, numeric(1)) * inv_shape scale * gamma_1k - t * exp(-scaled_t) } From a6e36eebd151ea9f54fc2287d31a87eb1866c528 Mon Sep 17 00:00:00 2001 From: Sam Date: Tue, 8 Oct 2024 10:54:06 +0100 Subject: [PATCH 08/24] rephrase --- R/primary_censored_dist.R | 5 ++++- .../stan/functions/primary_censored_dist_analytical_cdf.stan | 4 +--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R index ea24493..e26e2f8 100644 --- a/R/primary_censored_dist.R +++ b/R/primary_censored_dist.R @@ -285,9 +285,12 @@ primary_censored_cdf.pcens_pweibull_dunif <- function( # Use the lower incomplete gamma function scaled_t <- (t * inv_scale)^shape gamma_1k <- vapply(scaled_t, function(x) { + if (x <= 0) { + return(0) + } pracma::gammainc(1 + inv_shape, x)["lowinc"] }, numeric(1)) * inv_shape - scale * gamma_1k - t * exp(-scaled_t) + scale * gamma_1k } # Adjust q so that we have [q-pwindow, q] diff --git a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan b/inst/stan/functions/primary_censored_dist_analytical_cdf.stan index 20ff130..28041e7 100644 --- a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan +++ b/inst/stan/functions/primary_censored_dist_analytical_cdf.stan @@ -141,9 +141,7 @@ real primary_censored_lognormal_uniform_lcdf(data real d, real q, array[] real p real log_weibull_g(real t, real shape, real scale) { real x = pow(t / scale, shape); real log_gamma_term = log(gamma_p(1 + inv(shape), x)) - log(shape); - real log_exp_term = log(t) - log(scale) - x; - - return log_diff_exp(log_gamma_term, log_exp_term); + return log_gamma_term; } /** From f7c8a0d68ebf3bafdc5bab19fdb5e53787b839af Mon Sep 17 00:00:00 2001 From: Sam Date: Tue, 8 Oct 2024 11:09:20 +0100 Subject: [PATCH 09/24] update with @samuelbrand1 don't get confused feedback --- R/primary_censored_dist.R | 9 +++------ .../functions/primary_censored_dist_analytical_cdf.stan | 5 ++--- 2 files changed, 5 insertions(+), 9 deletions(-) diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R index e26e2f8..a95c492 100644 --- a/R/primary_censored_dist.R +++ b/R/primary_censored_dist.R @@ -285,12 +285,9 @@ primary_censored_cdf.pcens_pweibull_dunif <- function( # Use the lower incomplete gamma function scaled_t <- (t * inv_scale)^shape gamma_1k <- vapply(scaled_t, function(x) { - if (x <= 0) { - return(0) - } pracma::gammainc(1 + inv_shape, x)["lowinc"] - }, numeric(1)) * inv_shape - scale * gamma_1k + }, numeric(1)) + return(gamma_1k) } # Adjust q so that we have [q-pwindow, q] @@ -311,7 +308,7 @@ primary_censored_cdf.pcens_pweibull_dunif <- function( g_q_pwindow <- g(non_zero_q + pwindow) Q_T <- 1 - pweibull_q_pwindow - Delta_g <- g_q_pwindow - g_q + Delta_g <- scale * (g_q_pwindow - g_q) Delta_F_T <- pweibull_q_pwindow - pweibull_q # Calculate Q_Splus using the analytical formula diff --git a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan b/inst/stan/functions/primary_censored_dist_analytical_cdf.stan index 28041e7..6030fa2 100644 --- a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan +++ b/inst/stan/functions/primary_censored_dist_analytical_cdf.stan @@ -136,12 +136,11 @@ real primary_censored_lognormal_uniform_lcdf(data real d, real q, array[] real p * @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) / k - (t/λ) * exp(-(t/λ)^k) + * @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 / scale, shape); - real log_gamma_term = log(gamma_p(1 + inv(shape), x)) - log(shape); - return log_gamma_term; + return log(gamma_p(1 + inv(shape), x)); } /** From 422b23607181f23f6dc2489ef3a53a6600bb6ad8 Mon Sep 17 00:00:00 2001 From: Sam Date: Tue, 8 Oct 2024 11:56:21 +0100 Subject: [PATCH 10/24] keep testing we expect to be numeric numeric --- R/primary_censored_dist.R | 9 ++++----- tests/testthat/test-primary_censored_dist.R | 2 +- tests/testthat/test-rpd-primarycensoreddist.R | 6 ++++++ 3 files changed, 11 insertions(+), 6 deletions(-) diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R index a95c492..60ef7eb 100644 --- a/R/primary_censored_dist.R +++ b/R/primary_censored_dist.R @@ -284,10 +284,9 @@ primary_censored_cdf.pcens_pweibull_dunif <- function( g <- function(t) { # Use the lower incomplete gamma function scaled_t <- (t * inv_scale)^shape - gamma_1k <- vapply(scaled_t, function(x) { - pracma::gammainc(1 + inv_shape, x)["lowinc"] + vapply(scaled_t, function(x) { + pracma::gammainc(x, 1 + inv_shape)["lowinc"] }, numeric(1)) - return(gamma_1k) } # Adjust q so that we have [q-pwindow, q] @@ -308,12 +307,12 @@ primary_censored_cdf.pcens_pweibull_dunif <- function( g_q_pwindow <- g(non_zero_q + pwindow) Q_T <- 1 - pweibull_q_pwindow - Delta_g <- scale * (g_q_pwindow - g_q) + 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 + - (1 / pwindow) * Delta_g - + (scale / pwindow) * Delta_g - (non_zero_q / pwindow) * Delta_F_T # Compute the CDF as 1 - Q_Splus diff --git a/tests/testthat/test-primary_censored_dist.R b/tests/testthat/test-primary_censored_dist.R index d97d258..f6717aa 100644 --- a/tests/testthat/test-primary_censored_dist.R +++ b/tests/testthat/test-primary_censored_dist.R @@ -67,7 +67,7 @@ test_that( 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) 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) From 82ddc33a9aecf519355cfaf664f94dbc3e347d22 Mon Sep 17 00:00:00 2001 From: Sam Date: Tue, 8 Oct 2024 12:00:15 +0100 Subject: [PATCH 11/24] add some zero safety --- R/primary_censored_dist.R | 3 +++ inst/stan/functions/primary_censored_dist_analytical_cdf.stan | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R index 60ef7eb..299499b 100644 --- a/R/primary_censored_dist.R +++ b/R/primary_censored_dist.R @@ -285,6 +285,9 @@ primary_censored_cdf.pcens_pweibull_dunif <- function( # Use the lower incomplete gamma function scaled_t <- (t * inv_scale)^shape vapply(scaled_t, function(x) { + if (x <= 0) { + return(0) + } pracma::gammainc(x, 1 + inv_shape)["lowinc"] }, numeric(1)) } diff --git a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan b/inst/stan/functions/primary_censored_dist_analytical_cdf.stan index 6030fa2..ff4d715 100644 --- a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan +++ b/inst/stan/functions/primary_censored_dist_analytical_cdf.stan @@ -126,7 +126,7 @@ real primary_censored_lognormal_uniform_lcdf(data real d, real q, array[] real p } /** - * Compute the log of the integral of the Weibull PDF from 0 to t + * 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 From 0b42089692e66507f73ecade3c8cfc6fa7b30299 Mon Sep 17 00:00:00 2001 From: Sam Date: Tue, 8 Oct 2024 12:03:53 +0100 Subject: [PATCH 12/24] add tangent testing gap --- tests/testthat/test-fitdistdoublecens.R | 72 +++++++++++++++++++------ 1 file changed, 56 insertions(+), 16 deletions(-) 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" + ) + } +) From bfc15f836760600c32053eb72d1bce13635828ff Mon Sep 17 00:00:00 2001 From: Sam Date: Tue, 8 Oct 2024 13:33:49 +0100 Subject: [PATCH 13/24] add t > 0 enforcement for lower incomplete gamma --- R/primary_censored_dist.R | 4 +--- tests/testthat/test-primary_censored_dist.R | 4 ++-- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R index 299499b..98fdb92 100644 --- a/R/primary_censored_dist.R +++ b/R/primary_censored_dist.R @@ -283,11 +283,9 @@ primary_censored_cdf.pcens_pweibull_dunif <- function( g <- function(t) { # Use the lower incomplete gamma function + t <- pmax(t, 0) scaled_t <- (t * inv_scale)^shape vapply(scaled_t, function(x) { - if (x <= 0) { - return(0) - } pracma::gammainc(x, 1 + inv_shape)["lowinc"] }, numeric(1)) } diff --git a/tests/testthat/test-primary_censored_dist.R b/tests/testthat/test-primary_censored_dist.R index f6717aa..7f022f1 100644 --- a/tests/testthat/test-primary_censored_dist.R +++ b/tests/testthat/test-primary_censored_dist.R @@ -258,8 +258,8 @@ test_that( dprimary <- dunif shapes <- c(0.5, 1, 2, 3) - scales <- c(0.5, 1, 2, 5) - pwindows <- c(1, 2, 5, 10) + scales <- c(0.5, 1, 2) + pwindows <- c(1, 2, 5) for (shape in shapes) { for (scale in scales) { From 9bc66ce9ede1a810d341663b0db04c1737c52b56 Mon Sep 17 00:00:00 2001 From: Sam Date: Tue, 8 Oct 2024 13:35:41 +0100 Subject: [PATCH 14/24] flip args --- R/primary_censored_dist.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R index 98fdb92..eaff243 100644 --- a/R/primary_censored_dist.R +++ b/R/primary_censored_dist.R @@ -286,7 +286,7 @@ primary_censored_cdf.pcens_pweibull_dunif <- function( t <- pmax(t, 0) scaled_t <- (t * inv_scale)^shape vapply(scaled_t, function(x) { - pracma::gammainc(x, 1 + inv_shape)["lowinc"] + pracma::gammainc(1 + inv_shape, x)["lowinc"] }, numeric(1)) } From edbf8822d99e2460d59b1b2ee12d797cf31f29e6 Mon Sep 17 00:00:00 2001 From: Sam Date: Tue, 8 Oct 2024 13:44:33 +0100 Subject: [PATCH 15/24] all R code fix weibull zero bounding --- R/primary_censored_dist.R | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R index eaff243..7fb1db2 100644 --- a/R/primary_censored_dist.R +++ b/R/primary_censored_dist.R @@ -283,7 +283,6 @@ primary_censored_cdf.pcens_pweibull_dunif <- function( g <- function(t) { # Use the lower incomplete gamma function - t <- pmax(t, 0) scaled_t <- (t * inv_scale)^shape vapply(scaled_t, function(x) { pracma::gammainc(1 + inv_shape, x)["lowinc"] @@ -300,12 +299,14 @@ primary_censored_cdf.pcens_pweibull_dunif <- function( # 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 <- pmax(non_zero_q, 0) # Compute necessary survival and distribution functions pweibull_q <- partial_pweibull(non_zero_q) - pweibull_q_pwindow <- partial_pweibull(non_zero_q + pwindow) + pweibull_q_pwindow <- partial_pweibull(q_pwindow) g_q <- g(non_zero_q) - g_q_pwindow <- g(non_zero_q + pwindow) + g_q_pwindow <- g(q_pwindow) Q_T <- 1 - pweibull_q_pwindow Delta_g <- (g_q_pwindow - g_q) From 5ed19a9d363fb641c8cca8506161a736663e8c40 Mon Sep 17 00:00:00 2001 From: Sam Date: Tue, 8 Oct 2024 13:51:25 +0100 Subject: [PATCH 16/24] use unconstrained time outside of integrals and cdfs --- R/primary_censored_dist.R | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R index 7fb1db2..0a1b43d 100644 --- a/R/primary_censored_dist.R +++ b/R/primary_censored_dist.R @@ -285,7 +285,7 @@ primary_censored_cdf.pcens_pweibull_dunif <- function( # Use the lower incomplete gamma function scaled_t <- (t * inv_scale)^shape vapply(scaled_t, function(x) { - pracma::gammainc(1 + inv_shape, x)["lowinc"] + pracma::gammainc(x, 1 + inv_shape)["lowinc"] }, numeric(1)) } @@ -300,16 +300,16 @@ primary_censored_cdf.pcens_pweibull_dunif <- function( if (!all(zero_cases)) { non_zero_q <- q[!zero_cases] q_pwindow <- pmax(non_zero_q + pwindow, 0) - non_zero_q <- pmax(non_zero_q, 0) + non_zero_q_pos <- pmax(non_zero_q, 0) # Compute necessary survival and distribution functions - pweibull_q <- partial_pweibull(non_zero_q) + pweibull_q <- partial_pweibull(non_zero_q_pos) pweibull_q_pwindow <- partial_pweibull(q_pwindow) - g_q <- g(non_zero_q) + 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_g <- g_q_pwindow - g_q Delta_F_T <- pweibull_q_pwindow - pweibull_q # Calculate Q_Splus using the analytical formula From e700efccd732f6383563f1f788538137029e8433 Mon Sep 17 00:00:00 2001 From: Sam Date: Tue, 8 Oct 2024 18:37:15 +0100 Subject: [PATCH 17/24] fall back to the numerical solution when pwindow is large --- R/primary_censored_dist.R | 17 +++++++++++++++-- tests/testthat/test-primary_censored_dist.R | 4 ++-- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R index 0a1b43d..b088895 100644 --- a/R/primary_censored_dist.R +++ b/R/primary_censored_dist.R @@ -263,6 +263,12 @@ primary_censored_cdf.pcens_pweibull_dunif <- function( ) } + 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 @@ -284,9 +290,16 @@ primary_censored_cdf.pcens_pweibull_dunif <- function( g <- function(t) { # Use the lower incomplete gamma function scaled_t <- (t * inv_scale)^shape - vapply(scaled_t, function(x) { - pracma::gammainc(x, 1 + inv_shape)["lowinc"] + 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] diff --git a/tests/testthat/test-primary_censored_dist.R b/tests/testthat/test-primary_censored_dist.R index 7f022f1..8fcdb0e 100644 --- a/tests/testthat/test-primary_censored_dist.R +++ b/tests/testthat/test-primary_censored_dist.R @@ -257,9 +257,9 @@ test_that( dprimary_name <- "dunif" dprimary <- dunif - shapes <- c(0.5, 1, 2, 3) + shapes <- c(0.5, 1, 2) scales <- c(0.5, 1, 2) - pwindows <- c(1, 2, 5) + pwindows <- c(1, 2, 3, 4, 5) for (shape in shapes) { for (scale in scales) { From f1e2606bd787c74af4939644fd2f3da439862b59 Mon Sep 17 00:00:00 2001 From: Sam Date: Tue, 8 Oct 2024 18:57:20 +0100 Subject: [PATCH 18/24] add on stan weibull testing --- .../stan/functions/primary_censored_dist_analytical_cdf.stan | 5 ++++- .../test-stan-primary_censored_dist_analytical_cdf.R | 5 +++++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan b/inst/stan/functions/primary_censored_dist_analytical_cdf.stan index ff4d715..48f88ab 100644 --- a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan +++ b/inst/stan/functions/primary_censored_dist_analytical_cdf.stan @@ -140,7 +140,7 @@ real primary_censored_lognormal_uniform_lcdf(data real d, real q, array[] real p */ real log_weibull_g(real t, real shape, real scale) { real x = pow(t / scale, shape); - return log(gamma_p(1 + inv(shape), x)); + return log(gamma_p(x, 1 + inv(shape))); } /** @@ -229,6 +229,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/tests/testthat/test-stan-primary_censored_dist_analytical_cdf.R b/tests/testthat/test-stan-primary_censored_dist_analytical_cdf.R index 8993d86..e7dcd70 100644 --- a/tests/testthat/test-stan-primary_censored_dist_analytical_cdf.R +++ b/tests/testthat/test-stan-primary_censored_dist_analytical_cdf.R @@ -27,6 +27,8 @@ test_that( 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) } } @@ -63,6 +65,8 @@ test_that( 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) } } @@ -100,6 +104,7 @@ test_that( ) }, numeric(1)) + stan_result <- ifelse(is.nan(stan_result), 1, stan_result) expect_equal(r_result, stan_result, tolerance = 1e-6) } } From c150d2e16a147ff0c62817e6e498433dfb1cf468 Mon Sep 17 00:00:00 2001 From: Sam Date: Tue, 8 Oct 2024 18:58:15 +0100 Subject: [PATCH 19/24] turn it on everywhere --- inst/stan/functions/primary_censored_dist_analytical_cdf.stan | 1 + 1 file changed, 1 insertion(+) diff --git a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan b/inst/stan/functions/primary_censored_dist_analytical_cdf.stan index 48f88ab..a929abc 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 } From 3bf842c49db35a5fbb76ca27a5d734f00a23d3ed Mon Sep 17 00:00:00 2001 From: Sam Date: Wed, 9 Oct 2024 12:36:41 +0100 Subject: [PATCH 20/24] make it easier to find problems in testing --- .../primary_censored_dist_analytical_cdf.stan | 2 +- ...tan-primary_censored_dist_analytical_cdf.R | 27 ++++++++++++++++--- 2 files changed, 25 insertions(+), 4 deletions(-) diff --git a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan b/inst/stan/functions/primary_censored_dist_analytical_cdf.stan index a929abc..2877595 100644 --- a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan +++ b/inst/stan/functions/primary_censored_dist_analytical_cdf.stan @@ -141,7 +141,7 @@ real primary_censored_lognormal_uniform_lcdf(data real d, real q, array[] real p */ real log_weibull_g(real t, real shape, real scale) { real x = pow(t / scale, shape); - return log(gamma_p(x, 1 + inv(shape))); + return log(gamma_p(1 + inv(shape), x)); } /** diff --git a/tests/testthat/test-stan-primary_censored_dist_analytical_cdf.R b/tests/testthat/test-stan-primary_censored_dist_analytical_cdf.R index e7dcd70..7cd4085 100644 --- a/tests/testthat/test-stan-primary_censored_dist_analytical_cdf.R +++ b/tests/testthat/test-stan-primary_censored_dist_analytical_cdf.R @@ -29,7 +29,14 @@ test_that( }, numeric(1)) stan_result <- ifelse(is.nan(stan_result), 1, stan_result) - expect_equal(r_result, stan_result, tolerance = 1e-6) + expect_equal( + r_result, stan_result, + tolerance = 1e-6, + info = sprintf( + "Mismatch for shape = %s, rate = %s, pwindow = %s", + shape, rate, pwindow + ) + ) } } } @@ -67,7 +74,14 @@ test_that( }, numeric(1)) stan_result <- ifelse(is.nan(stan_result), 1, stan_result) - expect_equal(r_result, stan_result, tolerance = 1e-6) + expect_equal( + r_result, stan_result, + tolerance = 1e-6, + info = sprintf( + "Mismatch for meanlog = %s, sdlog = %s, pwindow = %s", + meanlog, sdlog, pwindow + ) + ) } } } @@ -105,7 +119,14 @@ test_that( }, numeric(1)) stan_result <- ifelse(is.nan(stan_result), 1, stan_result) - expect_equal(r_result, stan_result, tolerance = 1e-6) + expect_equal( + r_result, stan_result, + tolerance = 1e-6, + info = sprintf( + "Mismatch for shape = %s, scale = %s, pwindow = %s", + shape, scale, pwindow + ) + ) } } } From 2e35942b2bd71ee17d970467b9b2f70b32e69a0e Mon Sep 17 00:00:00 2001 From: Sam Date: Wed, 9 Oct 2024 12:59:19 +0100 Subject: [PATCH 21/24] unnormalise the lower incomplete gamma --- .../stan/functions/primary_censored_dist_analytical_cdf.stan | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan b/inst/stan/functions/primary_censored_dist_analytical_cdf.stan index 2877595..ae96298 100644 --- a/inst/stan/functions/primary_censored_dist_analytical_cdf.stan +++ b/inst/stan/functions/primary_censored_dist_analytical_cdf.stan @@ -140,8 +140,9 @@ real primary_censored_lognormal_uniform_lcdf(data real d, real q, array[] real p * @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 / scale, shape); - return log(gamma_p(1 + inv(shape), x)); + real x = pow(t * inv(scale), shape); + real a = 1 + inv(shape); + return log(gamma_p(a, x)) + lgamma(a); } /** From 42d5dcf695ef69937105204cdd385682cc69375d Mon Sep 17 00:00:00 2001 From: Sam Date: Wed, 9 Oct 2024 13:12:02 +0100 Subject: [PATCH 22/24] test invalid weibull passing --- tests/testthat/test-primary_censored_dist.R | 25 +++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/tests/testthat/test-primary_censored_dist.R b/tests/testthat/test-primary_censored_dist.R index 8fcdb0e..5a0bdf7 100644 --- a/tests/testthat/test-primary_censored_dist.R +++ b/tests/testthat/test-primary_censored_dist.R @@ -136,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" + ) } ) From 88b5fc0cbb08da2efaf06cd531808ab1c5a2825a Mon Sep 17 00:00:00 2001 From: Sam Date: Wed, 9 Oct 2024 13:14:10 +0100 Subject: [PATCH 23/24] add our required skipping for stan tests --- .../test-stan-primary_censored_dist_analytical_cdf.R | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/tests/testthat/test-stan-primary_censored_dist_analytical_cdf.R b/tests/testthat/test-stan-primary_censored_dist_analytical_cdf.R index 7cd4085..7698746 100644 --- a/tests/testthat/test-stan-primary_censored_dist_analytical_cdf.R +++ b/tests/testthat/test-stan-primary_censored_dist_analytical_cdf.R @@ -1,3 +1,9 @@ +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", From 60a29607299e09eb427984c1402287eca9754d2b Mon Sep 17 00:00:00 2001 From: Sam Date: Wed, 9 Oct 2024 13:17:44 +0100 Subject: [PATCH 24/24] clean up the news --- NEWS.md | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/NEWS.md b/NEWS.md index c8aa81d..08e6fe0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,15 +4,11 @@ This is the development version of `primarycensoreddist` and is not yet ready fo ## Package -* Added R and stan implementations of the primary censored cdf for the weibull distribution with uniform primary censoring. - -## Documentation - -* 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 @@ -20,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