Issue 85: Add infrastructure for analytical solutions in R

NAMESPACE

@@ -1,11 +1,16 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(primary_censored_cdf,default)
+S3method(primary_censored_cdf,pcens_numeric)
+S3method(primary_censored_cdf,pcens_pgamma_dunif)
+S3method(primary_censored_cdf,pcens_plnorm_dunif)
 export(check_dprimary)
 export(check_pdist)
 export(dexpgrowth)
 export(dpcens)
 export(dprimarycensoreddist)
 export(fitdistdoublecens)
+export(new_primary_censored_dist)
 export(pcd_as_cmdstan_data)
 export(pcd_cmdstan_model)
 export(pcd_load_stan_functions)
@@ -15,6 +20,7 @@ export(pcd_stan_path)
 export(pexpgrowth)
 export(ppcens)
 export(pprimarycensoreddist)
+export(primary_censored_cdf)
 export(rexpgrowth)
 export(rpcens)
 export(rprimarycensoreddist)

NEWS.md

@@ -6,6 +6,7 @@ This is the current development version.
 
 * Add `{touchstone}` based benchmarks for benchmarking R utility functions, and fitting the `stan` and `fitdistplus` models.
 * Added a "How it works" vignette.
+* Added R infrastructure for analytical solutions via the `primary_censored_dist` S3 class.
 
 # primarycensoreddist 0.4.0
 

R/dprimarycensoreddist.R

@@ -70,7 +70,8 @@
 dprimarycensoreddist <- function(
     x, pdist, pwindow = 1, swindow = 1,
     D = Inf, dprimary = stats::dunif,
-    dprimary_args = list(), log = FALSE, ...) {
+    dprimary_args = list(), log = FALSE,
+    pdist_name = NULL, dprimary_name = NULL, ...) {
   check_pdist(pdist, D, ...)
   check_dprimary(dprimary, pwindow, dprimary_args)
 
@@ -82,6 +83,13 @@ dprimarycensoreddist <- function(
     )
   }
 
+  if (is.null(pdist_name)) {
+    pdist_name <- .extract_function_name(substitute(pdist))
+  }
+  if (is.null(dprimary_name)) {
+    dprimary_name <- .extract_function_name(substitute(dprimary))
+  }
+
   # Compute CDFs for all unique points
   unique_points <- sort(unique(c(x, x + swindow)))
   unique_points <- unique_points[unique_points > 0]
@@ -90,7 +98,8 @@ dprimarycensoreddist <- function(
   }
 
   cdfs <- pprimarycensoreddist(
-    unique_points, pdist, pwindow, Inf, dprimary, dprimary_args, ...
+    unique_points, pdist, pwindow, Inf, dprimary, dprimary_args,
+    pdist_name = pdist_name, dprimary_name = dprimary_name, ...
   )
 
   # Create a lookup table for CDFs

R/pprimarycensoreddist.R

@@ -26,7 +26,19 @@
 #' @param dprimary_args List of additional arguments to be passed to
 #' dprimary. For example, when using `dexpgrowth`, you would
 #' pass `list(min = 0, max = pwindow, r = 0.2)` to set the minimum, maximum,
-#' and rate parameters.
+#' and rate parameters
+#'
+#' @param pdist_name A string specifying the name of the delay distribution
+#' function. If NULL, the function name is extracted using
+#' [.extract_function_name()]. Used to determine if a analytical solution
+#' exists for the primary censored distribution. Must be set if `pdist` is
+#' passed a pre-assigned variable rather than a function name.
+#'
+#' @param dprimary_name A string specifying the name of the primary event
+#' distribution function. If NULL, the function name is extracted using
+#' [.extract_function_name()]. Used to determine if a analytical solution
+#' exists for the primary censored distribution. Must be set if `dprimary` is
+#' passed a pre-assigned variable rather than a function name.
 #'
 #' @param ... Additional arguments to be passed to pdist
 #'
@@ -61,7 +73,19 @@
 #' }
 #' where \eqn{F_{\text{cens,norm}}(q)} is the normalized CDF.
 #'
+#' This function creates a `primary_censored_dist` object using
+#' [new_primary_censored_dist()] and then computes the primary event
+#' censored CDF using [primary_censored_cdf()]. This abstraction allows
+#' for automatic use of analytical solutions when available, while
+#' seamlessly falling back to numerical integration when necessary.
+#'
+#' Note: For analytical detection to work correctly, `pdist` and `dprimary`
+#' must be directly passed as distribution functions, not via assignment or
+#' `pdist_name` and `dprimary_name` must be used to override the default
+#' extraction of the function name.
+#'
 #' @family primarycensoreddist
+#' @seealso [new_primary_censored_dist()] and [primary_censored_cdf()]
 #'
 #' @examples
 #' # Example: Lognormal distribution with uniform primary events
@@ -75,25 +99,29 @@
 #' )
 pprimarycensoreddist <- function(
     q, pdist, pwindow = 1, D = Inf, dprimary = stats::dunif,
-    dprimary_args = list(), ...) {
+    dprimary_args = list(), pdist_name = NULL, dprimary_name = NULL, ...) {
   check_pdist(pdist, D, ...)
   check_dprimary(dprimary, pwindow, dprimary_args)
 
-  result <- vapply(q, function(d) {
-    if (d <= 0) {
-      return(0) # Return 0 for non-positive delays
-    } else {
-      integrand <- function(p) {
-        d_adj <- d - p
-        pdist(d_adj, ...) *
-          do.call(
-            dprimary, c(list(x = p, min = 0, max = pwindow), dprimary_args)
-          )
-      }
+  if (is.null(pdist_name)) {
+    pdist_name <- .extract_function_name(substitute(pdist))
+  }
+  if (is.null(dprimary_name)) {
+    dprimary_name <- .extract_function_name(substitute(dprimary))
+  }
 
-      stats::integrate(integrand, lower = 0, upper = pwindow)$value
-    }
-  }, numeric(1))
+  # Create a new primary_censored_dist object
+  pcens_obj <- new_primary_censored_dist(
+    pdist,
+    dprimary,
+    dprimary_args,
+    pdist_name = pdist_name,
+    dprimary_name = dprimary_name,
+    ...
+  )
+
+  # Compute the CDF using the S3 method
+  result <- primary_censored_cdf(pcens_obj, q, pwindow)
 
   if (!is.infinite(D)) {
     # Compute normalization factor for finite D

R/primary_censored_dist.R

@@ -0,0 +1,162 @@
+#' S3 class for primary event censored distribution computation
+#'
+#' @inheritParams pprimarycensoreddist
+#'
+#' @return An object of class primary_censored_cdf
+#'
+#' @family primary_censored_dist
+#'
+#' @export
+new_primary_censored_dist <- function(
+    pdist, dprimary, dprimary_args,
+    pdist_name = NULL,
+    dprimary_name = NULL, ...) {
+  if (is.null(pdist_name)) {
+    pdist_name <- .extract_function_name(substitute(pdist))
+  }
+  if (is.null(dprimary_name)) {
+    dprimary_name <- .extract_function_name(substitute(dprimary))
+  }
+
+  structure(
+    list(
+      pdist = pdist,
+      dprimary = dprimary,
+      dprimary_args = dprimary_args,
+      args = list(...)
+    ),
+    class = c(
+      paste0(
+        "pcens_",
+        pdist_name, "_",
+        dprimary_name
+      )
+    )
+  )
+}
+
+#' Compute primary event censored CDF
+#'
+#' @inheritParams pprimarycensoreddist
+#'
+#' @param object A `primary_censored_dist` object as created by
+#' [new_primary_censored_dist()].
+#'
+#' @param 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.
+#'
+#' @return Vector of primary event censored CDFs
+#'
+#' @family primary_censored_dist
+#'
+#' @export
+primary_censored_cdf <- function(
+    object, q, pwindow, use_numeric = FALSE) {
+  UseMethod("primary_censored_cdf")
+}
+
+#' Default method for computing primary event censored CDF
+#'
+#' This method serves as a fallback for combinations of delay and primary
+#' event distributions that don't have specific implementations. It uses
+#' the numeric integration method.
+#'
+#' @inheritParams primary_censored_cdf
+#'
+#' @family primary_censored_dist
+#'
+#' @export
+primary_censored_cdf.default <- function(
+    object, q, pwindow, use_numeric = FALSE) {
+  primary_censored_cdf.pcens_numeric(object, q, pwindow, use_numeric)
+}
+
+#' Numeric method for computing primary event censored CDF
+#'
+#' This method uses numerical integration to compute the primary event censored
+#' CDF for any combination of delay distribution and primary event distribution.
+#'
+#' @inheritParams primary_censored_cdf
+#' @inheritParams pprimarycensoreddist
+#'
+#' @details
+#' This method implements the numerical integration approach for computing
+#' the primary event censored CDF. It uses the same mathematical formulation
+#' as described in the details section of [pprimarycensoreddist()], but
+#' applies numerical integration instead of analytical solutions.
+#'
+#' @seealso [pprimarycensoreddist()] for the mathematical details of the
+#' primary event censored CDF computation.
+#'
+#' @family primary_censored_dist
+#'
+#' @export
+primary_censored_cdf.pcens_numeric <- function(
+    object, q, pwindow, use_numeric = FALSE) {
+  result <- vapply(q, function(d) {
+    if (d <= 0) {
+      return(0) # Return 0 for non-positive delays
+    } else {
+      integrand <- function(p) {
+        d_adj <- d - p
+        do.call(object$pdist, c(list(q = d_adj), object$args)) *
+          do.call(
+            object$dprimary,
+            c(list(x = p, min = 0, max = pwindow), object$dprimary_args)
+          )
+      }
+
+      stats::integrate(integrand, lower = 0, upper = pwindow)$value
+    }
+  }, numeric(1))
+
+  return(result)
+}
+
+#' Method for Gamma delay with uniform primary
+#'
+#' @inheritParams primary_censored_cdf
+#'
+#' @family primary_censored_dist
+#'
+#' @export
+primary_censored_cdf.pcens_pgamma_dunif <- function(
+    object, q, pwindow, use_numeric = FALSE) {
+  use_numeric <- TRUE
+  if (isTRUE(use_numeric)) {
+    return(
+      primary_censored_cdf.pcens_numeric(object, q, pwindow, use_numeric)
+    )
+  }
+
+  result <- vapply(q, function(n) {
+    # Implement analytical solution here
+  }, numeric(1))
+
+  return(result)
+}
+
+#' Method for Log-Normal delay with uniform primary
+#'
+#' @inheritParams primary_censored_cdf
+#'
+#' @family primary_censored_dist
+#'
+#' @export
+primary_censored_cdf.pcens_plnorm_dunif <- function(
+    object, q, pwindow, use_numeric = FALSE) {
+  use_numeric <- TRUE
+  if (isTRUE(use_numeric)) {
+    return(
+      primary_censored_cdf.pcens_numeric(object, q, pwindow, use_numeric)
+    )
+  }
+
+  result <- vapply(q, function(n) {
+    # Implement analytical solution here
+  }, numeric(1))
+
+  return(result)
+}

R/utils.R

@@ -0,0 +1,15 @@
+#' Extract Base Function Name
+#'
+#' This helper function extracts the base name of a function, removing an
+#' namespace prefixes.
+#'
+#' @param func The output of `substitute` on a function.
+#'
+#' @return A character string representing the base name of the function.
+#'
+#' @keywords internal
+.extract_function_name <- function(func) {
+  func_name <- deparse(func)
+  base_name <- sub("^.*::", "", func_name)
+  return(base_name)
+}

_pkgdown.yml

@@ -43,6 +43,10 @@ reference:
   desc: Probability density and random generation functions for primary event distributions
   contents:
   - has_concept("primaryeventdistributions")
+- title: Primary censored distribution class and methods
+  desc: S3 class and methods for computing primary event censored distributions, focusing on the internal machinery used by the package. Unlike the primary event distributions section which deals with specific distribution functions, this section covers the general framework for handling censored distributions.
+  contents:
+  - has_concept("primary_censored_dist")
 - title: Distribution checking functions
   desc: Functions to validate cumulative distribution functions (CDFs) and probability density functions (PDFs)
   contents:

inst/WORDLIST

@@ -15,6 +15,7 @@ PMFs
 Poisson
 SamuelBrand
 athowes
+cdf
 cdot
 cens
 cmdstan
@@ -23,6 +24,7 @@ disccensprob
 disccensunifprim
 discretisation
 discretising
+dist
 dp
 dprimary
 dprimarycensoreddist