Clean up v0.5.0: extract helpers, remove dead eps field, delete stale artifacts

- Add require_params(), get_delta() helpers and parameter length
  validation to R/utils.R; replace inline boilerplate across
  dfr_dist.R, diagnostics.R, and distributions.R
- Remove dead eps field from dfr_dist constructor and tests
- Simplify inv_cdf (remove unnecessary do.call wrapper)
- Update sampler roxygen to reflect inverse CDF implementation
- Add left-censoring note to distributions.Rd
- Add test-utils.R with 14 tests for the new helpers
- Delete stale vignette .R tangling artifacts and .html builds

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Author: queelius

R/dfr_dist.R

@@ -7,7 +7,6 @@
 #' @param rate A function that computes the hazard rate at time `t`.
 #' @param par The parameters of the distribution. Defaults to `NULL`,
 #'            which means that the parameters are unknown.
-#' @param eps The epsilon update for numerical integration. Defaults to 0.01.
 #' @param ob_col The column name for observation times in data frames.
 #'               Defaults to "t".
 #' @param delta_col The column name for event indicators in data frames.
@@ -25,14 +24,13 @@
 #'                If NULL, falls back to numerical Hessian via numDeriv::hessian.
 #' @return A `dfr_dist` object that inherits from `likelihood_model`.
 #' @export
-dfr_dist <- function(rate, par = NULL, eps = 0.01,
+dfr_dist <- function(rate, par = NULL,
                      ob_col = "t", delta_col = "delta",
                      cum_haz_rate = NULL, score_fn = NULL,
                      hess_fn = NULL) {
     structure(
         list(rate = rate,
              par = par,
-             eps = eps,
              ob_col = ob_col,
              delta_col = delta_col,
              cum_haz_rate = cum_haz_rate,
@@ -87,13 +85,11 @@ inv_cdf.dfr_dist <- function(x, ...) {
     cdf_fn <- cdf(x, ...)
     function(p, par = NULL, ...) {
         par <- get_params(par, x$par)
-        uniroot_args <- list(
-            f = function(t) {
-                cdf_fn(t, par, ...) - p
-            },
+        uniroot(
+            f = function(t) cdf_fn(t, par, ...) - p,
             interval = c(0, 1e3),
-            extendInt = "upX")
-        do.call(uniroot, uniroot_args)$root
+            extendInt = "upX"
+        )$root
     }
 }
 
@@ -109,23 +105,16 @@ params.dfr_dist <- function(x, ...) {
 }
 
 #' Sampling function for `dfr_dist` objects.
-#' 
-#' Since S(t,par) = exp(-cum_hz(t,par)), we can sample from the
-#' distribution by letting t = 0 (or some other positive number if
-#' we want to condition on T > t_min), sampling from an exponential
-#' distribution with `lambda = rate(t, par)`, and then rejecting
-#' the sample if `runif(1) > S(t, par)`. If accepted, add that
-#' observation to the sample, otherwise reject it, let `t = t + eps`
-#' where `eps` is some small number, and repeat. We continue this
-#' process until we have `n` observations for the sample.
+#'
+#' Uses inverse CDF sampling: generates uniform random values and
+#' transforms them through the quantile function (inverse CDF).
 #'
 #' @param x The object to obtain the sampler of.
-#' @param ... Additional arguments to pass into the survival function
+#' @param ... Additional arguments to pass into the inverse CDF constructor.
 #' @return A function that samples from the distribution. It accepts
-#' `n`, the number of samples to take, `t` is the time at which to start
-#' sampling, `par` are the parameters of the distribution, and `eps` is
-#' the  update for numerical integration. Finally, we pass additional
-#' arguments `...` into the hazard function.
+#' `n`, the number of samples to take, `par`, the parameters of the
+#' distribution, and `...`, additional arguments passed to the quantile
+#' function.
 #' @importFrom algebraic.dist params surv sampler
 #' @importFrom stats runif
 #' @export
@@ -313,12 +302,7 @@ loglik.dfr_dist <- function(model, ...) {
         par <- get_params(par, model$par)
 
         t <- df[[model$ob_col]]
-
-        if (model$delta_col %in% names(df)) {
-            delta <- df[[model$delta_col]]
-        } else {
-            delta <- rep(1, length(t))
-        }
+        delta <- get_delta(df, model$delta_col)
 
         ll <- 0
 
@@ -476,13 +460,7 @@ fit.dfr_dist <- function(object, ...) {
     function(df, par = NULL,
              method = c("BFGS", "Nelder-Mead", "L-BFGS-B", "CG", "SANN"),
              control = list(), ...) {
-        if (is.null(par)) {
-            par <- params(object)
-            if (is.null(par)) {
-                stop("Initial parameters required: specify 'par' argument ",
-                     "or set parameters in dfr_dist()")
-            }
-        }
+        par <- require_params(par, params(object))
 
         control <- modifyList(list(fnscale = -1), control)
         method <- match.arg(method)

R/diagnostics.R

@@ -75,12 +75,7 @@ residuals.dfr_dist <- function(object, data, par = NULL,
     type <- match.arg(type)
     H <- cum_haz(object, ...)
 
-    if (is.null(par)) {
-        par <- object$par
-        if (is.null(par)) {
-            stop("Parameters required: provide via 'par' argument or in distribution object")
-        }
-    }
+    par <- require_params(par, object$par)
 
     if (!object$ob_col %in% names(data)) {
         stop(sprintf("Time column '%s' not found in data", object$ob_col))
@@ -93,13 +88,7 @@ residuals.dfr_dist <- function(object, data, par = NULL,
         return(H_vals)
     }
 
-    if (object$delta_col %in% names(data)) {
-        delta <- data[[object$delta_col]]
-    } else {
-        delta <- rep(1, nrow(data))
-    }
-
-    delta - H_vals
+    get_delta(data, object$delta_col) - H_vals
 }
 
 # =============================================================================
@@ -172,12 +161,7 @@ plot.dfr_dist <- function(x, data = NULL, par = NULL,
                            empirical_col = "steelblue", ...) {
     what <- match.arg(what)
 
-    if (is.null(par)) {
-        par <- x$par
-        if (is.null(par)) {
-            stop("Parameters required: provide via 'par' argument or in distribution object")
-        }
-    }
+    par <- require_params(par, x$par)
 
     if (is.null(xlim)) {
         xlim <- if (!is.null(data) && x$ob_col %in% names(data)) {
@@ -210,11 +194,7 @@ plot.dfr_dist <- function(x, data = NULL, par = NULL,
 
     if (!is.null(data) && empirical && !add && x$ob_col %in% names(data)) {
         t <- data[[x$ob_col]]
-        delta <- if (x$delta_col %in% names(data)) {
-            data[[x$delta_col]]
-        } else {
-            rep(1, length(t))
-        }
+        delta <- get_delta(data, x$delta_col)
 
         km <- kaplan_meier(t, delta)
 

R/distributions.R

@@ -4,6 +4,14 @@
 #' Each provides the complete specification (rate, cum_haz_rate, score_fn,
 #' and where practical, hess_fn) for optimal performance.
 #'
+#' @section Left-Censoring Note:
+#' The analytical score and Hessian functions provided by these constructors
+#' assume event indicators in \{0, 1\} (right-censored and exact observations).
+#' For left-censored data (delta = -1), these functions are not applicable and
+#' the package automatically falls back to numerical differentiation via
+#' \code{numDeriv::grad} and \code{numDeriv::hessian} through the log-likelihood,
+#' which handles all censoring types correctly.
+#'
 #' @name distributions
 #' @family distributions
 NULL
@@ -67,11 +75,11 @@ dfr_exponential <- function(lambda = NULL) {
             par[[1]] * t
         },
         score_fn = function(df, par, ...) {
-            delta <- if ("delta" %in% names(df)) df$delta else rep(1, nrow(df))
+            delta <- get_delta(df)
             c(sum(delta == 1) / par[[1]] - sum(df$t))
         },
         hess_fn = function(df, par, ...) {
-            delta <- if ("delta" %in% names(df)) df$delta else rep(1, nrow(df))
+            delta <- get_delta(df)
             matrix(-sum(delta == 1) / par[[1]]^2, nrow = 1, ncol = 1)
         },
         par = lambda
@@ -158,7 +166,7 @@ dfr_weibull <- function(shape = NULL, scale = NULL) {
             k <- par[[1]]
             sigma <- par[[2]]
             t <- df$t
-            delta <- if ("delta" %in% names(df)) df$delta else rep(1, nrow(df))
+            delta <- get_delta(df)
 
             n_events <- sum(delta == 1)
             t_ratio <- t / sigma
@@ -175,7 +183,7 @@ dfr_weibull <- function(shape = NULL, scale = NULL) {
             k <- par[[1]]
             sigma <- par[[2]]
             t <- df$t
-            delta <- if ("delta" %in% names(df)) df$delta else rep(1, nrow(df))
+            delta <- get_delta(df)
 
             n_events <- sum(delta == 1)
             t_ratio <- t / sigma
@@ -260,7 +268,7 @@ dfr_gompertz <- function(a = NULL, b = NULL) {
             a <- par[[1]]
             b <- par[[2]]
             t <- df$t
-            delta <- if ("delta" %in% names(df)) df$delta else rep(1, nrow(df))
+            delta <- get_delta(df)
 
             exp_bt <- exp(b * t)
             n_events <- sum(delta == 1)
@@ -345,7 +353,7 @@ dfr_loglogistic <- function(alpha = NULL, beta = NULL) {
             alpha <- par[[1]]
             beta <- par[[2]]
             t <- df$t
-            delta <- if ("delta" %in% names(df)) df$delta else rep(1, nrow(df))
+            delta <- get_delta(df)
 
             t_ratio <- t / alpha
             t_ratio_beta <- t_ratio^beta

R/utils.R

@@ -2,8 +2,26 @@ get_params <- function(par, default = NULL) {
     if (is.null(par) || all(is.na(par))) {
         return(default)
     }
+    if (!is.null(default) && length(par) != length(default)) {
+        stop(sprintf(
+            "Parameter length mismatch: got %d, expected %d",
+            length(par), length(default)
+        ))
+    }
     if (!is.null(default)) {
         par[is.na(par)] <- default[is.na(par)]
     }
     par
 }
+
+require_params <- function(par, default) {
+    par <- get_params(par, default)
+    if (is.null(par)) {
+        stop("Parameters required: provide via 'par' argument or in distribution object")
+    }
+    par
+}
+
+get_delta <- function(df, delta_col = "delta") {
+    if (delta_col %in% names(df)) df[[delta_col]] else rep(1, nrow(df))
+}

man/dfr_dist.Rd

@@ -37,25 +37,6 @@ test_that("dfr_dist constructor stores parameters", {
   expect_equal(dist$par, c(lambda = 0.5))
 })
 
-test_that("dfr_dist constructor uses default eps", {
-  dist <- dfr_dist(
-    rate = function(t, par, ...) par[1],
-    par = c(1)
-  )
-
-  expect_equal(dist$eps, 0.01)
-})
-
-test_that("dfr_dist constructor accepts custom eps", {
-  dist <- dfr_dist(
-    rate = function(t, par, ...) par[1],
-    par = c(1),
-    eps = 0.001
-  )
-
-  expect_equal(dist$eps, 0.001)
-})
-
 test_that("is_dfr_dist returns TRUE for dfr_dist objects", {
   dist <- make_exponential_dfr(lambda = 1)
   expect_true(is_dfr_dist(dist))

man/distributions.Rd

@@ -584,7 +584,7 @@ test_that("fit throws error when no parameters available", {
 
   solver <- fit(dist)
   # Should error when calling without par
-  expect_error(solver(df), "Initial parameters required")
+  expect_error(solver(df), "Parameters required")
 })
 
 # =============================================================================