diff --git a/.Rbuildignore b/.Rbuildignore
index 60657a7..808480e 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -5,7 +5,6 @@
 ^\_config\.yml$
 ^\.github$
 ^CONTRIBUTING\.md$
-^NEWS\.md$
 ^_pkgdown\.yml$
 pkgdown/
 .DS_Store
diff --git a/DESCRIPTION b/DESCRIPTION
index e837edc..6fb62ef 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: bmstate
 Type: Package
 Title: Bayesian multistate modeling
-Version: 0.3.3
+Version: 0.4.0
 Authors@R:
     c(person(given = "Juho",
            family = "Timonen",
diff --git a/NEWS.md b/NEWS.md
new file mode 100644
index 0000000..634cee7
--- /dev/null
+++ b/NEWS.md
@@ -0,0 +1,85 @@
+# bmstate 0.4.0
+
+* **Correct the baseline-hazard prior scale.** The public model API stores
+  baseline-hazard reference values as natural-scale rates, but these values
+  were passed directly to a Stan parameter used as the location of
+  `log_w0`. For example, a reference rate of 0.001 therefore centered the
+  hazard near `exp(0.001) = 1.001`, rather than 0.001. Stan data creation now
+  applies `log()` exactly once at the R--Stan boundary. The default Stan data
+  field is consequently renamed from `mu_w0` to `mu_log_w0`; custom Stan
+  programs selected with `options(bmstate_stan_file = ...)` must adopt the new
+  name.
+  Fits produced by affected versions must be rerun.
+
+* **Make transition support and baseline-hazard calibration explicit.** The
+  previous cumulative-hazard summary silently omitted modeled transitions with
+  no events, which could leave empirical prior locations undefined and defer
+  failure until Stan data were validated. Fit-time validation now stops with
+  the affected transition names if any transition in the fitted matrix has no
+  observed occurrence or no positive at-risk time. This check is deliberately
+  not applied to out-of-sample prediction data. For valid transition systems,
+  reference rates are computed from event counts divided by exact at-risk
+  time, pooled within transition type. Documentation now calls the calibrated
+  quantity a natural-scale reference rate, not an arithmetic prior mean.
+
+* **Preserve joint posterior-draw alignment in prediction.** Replicating the
+  baseline-hazard arrays previously produced a different draw/subject ordering
+  from the hazard multipliers and prediction metadata, so parameters from one
+  posterior draw could be combined with another. The package now constructs an
+  explicit draw-major row index and uses the same ordering for spline weights,
+  log intercepts, subject multipliers, path links, and occupancy calculations.
+
+* **Retain zero state-visit probabilities.** `p_state_visit()` previously built
+  its output only from observed event rows, silently dropping requested
+  subjects or groups with no visit to an event state. It now completes the
+  state-by-group grid, joins event counts and denominators separately, and
+  returns explicit zero counts and probabilities.
+
+* **Keep simulated paths inside the requested horizon.** Applying
+  `min_t_step` could move a genuine event beyond `t_max`; the final-row cleanup
+  then relabeled that row as censoring without clearing its transition index.
+  Event-time adjustment now preserves a genuine pre-horizon event when the
+  adjustment would cross the horizon, every censor row has transition index
+  zero, invalid reversed horizons are rejected, and a zero-length horizon
+  returns a valid one-row censored path.
+
+* **Stabilize the oral one-compartment PK solution near equal rates.** Direct
+  evaluation divided by `ka - ke` and could produce cancellation, `0 / 0`, or
+  `0 * Inf` when absorption and elimination rates were equal or very different.
+  The transient solution now uses sign-aware `expm1()` factorizations, and the
+  transient and steady-state solutions use their continuous equal-rate limits
+  in a small numerical guard region. The limit uses the `ka` prefactor so that
+  both the value and first derivative match at equality to first order.
+
+* **Construct exactly the requested integration grid.** For a noninteger
+  maximum time, the old `seq(..., ceiling(t_max), by = delta)` expression could
+  generate more than `N_grid` points and make the declared Stan dimensions
+  disagree with the supplied arrays. The grid is now formed directly from
+  exactly `N_grid` midpoint indices.
+
+* **Honor delayed starts in transition-probability prediction.**
+  `solve_trans_prob_matrix()` previously initialized the forward equation at
+  the first output time, even when an earlier `t_start` was requested. It now
+  includes `t_start` in the integration times, returns only the requested
+  output rows, and validates times and hazard-parameter dimensions even on the
+  zero-duration identity path.
+
+* **Return the requested PK credible-interval width.** The plotting helper
+  divided the tail probability by two twice, so an 80% request used the 5th and
+  95th percentiles and returned a 90% interval. It now uses
+  `(1 - ci_alpha) / 2` and its complement directly.
+
+* **Reduce reverse-mode work for time-invariant hazard predictors.** The Stan
+  program previously expanded subject covariates to every interval and, before
+  hoisting, rebuilt the complete multiplier matrix inside each transition
+  loop. It now differentiates one predictor per subject and transition type,
+  then gathers those values for the relevant intervals. This is algebraically
+  identical and reduced gradient time by about 38% in a representative
+  benchmark with many intervals, with no measured regression when every
+  interval belonged to a different subject.
+
+* **Tighten edge-case validation and output construction.** The release adds
+  checks for positive PK rates, dosing intervals, concentration bounds, and
+  exposure-normalization scales; rejects malformed delayed-start requests;
+  handles base data frames without accidental dimension dropping; and builds
+  simulated path output in a list before binding it once.
diff --git a/R/DosingData.R b/R/DosingData.R
index 5a76b73..789fb57 100644
--- a/R/DosingData.R
+++ b/R/DosingData.R
@@ -67,6 +67,9 @@ PSSDosingData <- R6::R6Class(
       self$set_sub_ids(subject_ids)
       N_sub <- self$num_subjects()
       checkmate::assert_number(tau_ss, lower = 0)
+      if (tau_ss <= 0) {
+        stop("tau_ss must be strictly positive")
+      }
       checkmate::assert_list(doses, len = N_sub)
       checkmate::assert_list(times, len = N_sub)
       self$doses <- doses
@@ -269,6 +272,9 @@ simulate_dosing <- function(df_subjects, tau = 24, p_miss = 0.2, t_jitter = 4) {
 pk_2cpt_pss <- function(t, dose_ss, times, doses, theta, tau, MAX_CONC) {
   ensure_exposed_stan_functions()
   checkmate::assert_number(tau, lower = 0)
+  if (tau <= 0) {
+    stop("tau must be strictly positive")
+  }
   checkmate::assert_number(MAX_CONC, lower = 0)
   checkmate::assert_numeric(dose_ss, lower = 0)
   N_sub <- length(dose_ss)
diff --git a/R/MultistateModel.R b/R/MultistateModel.R
index 020b4e1..dc3f158 100644
--- a/R/MultistateModel.R
+++ b/R/MultistateModel.R
@@ -110,7 +110,7 @@ MultistateModel <- R6::R6Class("MultistateModel",
 
       log_w0 <- matrix(rep(log_w0, N), N, S, byrow = TRUE)
       log_m <- private$simulate_log_hazard_multipliers(df_subjects, beta_haz)
-      paths <- self$system$simulate(w_all, log_w0, log_m, min_t_step = 0.1)
+      paths <- self$system$simulate(w_all, log_w0, log_m, min_t_step = 1e-6)
       as_tibble(paths)
     }
   ),
@@ -163,8 +163,11 @@ MultistateModel <- R6::R6Class("MultistateModel",
     #' @param loc Location
     #' @param scale Scale
     set_xpsr_normalizers = function(loc = 0, scale = 1) {
-      checkmate::assert_numeric(loc, lower = 0, len = 1)
-      checkmate::assert_numeric(scale, lower = 0, len = 1)
+      checkmate::assert_number(loc, finite = TRUE)
+      checkmate::assert_number(scale, finite = TRUE)
+      if (scale <= 0) {
+        stop("xpsr normalization scale must be strictly positive")
+      }
       message(
         "setting xpsr normalizers to loc = ",
         round(loc, 5), ", scale = ", round(scale, 5)
@@ -174,35 +177,43 @@ MultistateModel <- R6::R6Class("MultistateModel",
       invisible(NULL)
     },
 
-    #' @description Get assumed prior mean baseline hazard rates.
-    #' @return Numeric vector with length equal to number of transitions
+    #' @description Get baseline-hazard prior reference rates.
+    #' @return A positive numeric vector with length equal to the number of
+    #' transitions. These natural-scale rates are transformed to locations on
+    #' the log-hazard scale when Stan data are created; they are not arithmetic
+    #' prior means.
     get_prior_mean_h0 = function() {
       v <- private$prior_mean_h0
       if (is.null(v)) {
-        stop("prior mean h0 has not been set")
+        stop("prior reference h0 has not been set")
       }
       v
     },
 
-    #' @description Set assumed prior mean baseline hazard rates (side
-    #' effect).
-    #' @param mean_h0 Numeric vector with length equal to number of transitions
+    #' @description Set baseline-hazard prior reference rates (side effect).
+    #' @param mean_h0 A positive numeric vector of natural-scale reference rates
+    #' with length equal to the number of transitions. The argument name is
+    #' retained for compatibility; the values are not arithmetic prior means.
     set_prior_mean_h0 = function(mean_h0) {
       N_trans <- self$system$tm()$num_trans()
-      checkmate::assert_numeric(mean_h0, len = N_trans, lower = 0)
+      checkmate::assert_numeric(mean_h0, len = N_trans, finite = TRUE)
+      if (any(mean_h0 <= 0)) {
+        stop("prior baseline hazard rates must be strictly positive")
+      }
       private$prior_mean_h0 <- mean_h0
       invisible(NULL)
     },
 
-    #' @description Set assumed prior mean baseline hazard rates (side
-    #' effect) based on average hazards in data.
+    #' @description Set baseline-hazard prior reference rates (side effect)
+    #' using event counts and at-risk time pooled within transition type.
     #' @param data A \code{\link{JointData}} or \code{\link{PathData}} object.
     set_prior_mean_h0_data = function(data) {
       if (inherits(data, "PathData")) {
         data <- JointData$new(data, NULL)
       }
       checkmate::assert_class(data, "JointData")
-      df_ttype <- average_haz_per_ttype(data$paths) |> dplyr::arrange(.data$trans_idx)
+      df_ttype <- average_haz_per_ttype(data$paths) |>
+        dplyr::arrange(.data$trans_idx)
       self$set_prior_mean_h0(exp(df_ttype$log_h0_avg))
     },
 
diff --git a/R/MultistateModelFit.R b/R/MultistateModelFit.R
index c16d01c..bb3387a 100644
--- a/R/MultistateModelFit.R
+++ b/R/MultistateModelFit.R
@@ -78,6 +78,9 @@ MultistateModelFit <- R6::R6Class("MultistateModelFit",
     #'
     #' @param name Param/quantity name
     get_draws = function(name = NULL) {
+      if (is.null(name)) {
+        return(private$draws)
+      }
       private$draws[[name]]
     },
 
@@ -261,7 +264,7 @@ MultistateModelFit <- R6::R6Class("MultistateModelFit",
     #' @description
     #' Full names of parameters that start with \code{log_z_}.
     log_z_pars = function() {
-      nams <- names(self$draws())
+      nams <- self$draws_names()
       match <- grepl(nams, pattern = "log_z_")
       nams[which(match)]
     },
@@ -508,14 +511,6 @@ msmfit_log_hazard_multipliers <- function(fit, oos = FALSE, data = NULL) {
     beta_xpsr <- array(0, dim = c(S, 1, 0, sd$N_trans_types))
   }
 
-  # Create x_haz_long (long version of hazard covariates vector)
-  N_sub <- sd$N_sub
-  first_indices <- sapply(seq_len(N_sub), function(x) which(sd$idx_sub == x)[1])
-  if (sd$nc_haz > 0) {
-    x_haz_long <- sd$x_haz[, sd$idx_sub, drop = FALSE]
-  } else {
-    x_haz_long <- array(0, dim = c(0, sd$N_int))
-  }
   an <- fit$model$get_xpsr_normalizers()
 
   # Call exposed Stan function for each draw (not optimal)
@@ -523,7 +518,7 @@ msmfit_log_hazard_multipliers <- function(fit, oos = FALSE, data = NULL) {
   for (s in seq_len(S)) {
     if (sd$do_pk == 1) {
       ba <- list(beta_xpsr[s, 1, 1, ])
-      aa <- list(xpsr[[s]][sd$idx_sub])
+      aa <- list(xpsr[[s]])
     } else {
       ba <- NULL
       aa <- NULL
@@ -531,17 +526,17 @@ msmfit_log_hazard_multipliers <- function(fit, oos = FALSE, data = NULL) {
     if (sd$nc_haz == 0 && sd$do_pk == 0) {
       r <- matrix(0, sd$N_sub, sd$N_trans)
     } else {
-      r <- compute_log_hazard_multiplier(
-        sd$N_int,
+      r_by_type <- compute_log_hazard_multiplier(
+        sd$N_sub,
+        sd$N_trans_types,
         mat2list(t(array(beta_oth[s, 1, , ], dim = c(sd$nc_haz, sd$N_trans_types)))),
         ba,
-        mat2list(t(x_haz_long)),
+        mat2list(t(sd$x_haz)),
         aa,
         an$loc,
-        an$scale,
-        sd$ttype
+        an$scale
       )
-      r <- r[first_indices, , drop = FALSE]
+      r <- r_by_type[, sd$ttype, drop = FALSE]
     }
     out[[s]] <- r
   }
@@ -581,6 +576,21 @@ msmfit_log_baseline_hazard <- function(fit, t = NULL) {
 }
 
 
+# Repeat each joint posterior draw for every subject while preserving draw-major
+# order. Keeping this in one helper makes the required parameter alignment
+# explicit and independently testable.
+repeat_hazard_parameter_draws <- function(w, log_w0, n_subjects) {
+  checkmate::assert_array(w, d = 3)
+  checkmate::assert_matrix(log_w0, nrows = dim(w)[1])
+  checkmate::assert_integerish(n_subjects, len = 1, lower = 1)
+  draw_rows <- rep(seq_len(dim(w)[1]), each = n_subjects)
+  list(
+    w = w[draw_rows, , , drop = FALSE],
+    log_w0 = log_w0[draw_rows, , drop = FALSE]
+  )
+}
+
+
 #' Extract and reshape draws of instant hazard related parameters
 #'
 #' @export
@@ -600,8 +610,9 @@ msmfit_inst_hazard_param_draws <- function(fit, oos = FALSE, data = NULL) {
   log_w0 <- fit$get_draws_of("log_w0")
   w <- access_one_dim(w, dim_index = 2, value = 1)
   log_w0 <- access_one_dim(log_w0, dim_index = 2, value = 1)
-  w_rep <- abind::abind(replicate(N, w, simplify = FALSE), along = 1)
-  log_w0_rep <- abind::abind(replicate(N, log_w0, simplify = FALSE), along = 1)
+  repeated <- repeat_hazard_parameter_draws(w, log_w0, N)
+  w_rep <- repeated$w
+  log_w0_rep <- repeated$log_w0
   log_m_reshaped <- do.call(rbind, log_m)
 
   # Subject-draw df
@@ -646,8 +657,6 @@ generate_paths <- function(fit, oos = FALSE, t_start = 0, t_max = NULL, n_rep =
   checkmate::assert_class(fit, "MultistateModelFit")
   checkmate::assert_integerish(n_rep, lower = 1, len = 1)
   sd <- msmfit_stan_data(fit, data)
-  log_m <- msmfit_log_hazard_multipliers(fit, oos, data)
-
   # Get and reshape draws and state vector at t_start
   sys <- fit$model$system
   S <- fit$num_draws()
@@ -729,9 +738,13 @@ p_state_occupancy <- function(fit, oos = FALSE, t_start = 0, t_out = NULL,
 
   sys <- fit$model$system
   if (is.null(t_out)) {
-    t_out <- seq(t_start, sys$get_tmax(), length.out = 30)
+    if (t_start == sys$get_tmax()) {
+      t_out <- t_start
+    } else {
+      t_out <- seq(t_start, sys$get_tmax(), length.out = 30)
+    }
   } else {
-    checkmate::assert_numeric(t_out, min.len = 2)
+    checkmate::assert_numeric(t_out, min.len = 1, finite = TRUE)
   }
 
   # Get and reshape draws
diff --git a/R/MultistateSystem.R b/R/MultistateSystem.R
index 3078bf0..9b80bc0 100644
--- a/R/MultistateSystem.R
+++ b/R/MultistateSystem.R
@@ -16,6 +16,15 @@ MultistateSystem <- R6::R6Class("MultistateSystem",
         w <- matrix(w, nrow = 1)
       }
       checkmate::assert_array(w, d = 2)
+      if (t_start == t_max) {
+        return(cbind(
+          time = t_start,
+          state = init_state,
+          is_event = 0,
+          is_censor = 1,
+          trans_idx = 0
+        ))
+      }
 
       # Setup
       states <- init_state
@@ -33,8 +42,27 @@ MultistateSystem <- R6::R6Class("MultistateSystem",
           states[j], t_cur, t_max, w, log_w0, log_m
         )
         t_next <- trans$t
-        if (t_next - t_cur < min_t_step) {
-          t_next <- t_cur + min_t_step
+        if (trans$idx > 0) {
+          original_t_next <- t_next
+          if (t_next <= t_cur) {
+            # Ensure representable forward progress even when an exponential
+            # waiting time rounds to zero and min_t_step is zero.
+            numerical_step <- 2 * .Machine$double.eps * max(1, abs(t_cur))
+            t_next <- t_cur + max(min_t_step, numerical_step)
+          } else if (t_next - t_cur < min_t_step) {
+            t_next <- t_cur + min_t_step
+          }
+          if (t_next >= t_max) {
+            if (original_t_next > t_cur && original_t_next < t_max) {
+              # Do not move a genuine pre-horizon event to or beyond censoring.
+              t_next <- original_t_next
+            } else {
+              # No representable interior event time remains.
+              t_next <- t_max
+              trans$idx <- 0
+              trans$new_state <- 0
+            }
+          }
         }
 
         # Update time and state
@@ -57,6 +85,7 @@ MultistateSystem <- R6::R6Class("MultistateSystem",
       is_event[1] <- 0 # initial state is never an event
       is_event[L] <- 0 # last state is never an event (it is censoring time)
       states[L] <- states[L - 1]
+      tidx[L] <- 0
       cbind(
         time = times,
         state = states,
@@ -69,18 +98,28 @@ MultistateSystem <- R6::R6Class("MultistateSystem",
     # Generate transition given state and transition intensity functions
     generate_transition = function(state, t_init, t_max, w, log_w0, log_m) {
       tol <- 1.05
+      rare_hazard_cutoff <- 1e-9
       possible <- private$transmat$possible_transitions_from(state)
-      J <- private$transmat$num_trans()
-      UB <- rep(0, J)
-      for (j in seq_len(J)) {
-        UB[j] <- tol * self$max_inst_hazard(t_init, t_max, w[j, ], log_w0[j], log_m[j])
-      }
-      possible <- possible[which(UB[possible] > 1e-9)] # so rare that not possible
       if (length(possible) == 0) {
         # Absorbing state, no transitions possible
         return(list(t = t_max, new_state = 0, idx = 0))
       }
-      UB <- UB[possible]
+      UB <- tol * vapply(possible, function(j) {
+        self$max_inst_hazard(t_init, t_max, w[j, ], log_w0[j], log_m[j])
+      }, numeric(1))
+      if (any(!is.finite(UB))) {
+        stop("non-finite transition-hazard upper bound")
+      }
+      # This is an intentional simulation approximation: transitions below
+      # this absolute rate threshold are treated as impossible to avoid doing
+      # proposal work for events that are negligible on the package's intended
+      # time scale.
+      keep <- UB > rare_hazard_cutoff
+      possible <- possible[keep]
+      UB <- UB[keep]
+      if (length(possible) == 0) {
+        return(list(t = t_max, new_state = 0, idx = 0))
+      }
       w <- w[possible, , drop = FALSE]
       log_w0 <- log_w0[possible]
       log_m <- log_m[possible]
@@ -107,12 +146,6 @@ MultistateSystem <- R6::R6Class("MultistateSystem",
         new_state <- 0
       }
 
-      # Safeguard against infinite loop
-      dt_min <- 1e-9
-      if (t_min_found - t_init < dt_min) {
-        t_min_found <- t_init + dt_min
-      }
-
       # Return
       list(
         t = t_min_found, new_state = new_state, idx = trans_idx
@@ -130,7 +163,7 @@ MultistateSystem <- R6::R6Class("MultistateSystem",
         n_tries <- n_tries + 1
         u1 <- stats::runif(1)
         t <- t - 1 / lambda_ub * log(u1)
-        if (t > t_max) {
+        if (t >= t_max) {
           return(list(t = t_max, n_tries = n_tries, censor = TRUE))
         }
         lambda_t <- exp(self$log_inst_hazard(t, w, log_w0, log_m))
@@ -276,7 +309,11 @@ MultistateSystem <- R6::R6Class("MultistateSystem",
       knots <- self$get_knots()
       L <- length(knots)
       BK <- knots[c(1, L)]
-      knots <- knots[2:(L - 1)]
+      if (L > 2) {
+        knots <- knots[2:(L - 1)]
+      } else {
+        knots <- numeric(0)
+      }
       bspline_basis(t, private$spline_k, knots, BK)
     },
 
@@ -284,7 +321,8 @@ MultistateSystem <- R6::R6Class("MultistateSystem",
     #'
     #' @param t Time point(s). Not used if \code{SBF} is given.
     #' @param w Spline basis function weights (vector)
-    #' @param log_w0 Intercept (log)
+    #' @param log_w0 Intercept (log), or \code{-Inf} for an identically zero
+    #' baseline hazard
     #' @param log_m Hazard multiplier (log)
     #' @param SBF Pre-computed basis function matrix at \code{t}.
     log_inst_hazard = function(t, w, log_w0, log_m, SBF = NULL) {
@@ -331,20 +369,36 @@ MultistateSystem <- R6::R6Class("MultistateSystem",
       sum(diag(self$tm()$matrix)) > 0
     },
 
-    #' @description Max instant hazard on interval (t1, t2)
+    #' @description Approximate maximum instantaneous hazard on an interval,
+    #' evaluated on a 100-point grid. Path generation inflates this value by
+    #' five percent before thinning.
     #'
     #' @param t1 Start time point
     #' @param t2 End time point
     #' @param w Spline basis function weights (vector)
-    #' @param log_w0 Intercept (log)
+    #' @param log_w0 Intercept (log), or \code{-Inf} for an identically zero
+    #' baseline hazard
     #' @param log_m Hazard multiplier (log)
     max_inst_hazard = function(t1, t2, w, log_w0, log_m) {
-      ttt <- seq(t1, t2, length.out = 100)
-      log_haz <- log_m + self$log_baseline_hazard(ttt, log_w0, w)
-      return(exp(max(log_haz)))
+      checkmate::assert_number(t1, lower = min(self$get_knots()))
+      checkmate::assert_number(t2, lower = t1, upper = self$get_tmax())
+      checkmate::assert_numeric(w, len = self$num_weights(), finite = TRUE)
+      checkmate::assert_number(log_w0, upper = .Machine$double.xmax)
+      checkmate::assert_number(log_m, finite = TRUE)
+      if (is.infinite(log_w0)) {
+        return(0)
+      }
+
+      # Fast numerical envelope used by path generation. The caller inflates
+      # this grid maximum by five percent before thinning.
+      t_grid <- seq(t1, t2, length.out = 100)
+      log_haz <- log_m + self$log_baseline_hazard(t_grid, log_w0, w)
+      exp(max(log_haz))
     },
 
-    #' @description Generate paths
+    #' @description Generate paths. For computational efficiency, a transition
+    #' whose five-percent-inflated numerical hazard envelope is at most
+    #' \code{1e-9} is treated as impossible.
     #'
     #' @param w An array of shape \code{n_draws} x \code{n_trans} x
     #' \code{n_weights}
@@ -356,7 +410,8 @@ MultistateSystem <- R6::R6Class("MultistateSystem",
     #' @param t_max Max time. If \code{NULL}, the max
     #' time of the model is used.
     #' @param n_rep Number of repetitions to do for each draw.
-    #' @param min_t_step Minimal time step.
+    #' @param min_t_step Minimum separation applied to consecutive simulated
+    #' event times when this does not move an event to or beyond \code{t_max}.
     #' @return A data frame with \code{n_draws} x \code{n_rep} paths.
     simulate = function(w, log_w0, log_m, init_state = 1, t_start = 0,
                         t_max = NULL, n_rep = 1, min_t_step = 1e-6) {
@@ -379,11 +434,13 @@ MultistateSystem <- R6::R6Class("MultistateSystem",
       pb <- progress::progress_bar$new(total = n_paths)
 
       # Set max time
-      out <- NULL
+      out <- vector("list", n_paths)
       if (is.null(t_max)) {
         t_max <- self$get_tmax()
       }
-      checkmate::assert_number(t_max, lower = 0)
+      checkmate::assert_number(
+        t_max, lower = t_start, upper = self$get_tmax(), finite = TRUE
+      )
 
       # Should not be done in parallel as such because can mess order in link df
       cnt <- 0
@@ -397,10 +454,10 @@ MultistateSystem <- R6::R6Class("MultistateSystem",
             min_t_step
           )
           p <- cbind(p, rep(cnt, nrow(p)))
-          out <- rbind(out, p)
+          out[[cnt]] <- p
         }
       }
-      df <- data.frame(out)
+      df <- data.frame(do.call(rbind, out))
       colnames(df)[ncol(df)] <- "path_id"
       df
     }
@@ -469,9 +526,15 @@ solve_trans_prob_matrix <- function(system, t_out, log_w0, w = NULL,
                                     log_m = NULL, t_start = 0,
                                     ...) {
   checkmate::assert_class(system, "MultistateSystem")
-  checkmate::assert_numeric(t_out)
+  checkmate::assert_numeric(t_out, min.len = 1, finite = TRUE)
   K <- length(t_out)
   checkmate::assert_number(t_start, lower = 0)
+  if (any(t_out < t_start)) {
+    stop("all t_out values must be greater than or equal to t_start")
+  }
+  if (is.unsorted(t_out, strictly = TRUE)) {
+    stop("t_out must be strictly increasing")
+  }
   H <- system$num_trans()
   W <- system$num_weights()
   S <- system$num_states()
@@ -481,10 +544,19 @@ solve_trans_prob_matrix <- function(system, t_out, log_w0, w = NULL,
   if (is.null(log_m)) {
     log_m <- rep(0, H)
   }
-  kfe <- solve_time_evolution(system, t_out, log_w0, w, log_m, ...)
+  checkmate::assert_matrix(w, ncols = W, nrows = H)
+  checkmate::assert_numeric(log_w0, len = H)
+  checkmate::assert_numeric(log_m, len = H)
+  solve_times <- unique(c(t_start, t_out))
+  if (length(solve_times) == 1) {
+    P <- array(diag(1, S, S), dim = c(1, S, S))
+    return(P)
+  }
+  kfe <- solve_time_evolution(system, solve_times, log_w0, w, log_m, ...)
+  output_rows <- match(t_out, kfe[, 1])
   P <- array(0, dim = c(K, S, S))
   for (k in seq_len(K)) {
-    P[k, , ] <- matrix(kfe[k, 2:ncol(kfe)], S, S)
+    P[k, , ] <- matrix(kfe[output_rows[k], 2:ncol(kfe)], S, S)
   }
   P
 }
diff --git a/R/PKModel.R b/R/PKModel.R
index dac6723..290070e 100644
--- a/R/PKModel.R
+++ b/R/PKModel.R
@@ -23,6 +23,9 @@ PKModel <- R6::R6Class("PKModel",
     #' @param value Upper bound for concentration, to avoid numerical issues.
     set_max_conc = function(value) {
       checkmate::assert_number(value, lower = 0)
+      if (value <= 0) {
+        stop("concentration upper bound must be strictly positive")
+      }
       message("setting max conc = ", round(value, 5))
       private$MAX_CONC <- value
       invisible(NULL)
@@ -90,16 +93,23 @@ PKModel <- R6::R6Class("PKModel",
       checkmate::assert_number(theta$V2, lower = 0)
       checkmate::assert_number(dose, lower = 0)
       checkmate::assert_number(tau, lower = 0)
+      if (theta$ka <= 0 || theta$CL <= 0 || theta$V2 <= 0 || tau <= 0) {
+        stop("ka, CL, V2, and tau must be strictly positive")
+      }
       N <- length(t)
       ka <- theta$ka
       ke <- theta$CL / theta$V2
-      A <- (dose / theta$V2) * (ka / (ka - ke))
-      conc <- rep(0, N)
-      ma <- A / (-expm1(-ka * tau))
-      me <- A / (-expm1(-ke * tau))
       tt <- t %% tau
-      conc <- me * exp(-ke * tt) - ma * exp(-ka * tt)
-      conc
+      d <- ka - ke
+      rel_tol <- 1e-8 * max(abs(ka), abs(ke))
+      if (abs(d) <= rel_tol) {
+        k <- 0.5 * (ka + ke)
+        f <- exp(-k * tt) / (-expm1(-k * tau))
+        return((dose / theta$V2) * ka * f * (tt + tau / expm1(k * tau)))
+      }
+      f_ke <- exp(-ke * tt) / (-expm1(-ke * tau))
+      f_ka <- exp(-ka * tt) / (-expm1(-ka * tau))
+      (dose / theta$V2) * ka / d * (f_ke - f_ka)
     },
 
     #' @description Compute exposure, which is the steady-state
diff --git a/R/PathData.R b/R/PathData.R
index 36d52e5..68bf9db 100644
--- a/R/PathData.R
+++ b/R/PathData.R
@@ -85,7 +85,7 @@ PathData <- R6::R6Class(
       checkmate::assert_integerish(path_df$path_id)
 
       # As tibbles
-      subject_df <- as_tibble(subject_df[, cols1])
+      subject_df <- as_tibble(subject_df[, cols1, drop = FALSE])
       path_df <- as_tibble(path_df[, cols2])
       link_df <- as_tibble(link_df[, cols3])
 
@@ -444,8 +444,8 @@ msfit_average_hazard <- function(msfit) {
   msfit$Haz |>
     dplyr::group_by(.data$trans) |>
     summarise(
-      avg_haz = (dplyr::last(.data$Haz) - dplyr::first(.data$Haz)) /
-        (dplyr::last(.data$time) - dplyr::first(.data$time))
+      avg_haz = dplyr::last(.data$Haz) / dplyr::last(.data$time),
+      .groups = "drop"
     )
 }
 
@@ -650,29 +650,50 @@ p_state_visit <- function(pd, t = NULL, by = NULL) {
   }
   checkmate::assert_numeric(t, lower = 0)
 
-  if (!is.null(by)) {
-    checkmate::assert_character(by, len = 1)
-    c <- pd$as_data_frame(covariates = by) |>
-      dplyr::group_by(.data$state, .data[[by]])
-  } else {
-    c <- pd$as_data_frame() |>
-      dplyr::group_by(.data$state)
-  }
   estates <- which(pd$state_names() %in% pd$get_event_state_names())
-  df <- count_paths_with_event(c, t, S) |> dplyr::filter(.data$state %in% estates)
+  events <- pd$as_data_frame(covariates = by) |>
+    dplyr::filter(
+      .data$time <= t,
+      .data$state %in% estates
+    )
+
   if (!is.null(by)) {
-    df_all <- c |> dplyr::ungroup()
-    df_all <- df_all |> dplyr::group_by(.data[[by]])
-    df_all <- df_all |>
-      dplyr::distinct(.data$path_id) |>
-      dplyr::count()
-    df <- df |> dplyr::left_join(df_all, by = by)
+    checkmate::assert_character(by, len = 1)
+    links <- pd$full_link(covariates = by) |>
+      dplyr::distinct(.data$path_id, .data[[by]])
+    groups <- links |>
+      dplyr::distinct(.data[[by]])
+    grid <- tidyr::expand_grid(
+      state_idx = estates,
+      group_value = groups[[by]]
+    )
+    names(grid)[2] <- by
+    counts <- events |>
+      dplyr::transmute(
+        path_id = .data$path_id,
+        state_idx = .data$state,
+        group_value = .data[[by]]
+      ) |>
+      dplyr::distinct() |>
+      dplyr::count(.data$state_idx, .data$group_value, name = "n_event")
+    names(counts)[2] <- by
+    totals <- links |>
+      dplyr::count(.data[[by]], name = "n")
+    df <- grid |>
+      dplyr::left_join(counts, by = c("state_idx", by)) |>
+      dplyr::left_join(totals, by = by)
+    df$n_event[is.na(df$n_event)] <- 0L
     df$prob <- df$n_event / df$n
   } else {
+    counts <- events |>
+      dplyr::transmute(path_id = .data$path_id, state_idx = .data$state) |>
+      dplyr::distinct() |>
+      dplyr::count(.data$state_idx, name = "n_event")
+    df <- data.frame(state_idx = estates) |>
+      dplyr::left_join(counts, by = "state_idx")
+    df$n_event[is.na(df$n_event)] <- 0L
     df$prob <- df$n_event / pd$n_paths()
   }
-  df$state_idx <- df$state
-  df$state <- NULL
   sdf <- pd$transmat$states_df()
   df |> dplyr::left_join(sdf, by = "state_idx")
 }
@@ -783,6 +804,7 @@ df_to_paths_df_part2 <- function(pdf, tm) {
 #' @return A \code{\link{PathData}} object
 df_to_pathdata <- function(df, tm, covs = NULL, validate = TRUE) {
   check_columns(df, c("state", "time", "subject_id", "is_transition"))
+  df <- tibble::as_tibble(df)
   checkmate::assert_integerish(df$state)
   checkmate::assert_numeric(df$time)
   checkmate::assert_character(df$subject_id)
diff --git a/R/stan-data.R b/R/stan-data.R
index b67bd64..d84cca4 100644
--- a/R/stan-data.R
+++ b/R/stan-data.R
@@ -31,7 +31,9 @@ create_stan_data_model <- function(model) {
     N_trans = tm$num_trans(),
     N_trans_types = tm$num_trans_types(),
     ttype = tm$trans_df()$trans_type,
-    mu_w0 = model$get_prior_mean_h0(),
+    # The public API specifies prior locations as hazard rates. Stan models
+    # log hazards, so transform exactly once at this boundary.
+    mu_log_w0 = log(model$get_prior_mean_h0()),
     nc_haz = length(model$data_covs("haz")),
     N_sbf = model$system$num_weights()
   )
@@ -85,7 +87,7 @@ create_stan_data_timegrid <- function(model) {
   if (delta_grid > 0.25 * t_max) {
     stop("delta_grid is very large compared to t_max")
   }
-  t_grid <- seq(delta_grid / 2, ceiling(t_max), by = delta_grid) # midpoints
+  t_grid <- (seq_len(G) - 0.5) * delta_grid # exactly G midpoints
 
   # Return
   list(
@@ -101,6 +103,13 @@ create_stan_data_timegrid <- function(model) {
 create_stan_data_transitions <- function(pd) {
   tm <- pd$transmat
   dat <- pd$as_transitions()
+  interval_length <- dat$time - dat$time_prev
+  if (any(interval_length < 0)) {
+    stop("transition data contain a negative at-risk interval")
+  }
+  if (any(dat$trans_idx > 0 & interval_length <= 0)) {
+    stop("an observed transition has no positive at-risk time")
+  }
   N_int <- nrow(dat)
   N_trans <- tm$num_trans()
   transition <- matrix(0, N_trans, N_int)
@@ -227,11 +236,17 @@ standata_scaled_covariates <- function(pd, model, name) {
   for (cn in covs) {
     j <- j + 1
     xx <- sub_df[[cn]]
+    if (!is.numeric(xx)) {
+      stop("covariate '", cn, "' must be numeric before normalization")
+    }
     xj_loc <- norms$locations[[cn]]
     xj_scale <- norms$scales[[cn]]
     if (is.null(xj_loc) || is.null(xj_scale)) {
       stop("no normalizers set for ", cn)
     }
+    if (!is.finite(xj_loc) || !is.finite(xj_scale) || xj_scale <= 0) {
+      stop("normalization scale for '", cn, "' must be finite and positive")
+    }
     x_norm <- (xx - xj_loc) / xj_scale
     check_normalized_covariate(x_norm, cn)
     x[[j]] <- x_norm
@@ -295,22 +310,63 @@ create_stan_data_intervalidx <- function(t_start, t_end, t_grid, delta_grid) {
   )
 }
 
-# Average hazard per transition type, ignoring transitions that did not
-# occur
+# Estimate a constant hazard per transition type by pooling event counts and
+# at-risk time over all transitions of that type. Every modeled transition must
+# have at least one observed event and positive at-risk time.
 average_haz_per_ttype <- function(pd) {
-  msfit <- pd$fit_mstate()
-  h0 <- msfit_average_hazard(msfit) |> dplyr::arrange(.data$trans)
-  df_trans <- pd$transmat$trans_df()
-  df_ttype <- df_trans |> dplyr::select("trans_idx", "trans_type")
-  df_ttype$trans <- df_ttype$trans_idx
-  h0 <- h0 |> left_join(df_ttype, by = "trans")
-  df_mean_log_h0 <- h0 |>
-    dplyr::filter(.data$avg_haz > 0) |>
+  checkmate::assert_class(pd, "PathData")
+  dat <- pd$as_transitions(truncate = TRUE)
+  trans <- pd$transmat$trans_df()
+  exposure <- numeric(nrow(trans))
+  n_event <- integer(nrow(trans))
+
+  for (h in seq_len(nrow(trans))) {
+    at_risk <- dat$from == trans$prev_state[h]
+    exposure[h] <- sum(dat$time[at_risk] - dat$time_prev[at_risk])
+    n_event[h] <- sum(dat$trans_idx == trans$trans_idx[h])
+    if (exposure[h] <= 0 && n_event[h] > 0) {
+      stop("transition ", h, " has events but no positive at-risk time")
+    }
+  }
+
+  unobserved <- which(n_event == 0)
+  if (length(unobserved) > 0) {
+    stop(
+      "cannot calibrate or fit transition(s) with no observed events: ",
+      paste(trans$trans_char[unobserved], collapse = ", "),
+      ". Remove those transitions from the model or provide informative data."
+    )
+  }
+  unexposed <- which(exposure <= 0)
+  if (length(unexposed) > 0) {
+    stop(
+      "cannot calibrate or fit transition(s) with no positive at-risk time: ",
+      paste(trans$trans_char[unexposed], collapse = ", "),
+      ". Remove those transitions from the model or provide informative data."
+    )
+  }
+
+  rate_df <- data.frame(
+    trans_idx = trans$trans_idx,
+    trans_type = trans$trans_type,
+    exposure = exposure,
+    n_event = n_event
+  )
+  type_df <- rate_df |>
     dplyr::group_by(.data$trans_type) |>
-    mutate(log_haz = log(.data$avg_haz)) |>
-    summarize(log_h0_avg = mean(.data$log_haz)) |>
-    dplyr::arrange(.data$trans_type)
-  df_ttype |> left_join(df_mean_log_h0, by = "trans_type")
+    dplyr::summarise(
+      exposure = sum(.data$exposure),
+      n_event = sum(.data$n_event),
+      .groups = "drop"
+    ) |>
+    dplyr::mutate(
+      rate = .data$n_event / .data$exposure,
+      log_h0_avg = log(.data$rate)
+    ) |>
+    dplyr::select("trans_type", "log_h0_avg")
+  rate_df |>
+    dplyr::select("trans_idx", "trans_type") |>
+    dplyr::left_join(type_df, by = "trans_type")
 }
 
 # Edit Stan data, sort of computes which() for the binary vectors of each
diff --git a/R/stan.R b/R/stan.R
index 0403066..7720669 100644
--- a/R/stan.R
+++ b/R/stan.R
@@ -50,7 +50,7 @@ ensure_exposed_stan_functions <- function(...) {
 #'
 #' @description
 #' \emph{NOTE:} This function has a side effect of setting covariate
-#' normalizers, prior assumed mean baseline hazard, and concentration upper
+#' normalizers, baseline-hazard prior reference rates, and concentration upper
 #' bound (PK) based on data.
 #'
 #' @export
@@ -59,8 +59,8 @@ ensure_exposed_stan_functions <- function(...) {
 #' @param return_stanfit Return also the raw 'Stan' fit object?
 #' @inheritParams create_stan_model
 #' @param set_normalizers Set covariate normalization automatically?
-#' @param set_prior_h0 Set prior mean average baseline hazard levels based
-#' on data?
+#' @param set_prior_h0 Set baseline-hazard prior reference rates from event
+#' counts and at-risk time in the data?
 #' @param max_conc_factor Factor to multiply observed max concentration by
 #' to get concentration upper bound.
 #' @param method Must be one of \code{"sample"} (default),
@@ -92,7 +92,7 @@ fit_stan <- function(model, data,
   # Get Stan model object
   stan_model <- create_stan_model(filepath = filepath)
 
-  # Set prior mean baseline hazard rates (side effect)
+  # Set baseline-hazard prior reference rates (side effect)
   if (set_prior_h0) {
     model$set_prior_mean_h0_data(data)
   }
@@ -100,11 +100,24 @@ fit_stan <- function(model, data,
   # Set covariate normalizing locations and scales (side effect)
   if (set_normalizers) {
     model$set_normalizers(data)
-    if (!is.null(data$dosing)) {
+    if (model$has_pk()) {
+      if (is.null(data$dosing)) {
+        stop("a model with a PK component requires dosing data")
+      }
       log_mu_CL <- -2 # should match msm.stan
+      if (any(!is.finite(data$dosing$dose_ss)) || any(data$dosing$dose_ss <= 0)) {
+        stop("steady-state doses must be finite and positive for log-exposure modeling")
+      }
       aaa <- log(data$dosing$dose_ss) - log_mu_CL
       loc <- mean(aaa)
       sca <- stats::sd(aaa)
+      if (!is.finite(sca) || sca <= sqrt(.Machine$double.eps)) {
+        sca <- 0.5
+        warning(
+          "steady-state dose has no usable between-subject variation; ",
+          "using the default xpsr normalization scale of 0.5"
+        )
+      }
       model$set_xpsr_normalizers(loc, sca)
     }
   }
diff --git a/R/utils.R b/R/utils.R
index 3c97474..7942caf 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -1,11 +1,11 @@
 # Helper
 pksim_to_quantiles <- function(sim, ci_alpha) {
   checkmate::assert_number(ci_alpha, lower = 0, upper = 1)
-  av <- (1 - ci_alpha) / 2
+  tail_prob <- (1 - ci_alpha) / 2
   sim <- sim |>
     dplyr::group_by(.data$subject_id, .data$time) |>
     dplyr::summarise(q = list(
-      stats::quantile(.data$val, probs = c(av / 2, 0.5, 1 - av / 2))
+      stats::quantile(.data$val, probs = c(tail_prob, 0.5, 1 - tail_prob))
     ), .groups = "drop") |>
     tidyr::unnest_wider(q, names_sep = "_")
   colnames(sim)[3:5] <- c("lower", "val", "upper")
@@ -117,10 +117,14 @@ sim_subject_ids <- function(N) {
 
 # Truncate path data frame
 truncate_after_terminal_events <- function(df, term_state_inds) {
-  term_events <- df |>
-    dplyr::filter(.data$state %in% term_state_inds, .data$trans_idx > 0) |>
+  term_rows <- df |>
+    dplyr::filter(.data$state %in% term_state_inds, .data$trans_idx > 0)
+  if (nrow(term_rows) == 0) {
+    return(df)
+  }
+  term_events <- term_rows |>
     dplyr::group_by(.data$path_id) |>
-    summarise(term_time = min(.data$time, na.rm = T)) |>
+    summarise(term_time = min(.data$time, na.rm = TRUE)) |>
     dplyr::ungroup()
   no_terms <- df |>
     dplyr::anti_join(term_events, by = "path_id")
@@ -136,6 +140,12 @@ truncate_after_terminal_events <- function(df, term_state_inds) {
 # make sense for given data
 prefit_checks <- function(model, data) {
   pd <- data$paths
+  if (!model$pk_only) {
+    # This validates that every transition retained in the fitted hazard model
+    # has at least one event and positive at-risk time. Prediction data are not
+    # subject to this requirement.
+    average_haz_per_ttype(pd)
+  }
   lens <- pd$as_transitions() |>
     dplyr::mutate(time_len = .data$time - .data$time_prev) |>
     dplyr::pull(.data$time_len)
diff --git a/inst/stan/msm.stan b/inst/stan/msm.stan
index 389dd88..2cedd2b 100644
--- a/inst/stan/msm.stan
+++ b/inst/stan/msm.stan
@@ -7,29 +7,27 @@ functions {
 
   // log of hazard multiplier
   matrix compute_log_hazard_multiplier(
-      data int N_int,
+      data int N_row,
+      data int N_trans_types,
       array[] vector beta_oth,
       array[] vector beta_xpsr,
       data array[] vector x_haz,
       array[] vector x_xpsr,
       real xpsr_loc,
-      real xpsr_scale,
-      data array[] int ttype
+      real xpsr_scale
   ) {
-    int N_trans = size(ttype);
     int nc_haz = size(x_haz);
-    matrix[N_int, N_trans] log_C_haz = rep_matrix(0.0, N_int, N_trans);
-    for(j in 1:N_trans){
-      int h = ttype[j];
-      if(nc_haz > 0){
-        for(k in 1:nc_haz){
-          log_C_haz[,j] += beta_oth[k][h] * x_haz[k];
-        }
-      }
-      if(size(beta_xpsr)==1){
-        log_C_haz[,j] += beta_xpsr[1][h] * ((x_xpsr[1] - xpsr_loc) / xpsr_scale);
+    matrix[N_row, N_trans_types] log_C_haz =
+      rep_matrix(0.0, N_row, N_trans_types);
+    if(nc_haz > 0){
+      for(k in 1:nc_haz){
+        log_C_haz += x_haz[k] * beta_oth[k]';
       }
     }
+    if(size(beta_xpsr)==1){
+      log_C_haz += ((x_xpsr[1] - xpsr_loc) / xpsr_scale) *
+        beta_xpsr[1]';
+    }
     return(log_C_haz);
   }
 
@@ -86,7 +84,22 @@ functions {
       reject("t = ", t,", should be non-negative");
     }
     real ke = CL / V2;
-    return(A0 * ka/(ka-ke) * (exp(-ke*t) - exp(-ka*t)) + C0 * exp(-ke*t));
+    real d = ka - ke;
+    real scale = fmax(abs(ka), abs(ke));
+    real absorbed;
+    // Parameter-dependent numerical guard. The branch approximation is
+    // first-order matched, so its value and gradient agree at ka = ke.
+    if(abs(d) <= 1e-8 * scale){
+      real k = 0.5 * (ka + ke);
+      absorbed = A0 * ka * t * exp(-k*t);
+    } else if(d > 0){
+      absorbed = A0 * ka * exp(-ke*t) * (-expm1(-d*t)) / d;
+    } else {
+      // When ka < ke, factoring by exp(-ke*t) can create 0 * Inf.
+      // Factor by the more slowly decaying exp(-ka*t) instead.
+      absorbed = A0 * ka * exp(-ka*t) * expm1(d*t) / d;
+    }
+    return(absorbed + C0 * exp(-ke*t));
   }
 
   // Analytic solution with general initial condition A0
@@ -102,18 +115,24 @@ functions {
   real two_cpt_central_ss(real t, real tau, real dose, real ka, real CL, real V2)
   {
     real ke = CL / V2;
-    real A = dose * (ka / (ka-ke));
-    real tt = fmod(t, tau);
-    real ma = A * inv(-expm1(-ka*tau)); // 1/(1-exp(-ka*tau))
-    real me = A * inv(-expm1(-ke*tau));
-    return(me * exp(-ke*tt) - ma * exp(-ka * tt));
+    real d = ka - ke;
+    real tt = t - floor(t / tau) * tau;
+    real scale = fmax(abs(ka), abs(ke));
+    // As above, this guard is first-order matched at the equal-rate limit.
+    if(abs(d) <= 1e-8 * scale){
+      real k = 0.5 * (ka + ke);
+      real f = exp(-k*tt) * inv(-expm1(-k*tau));
+      return(dose * ka * f * (tt + tau / expm1(k*tau)));
+    }
+    return(dose * ka / d * (
+      exp(-ke*tt) * inv(-expm1(-ke*tau))
+      - exp(-ka*tt) * inv(-expm1(-ka*tau))
+    ));
   }
 
   // Two-cpt PK model (steady state at trough)
   real two_cpt_central_ss0(real tau, real dose, real ka, real CL, real V2) {
-    real ke = CL / V2;
-    real A = dose * (ka / (ka-ke));
-    return(A * (inv(-expm1(-ke*tau)) - inv(-expm1(-ka*tau))));
+    return(two_cpt_central_ss(0.0, tau, dose, ka, CL, V2));
   }
 
   // Two-cpt PK model (steady state at trough)
@@ -262,7 +281,6 @@ functions {
 
     // Find drug amounts in both compartments at dose_times
     int N_t = num_elements(t[1]);
-    vector[N_t] max_conc = rep_vector(MAX_CONC, N_t);
     int N_sub = num_elements(dose_ss);
     int D = num_elements(dose_times[1]);
     array[2, N_sub] vector[D] amounts = pop_2cpt_partly_ss_stage1(
@@ -294,7 +312,7 @@ data {
   int<lower=0, upper=1> omit_lik_haz; // flag
   int<lower=0, upper=1> omit_lik_pk; //flag
   int<lower=0> nc_haz; // number of hazard covariates
-  vector[N_trans] mu_w0; // Assumed mean h0
+  vector[N_trans] mu_log_w0; // Prior location for log baseline hazard
   int<lower=1> N_grid; // number of integration grid points
   real<lower=0> delta_grid; // grid step size
   matrix[N_grid, N_sbf] SBF_grid; // basis functions evaluated at t_grid
@@ -317,7 +335,7 @@ data {
   int<lower=1> N_int; // total number of intervals
   int<lower=1> D_trans; // max num of occurrences for a transition
   int<lower=1> D_risk; // max num of at-risk intervals for a transition
-  array[N_trans] int<lower=0> sum_trans; // total number of occurred
+  array[N_trans] int<lower=1> sum_trans; // total number of occurred
   array[N_trans] int<lower=1> sum_risk; // total number at risk
   array[N_trans, D_trans] int<lower=0,upper=N_int> which_trans;
   array[N_trans, D_risk] int<lower=0,upper=N_int> which_risk;
@@ -355,18 +373,28 @@ data {
 
 transformed data {
   array[N_sub] vector[1] t0_ss;
-  for(n in 1:N_sub){
-    t0_ss[n] = rep_vector(0.0, 1);
+  if(I_xpsr != do_pk){
+    reject("I_xpsr must equal do_pk");
+  }
+  if(I_xpsr == 1 && xpsr_scale <= 0){
+    reject("xpsr_scale must be strictly positive when exposure is modeled");
   }
-  // Set x corresponding to each interval
-  array[nc_haz] vector[N_int] x_haz_long;
-  if(nc_haz > 0){
-    for(j in 1:nc_haz){
-      for(n in 1:N_int){
-        x_haz_long[j][n] = x_haz[j][idx_sub[n]];
+  if(do_pk == 1){
+    if(tau_ss <= 0){
+      reject("tau_ss must be strictly positive when PK is modeled");
+    }
+    if(MAX_CONC <= 0){
+      reject("MAX_CONC must be strictly positive when PK is modeled");
+    }
+    for(n in 1:N_sub){
+      if(dose_ss[n] <= 0){
+        reject("dose_ss must be strictly positive when PK is modeled");
       }
     }
   }
+  for(n in 1:N_sub){
+    t0_ss[n] = rep_vector(0.0, 1);
+  }
 }
 
 parameters {
@@ -408,7 +436,7 @@ transformed parameters {
     for(j in 1:N_trans){
       weights[1,j] = mu_weights[1][:,ttype[j]] +
         sig_weights[1][:,ttype[j]] .* z_weights[1][:,j];
-      log_w0[1,j] = mu_w0[j] + sig_w0[1] * z_w0[1][j];
+      log_w0[1,j] = mu_log_w0[j] + sig_w0[1] * z_w0[1][j];
     }
   }
 
@@ -416,7 +444,6 @@ transformed parameters {
   array[do_pk] matrix[N_sub, 3] log_theta_pk;
   array[do_pk, N_sub] vector[2] conc_mu_pk;
   array[do_pk] vector[N_sub] ss_xpsr; // xpsr for each subject
-  array[do_pk] vector[N_int] x_xpsr_long; // xpsr for each interval
 
   if(do_pk == 1){
 
@@ -434,11 +461,6 @@ transformed parameters {
 
     // Concentration xpsr at steady state
     ss_xpsr[1] = log_ss_area_under_conc(dose_ss, log_theta_pk[1]);
-
-    // Set xpsr corresponding to each interval
-    for(n in 1:N_int){
-      x_xpsr_long[1][n] = ss_xpsr[1][idx_sub[n]];
-    }
   }
 }
 
@@ -468,7 +490,7 @@ model {
 
   // PK PARAM PRIOR
   if(do_pk == 1){
-    for(n in 1:size(log_z_pk[1])){
+    for(n in 1:N_sub){
       log_z_pk[1, n] ~ normal(0, 1);
     }
     log_mu_pk[1] ~ normal(0, 2);
@@ -482,34 +504,24 @@ model {
 
   // Hazard model likelihood
   if (omit_lik_haz == 0 && do_haz == 1) {
+    // Covariates are time invariant, so differentiate one predictor per
+    // subject and transition type, then gather it for each interval.
+    matrix[N_sub, N_trans_types] log_C_haz = compute_log_hazard_multiplier(
+      N_sub, N_trans_types, beta_oth[1], beta_xpsr[1], x_haz, ss_xpsr,
+      xpsr_loc, xpsr_scale
+    );
     for(h in 1:N_trans){
-
-      // log of hazard multiplier on each interval
-      matrix[N_int, N_trans] log_C_haz = compute_log_hazard_multiplier(
-        N_int, beta_oth[1], beta_xpsr[1], x_haz_long, x_xpsr_long, xpsr_loc, xpsr_scale,
-        ttype
-      );
-
-      // Ragged array access
-      array[sum_risk[h]] int idx_atr = which_risk[h, 1:sum_risk[h]];
       array[sum_trans[h]] int idx_occ = which_trans[h, 1:sum_trans[h]];
-
-      // Occurred transitions (log hazard at interval end time)
       target += log_hazard(
-        log_C_haz[idx_occ, h], SBF_end[idx_occ,:], weights[1,h], log_w0[1,h]
+        log_C_haz[idx_sub[idx_occ], ttype[h]], SBF_end[idx_occ,:],
+        weights[1,h], log_w0[1,h]
       );
-
-      // Evaluate baseline hazard at grid, prepad with zero
+      array[sum_risk[h]] int idx_atr = which_risk[h, 1:sum_risk[h]];
       vector[N_grid+1] h0_grid = rep_vector(0.0, N_grid+1);
+      vector[N_grid+1] h0_int;
       h0_grid[2:(N_grid+1)] = exp(log_basehaz(SBF_grid, weights[1,h], log_w0[1,h]));
-
-      // Integrated baseline hazard at grid
-      vector[N_grid+1] h0_int = cumulative_sum(h0_grid) * delta_grid;
-
-      // Transitions that were at risk (- integrated hazard over interval)
-      // Taking into account that indices are off by one because of
-      // The prepadded zero
-      target += - exp(log_C_haz[idx_atr, h]) .*
+      h0_int = cumulative_sum(h0_grid) * delta_grid;
+      target += - exp(log_C_haz[idx_sub[idx_atr], ttype[h]]) .*
         (h0_int[t_end_idx[idx_atr]] - h0_int[t_start_idx_m1[idx_atr]]) .*
         correction_multiplier[idx_atr];
     }
diff --git a/man/MultiStateModel.Rd b/man/MultiStateModel.Rd
index 5eaa27e..2ae9c68 100644
--- a/man/MultiStateModel.Rd
+++ b/man/MultiStateModel.Rd
@@ -126,21 +126,23 @@ Set normalization constants for exposure (side effect)
 \if{html}{\out{<a id="method-MultistateModel-get_prior_mean_h0"></a>}}
 \if{latex}{\out{\hypertarget{method-MultistateModel-get_prior_mean_h0}{}}}
 \subsection{Method \code{get_prior_mean_h0()}}{
-Get assumed prior mean baseline hazard rates.
+Get baseline-hazard prior reference rates.
 \subsection{Usage}{
 \if{html}{\out{<div class="r">}}\preformatted{MultistateModel$get_prior_mean_h0()}\if{html}{\out{</div>}}
 }
 
 \subsection{Returns}{
-Numeric vector with length equal to number of transitions
+A positive numeric vector with length equal to the number of
+transitions. These natural-scale rates are transformed to locations on
+the log-hazard scale when Stan data are created; they are not arithmetic
+prior means.
 }
 }
 \if{html}{\out{<hr>}}
 \if{html}{\out{<a id="method-MultistateModel-set_prior_mean_h0"></a>}}
 \if{latex}{\out{\hypertarget{method-MultistateModel-set_prior_mean_h0}{}}}
 \subsection{Method \code{set_prior_mean_h0()}}{
-Set assumed prior mean baseline hazard rates (side
-effect).
+Set baseline-hazard prior reference rates (side effect).
 \subsection{Usage}{
 \if{html}{\out{<div class="r">}}\preformatted{MultistateModel$set_prior_mean_h0(mean_h0)}\if{html}{\out{</div>}}
 }
@@ -148,7 +150,9 @@ effect).
 \subsection{Arguments}{
 \if{html}{\out{<div class="arguments">}}
 \describe{
-\item{\code{mean_h0}}{Numeric vector with length equal to number of transitions}
+\item{\code{mean_h0}}{A positive numeric vector of natural-scale reference rates
+with length equal to the number of transitions. The argument name is
+retained for compatibility; the values are not arithmetic prior means.}
 }
 \if{html}{\out{</div>}}
 }
@@ -157,8 +161,8 @@ effect).
 \if{html}{\out{<a id="method-MultistateModel-set_prior_mean_h0_data"></a>}}
 \if{latex}{\out{\hypertarget{method-MultistateModel-set_prior_mean_h0_data}{}}}
 \subsection{Method \code{set_prior_mean_h0_data()}}{
-Set assumed prior mean baseline hazard rates (side
-effect) based on average hazards in data.
+Set baseline-hazard prior reference rates (side effect)
+using event counts and at-risk time pooled within transition type.
 \subsection{Usage}{
 \if{html}{\out{<div class="r">}}\preformatted{MultistateModel$set_prior_mean_h0_data(data)}\if{html}{\out{</div>}}
 }
diff --git a/man/MultistateSystem.Rd b/man/MultistateSystem.Rd
index 9e550af..6d9931e 100644
--- a/man/MultistateSystem.Rd
+++ b/man/MultistateSystem.Rd
@@ -243,7 +243,8 @@ Evaluate log instant hazard
 
 \item{\code{w}}{Spline basis function weights (vector)}
 
-\item{\code{log_w0}}{Intercept (log)}
+\item{\code{log_w0}}{Intercept (log), or \code{-Inf} for an identically zero
+baseline hazard}
 
 \item{\code{log_m}}{Hazard multiplier (log)}
 
@@ -295,7 +296,9 @@ Boolean value
 \if{html}{\out{<a id="method-MultistateSystem-max_inst_hazard"></a>}}
 \if{latex}{\out{\hypertarget{method-MultistateSystem-max_inst_hazard}{}}}
 \subsection{Method \code{max_inst_hazard()}}{
-Max instant hazard on interval (t1, t2)
+Approximate maximum instantaneous hazard on an interval,
+evaluated on a 100-point grid. Path generation inflates this value by
+five percent before thinning.
 \subsection{Usage}{
 \if{html}{\out{<div class="r">}}\preformatted{MultistateSystem$max_inst_hazard(t1, t2, w, log_w0, log_m)}\if{html}{\out{</div>}}
 }
@@ -309,7 +312,8 @@ Max instant hazard on interval (t1, t2)
 
 \item{\code{w}}{Spline basis function weights (vector)}
 
-\item{\code{log_w0}}{Intercept (log)}
+\item{\code{log_w0}}{Intercept (log), or \code{-Inf} for an identically zero
+baseline hazard}
 
 \item{\code{log_m}}{Hazard multiplier (log)}
 }
@@ -320,7 +324,9 @@ Max instant hazard on interval (t1, t2)
 \if{html}{\out{<a id="method-MultistateSystem-simulate"></a>}}
 \if{latex}{\out{\hypertarget{method-MultistateSystem-simulate}{}}}
 \subsection{Method \code{simulate()}}{
-Generate paths
+Generate paths. For computational efficiency, a transition
+whose five-percent-inflated numerical hazard envelope is at most
+\code{1e-9} is treated as impossible.
 \subsection{Usage}{
 \if{html}{\out{<div class="r">}}\preformatted{MultistateSystem$simulate(
   w,
@@ -354,7 +360,8 @@ time of the model is used.}
 
 \item{\code{n_rep}}{Number of repetitions to do for each draw.}
 
-\item{\code{min_t_step}}{Minimal time step.}
+\item{\code{min_t_step}}{Minimum separation applied to consecutive simulated
+event times when this does not move an event to or beyond \code{t_max}.}
 }
 \if{html}{\out{</div>}}
 }
diff --git a/man/fit_stan.Rd b/man/fit_stan.Rd
index 8ba7968..1545611 100644
--- a/man/fit_stan.Rd
+++ b/man/fit_stan.Rd
@@ -23,8 +23,8 @@ fit_stan(
 
 \item{set_normalizers}{Set covariate normalization automatically?}
 
-\item{set_prior_h0}{Set prior mean average baseline hazard levels based
-on data?}
+\item{set_prior_h0}{Set baseline-hazard prior reference rates from event
+counts and at-risk time in the data?}
 
 \item{filepath}{Deprecated.}
 
@@ -45,7 +45,7 @@ A \code{\link{MultistateModelFit}} object.
 }
 \description{
 \emph{NOTE:} This function has a side effect of setting covariate
-normalizers, prior assumed mean baseline hazard, and concentration upper
+normalizers, baseline-hazard prior reference rates, and concentration upper
 bound (PK) based on data.
 }
 \seealso{
diff --git a/tests/testthat/test-correctness-regressions.R b/tests/testthat/test-correctness-regressions.R
new file mode 100644
index 0000000..b17aed8
--- /dev/null
+++ b/tests/testthat/test-correctness-regressions.R
@@ -0,0 +1,388 @@
+make_sparse_survival_paths <- function(event = TRUE) {
+  tm <- transmat_survival(c("Healthy", "Dead"))
+  subject_df <- data.frame(subject_id = c("a", "b"))
+  path_df <- data.frame(
+    path_id = c(1L, 1L, 2L, 2L),
+    state = c(1L, if (event) 2L else 1L, 1L, 1L),
+    time = c(0, 1, 0, 2),
+    trans_idx = c(0L, if (event) 1L else 0L, 0L, 0L)
+  )
+  link_df <- data.frame(
+    path_id = 1:2,
+    subject_id = c("a", "b"),
+    draw_idx = 1L,
+    rep_idx = 1L
+  )
+  PathData$new(subject_df, path_df, link_df, tm)
+}
+
+test_that("baseline-hazard prior uses log scale in Stan data", {
+  pd <- make_sparse_survival_paths()
+  mod <- create_msm(pd$transmat, t_max = 2, n_grid = 10)
+  mod$set_prior_mean_h0_data(pd)
+
+  expect_equal(mod$get_prior_mean_h0(), 1 / 3)
+  expect_equal(create_stan_data_model(mod)$mu_log_w0, log(1 / 3))
+  expect_error(mod$set_prior_mean_h0(0), "strictly positive")
+})
+
+test_that("a modeled transition with no events is rejected", {
+  pd <- make_sparse_survival_paths(event = FALSE)
+  mod <- create_msm(pd$transmat, t_max = 2, n_grid = 10)
+
+  expect_error(
+    prefit_checks(mod, JointData$new(pd, NULL)),
+    "no observed events: Healthy -> Dead"
+  )
+  expect_error(
+    mod$set_prior_mean_h0_data(pd),
+    "no observed events: Healthy -> Dead"
+  )
+  mod$prior_only <- TRUE
+  expect_error(
+    prefit_checks(mod, JointData$new(pd, NULL)),
+    "no observed events: Healthy -> Dead"
+  )
+  mod$set_prior_mean_h0(1 / 3)
+  sd <- create_stan_data(mod, pd)
+  expect_equal(sd$D_trans, 0)
+  expect_equal(sd$sum_trans, 0)
+  expect_equal(sd$sum_risk, 2)
+})
+
+test_that("empirical hazard calibration pools within transition type", {
+  tm <- transmat_illnessdeath(c("Healthy", "Ill", "Dead"))
+  subject_df <- data.frame(subject_id = c("a", "b", "c"))
+  path_df <- data.frame(
+    path_id = c(1L, 1L, 1L, 2L, 2L, 3L, 3L),
+    state = c(1L, 2L, 3L, 1L, 3L, 1L, 1L),
+    time = c(0, 1, 3, 0, 2, 0, 4),
+    trans_idx = c(0L, 1L, 3L, 0L, 2L, 0L, 0L)
+  )
+  link_df <- data.frame(
+    path_id = 1:3,
+    subject_id = c("a", "b", "c"),
+    draw_idx = 1L,
+    rep_idx = 1L
+  )
+  pd <- PathData$new(subject_df, path_df, link_df, tm)
+  rates <- average_haz_per_ttype(pd) |>
+    dplyr::arrange(.data$trans_idx)
+
+  expect_equal(exp(rates$log_h0_avg), c(1 / 7, 2 / 9, 2 / 9))
+})
+
+test_that("fit validation rejects unsupported transitions but prediction data do not", {
+  tm <- TransitionMatrix$new(
+    matrix(c(
+      0, 1, 0, 0,
+      0, 0, 0, 0,
+      0, 0, 0, 1,
+      0, 0, 0, 0
+    ), nrow = 4, byrow = TRUE),
+    c("A", "B", "C", "D")
+  )
+  subject_df <- data.frame(subject_id = c("a", "b"))
+  path_df <- data.frame(
+    path_id = c(1L, 1L, 2L, 2L),
+    state = c(1L, 2L, 1L, 1L),
+    time = c(0, 1, 0, 2),
+    trans_idx = c(0L, 1L, 0L, 0L)
+  )
+  link_df <- data.frame(
+    path_id = 1:2,
+    subject_id = c("a", "b"),
+    draw_idx = 1L,
+    rep_idx = 1L
+  )
+  pd <- PathData$new(subject_df, path_df, link_df, tm)
+  mod <- create_msm(tm, t_max = 2, n_grid = 10)
+
+  expect_error(
+    prefit_checks(mod, JointData$new(pd, NULL)),
+    "no observed events: C -> D"
+  )
+  expect_error(
+    average_haz_per_ttype(pd),
+    "no observed events: C -> D"
+  )
+
+  mod$set_prior_mean_h0(c(1 / 3, 1 / 3))
+  sd <- create_stan_data(mod, pd)
+  expect_equal(sd$sum_trans, c(1, 0))
+  expect_equal(sd$sum_risk, c(2, 0))
+})
+
+test_that("zero-duration rows do not count as hazard information", {
+  tm <- transmat_survival(c("Healthy", "Dead"))
+  subject_df <- data.frame(subject_id = "a")
+  link_df <- data.frame(
+    path_id = 1L, subject_id = "a", draw_idx = 1L, rep_idx = 1L
+  )
+  make_pd <- function(event) {
+    PathData$new(
+      subject_df,
+      data.frame(
+        path_id = c(1L, 1L),
+        state = c(1L, if (event) 2L else 1L),
+        time = c(0, 0),
+        trans_idx = c(0L, if (event) 1L else 0L)
+      ),
+      link_df,
+      tm
+    )
+  }
+
+  expect_error(
+    create_stan_data_transitions(make_pd(TRUE)),
+    "observed transition has no positive at-risk time"
+  )
+  sd <- create_stan_data_transitions(make_pd(FALSE))
+  expect_equal(sd$D_trans, 0)
+  expect_equal(sd$sum_trans, 0)
+})
+
+test_that("state-visit probabilities retain zero-event subjects", {
+  pd <- make_sparse_survival_paths()
+  p <- p_state_visit(pd, by = "subject_id") |>
+    dplyr::arrange(.data$subject_id)
+
+  expect_equal(p$subject_id, c("a", "b"))
+  expect_equal(p$prob, c(1, 0))
+  expect_equal(p$n_event, c(1, 0))
+})
+
+test_that("posterior hazard draws remain draw-major", {
+  w <- array(1:3, dim = c(3, 1, 1))
+  log_w0 <- matrix(11:13, nrow = 3)
+  out <- repeat_hazard_parameter_draws(w, log_w0, n_subjects = 2)
+
+  expect_equal(as.vector(out$w), c(1, 1, 2, 2, 3, 3))
+  expect_equal(as.vector(out$log_w0), c(11, 11, 12, 12, 13, 13))
+})
+
+test_that("noninteger time grids contain exactly the requested midpoints", {
+  mod <- create_msm(transmat_survival(), t_max = 365.25, n_grid = 1000)
+  grid <- create_stan_data_timegrid(mod)
+
+  expect_length(grid$t_grid, 1000)
+  expect_equal(grid$N_grid, 1000)
+  expect_true(all(grid$t_grid > 0 & grid$t_grid < mod$get_tmax()))
+  expect_equal(diff(range(grid$t_grid)), 999 * grid$delta_grid)
+})
+
+test_that("transition probabilities honor a delayed start time", {
+  sys <- MultistateSystem$new(transmat_survival())
+  sys$set_knots(c(0, 10))
+  lambda <- 0.2
+  p <- solve_trans_prob_matrix(
+    sys, t_out = c(6, 7), log_w0 = log(lambda), t_start = 5
+  )
+
+  expect_equal(p[, 1, 1], exp(-lambda * c(1, 2)), tolerance = 1e-6)
+  expect_equal(p[, 1, 2], 1 - exp(-lambda * c(1, 2)), tolerance = 1e-6)
+  expect_equal(
+    solve_trans_prob_matrix(sys, t_out = 5, log_w0 = log(lambda), t_start = 5)[1, , ],
+    diag(2)
+  )
+  expect_error(
+    solve_trans_prob_matrix(sys, t_out = c(7, 6), log_w0 = log(lambda), t_start = 5),
+    "strictly increasing"
+  )
+})
+
+test_that("identity transition probabilities still validate hazard parameters", {
+  sys <- MultistateSystem$new(transmat_survival())
+  sys$set_knots(c(0, 10))
+  log_lambda <- log(0.2)
+
+  expect_error(
+    solve_trans_prob_matrix(sys, 5, rep(log_lambda, 2), t_start = 5)
+  )
+  expect_error(
+    solve_trans_prob_matrix(
+      sys,
+      5,
+      log_lambda,
+      w = matrix(0, nrow = 2, ncol = sys$num_weights()),
+      t_start = 5
+    )
+  )
+  expect_error(
+    solve_trans_prob_matrix(sys, 5, log_lambda, log_m = c(0, 0), t_start = 5)
+  )
+})
+
+test_that("thinning uses the 100-point hazard envelope", {
+  sys <- MultistateSystem$new(transmat_survival())
+  sys$set_knots(c(0, 0.4999, 0.5001, 1))
+  w <- c(-3, 3, -3, 3, -3)
+  bound <- sys$max_inst_hazard(0, 1, w, log_w0 = -2, log_m = 0.4)
+  grid <- seq(0, 1, length.out = 100)
+  evaluated <- exp(sys$log_inst_hazard(grid, w, -2, 0.4))
+
+  expect_equal(bound, max(evaluated))
+})
+
+test_that("path generation skips hazards below the deliberate cutoff", {
+  sys <- MultistateSystem$new(transmat_survival())
+  sys$set_knots(c(0, 1e9))
+  w <- array(0, dim = c(1, 1, sys$num_weights()))
+  path <- suppressMessages(sys$simulate(
+    w = w,
+    log_w0 = matrix(log(0.5e-9), nrow = 1),
+    log_m = matrix(0, nrow = 1),
+    t_max = 1e9
+  ))
+
+  expect_equal(path$trans_idx, c(0, 0))
+})
+
+test_that("zero baseline hazard has zero bound and generates no event", {
+  sys <- MultistateSystem$new(transmat_survival())
+  sys$set_knots(c(0, 1))
+  w <- array(0, dim = c(1, 1, sys$num_weights()))
+
+  expect_equal(sys$max_inst_hazard(0, 1, w[1, 1, ], -Inf, 0), 0)
+
+  path <- suppressMessages(sys$simulate(
+    w = w,
+    log_w0 = matrix(-Inf, nrow = 1),
+    log_m = matrix(0, nrow = 1),
+    t_max = 1
+  ))
+  expect_equal(path$time, c(0, 1))
+  expect_false(any(path$is_event == 1))
+})
+
+test_that("path simulation respects its horizon and censor-row contract", {
+  sys <- MultistateSystem$new(transmat_survival())
+  sys$set_knots(c(0, 1))
+  w <- array(0, dim = c(1, 1, sys$num_weights()))
+  log_w0 <- matrix(log(100), nrow = 1)
+  log_m <- matrix(0, nrow = 1)
+
+  set.seed(1)
+  path <- suppressMessages(sys$simulate(
+    w,
+    log_w0,
+    log_m,
+    t_start = 0,
+    t_max = 0.5,
+    min_t_step = 1
+  ))
+  expect_lte(max(path$time), 0.5)
+  expect_equal(tail(path$time, 1), 0.5)
+  expect_equal(tail(path$is_event, 1), 0)
+  expect_equal(tail(path$is_censor, 1), 1)
+  expect_equal(tail(path$trans_idx, 1), 0)
+
+  zero_horizon <- suppressMessages(sys$simulate(
+    w,
+    log_w0,
+    log_m,
+    t_start = 0.5,
+    t_max = 0.5
+  ))
+  expect_equal(nrow(zero_horizon), 1)
+  expect_equal(zero_horizon$time, 0.5)
+  expect_equal(zero_horizon$state, 1)
+  expect_equal(zero_horizon$is_event, 0)
+  expect_equal(zero_horizon$is_censor, 1)
+  expect_equal(zero_horizon$trans_idx, 0)
+
+  expect_error(
+    suppressMessages(sys$simulate(
+      w,
+      log_w0,
+      log_m,
+      t_start = 0.6,
+      t_max = 0.5
+    ))
+  )
+})
+
+test_that("steady-state PK solution is finite when ka equals ke", {
+  pk <- PKModel$new(list(ka = NULL, CL = NULL, V2 = NULL))
+  times <- c(0, 1, 6, 12)
+  theta_equal <- list(ka = 0.2, CL = 0.2, V2 = 1)
+  theta_near <- list(ka = 0.2, CL = 0.2 * (1 + 1e-10), V2 = 1)
+  x <- pk$simulate_ss(times, theta_equal, dose = 30, tau = 12)
+  y <- pk$simulate_ss(times, theta_near, dose = 30, tau = 12)
+
+  expect_true(all(is.finite(x) & x >= 0))
+  expect_equal(x, y, tolerance = 1e-7)
+})
+
+test_that("near-equal steady-state PK limit retains the first-order ka term", {
+  pk <- PKModel$new(list(ka = NULL, CL = NULL, V2 = NULL))
+  times <- c(0, 1, 6, 12)
+  dose <- 30
+  tau <- 12
+  ke <- 0.2
+
+  for (sign in c(-1, 1)) {
+    ka <- ke * (1 + sign * 5e-9)
+    k <- 0.5 * (ka + ke)
+    tt <- times %% tau
+    f <- exp(-k * tt) / (-expm1(-k * tau))
+    expected <- dose * ka * f * (tt + tau / expm1(k * tau))
+    actual <- pk$simulate_ss(
+      times,
+      list(ka = ka, CL = ke, V2 = 1),
+      dose = dose,
+      tau = tau
+    )
+    expect_equal(actual, expected, tolerance = 1e-13)
+  }
+})
+
+test_that("small but distinct PK rates are not treated as equal", {
+  pk <- PKModel$new(list(ka = NULL, CL = NULL, V2 = NULL))
+  times <- c(0, 1, 6, 12)
+  dose <- 30
+  tau <- 12
+  ka <- 1e-9
+  ke <- 2e-9
+  tt <- times %% tau
+  f_ke <- exp(-ke * tt) / (-expm1(-ke * tau))
+  f_ka <- exp(-ka * tt) / (-expm1(-ka * tau))
+  expected <- dose * ka / (ka - ke) * (f_ke - f_ka)
+  actual <- pk$simulate_ss(
+    times,
+    list(ka = ka, CL = ke, V2 = 1),
+    dose = dose,
+    tau = tau
+  )
+
+  expect_equal(actual, expected, tolerance = 1e-8)
+})
+
+test_that("requested PK interval width uses the correct quantiles", {
+  sim <- data.frame(
+    subject_id = "a",
+    time = 1,
+    val = seq_len(100)
+  )
+  out <- pksim_to_quantiles(sim, ci_alpha = 0.8)
+  expect_equal(out$lower, unname(stats::quantile(sim$val, 0.1)))
+  expect_equal(out$upper, unname(stats::quantile(sim$val, 0.9)))
+})
+
+test_that("base data frames work without covariates", {
+  df <- data.frame(
+    subject_id = rep(c("a", "b"), each = 2),
+    time = c(0, 1, 0, 2),
+    state = c(1L, 2L, 1L, 1L),
+    is_transition = c(FALSE, TRUE, FALSE, FALSE)
+  )
+  pd <- df_to_pathdata(df, transmat_survival())
+  expect_s3_class(pd, "PathData")
+  expect_equal(pd$n_paths(), 2)
+})
+
+test_that("exposure normalization scale must be positive", {
+  mod <- create_msm(transmat_survival())
+  expect_error(mod$set_xpsr_normalizers(0, 0), "strictly positive")
+  expect_message(mod$set_xpsr_normalizers(-1, 0.5), "loc = -1")
+})
diff --git a/tests/testthat/test-stan.R b/tests/testthat/test-stan.R
index 569e603..fcbe9eb 100644
--- a/tests/testthat/test-stan.R
+++ b/tests/testthat/test-stan.R
@@ -69,3 +69,74 @@ test_that("fitting with Stan works (multi-transition)", {
   pp <- plot_state_occupancy(r)
   expect_true(is_ggplot(pp))
 })
+
+test_that("transient PK solution is stable when elimination is much faster", {
+  ensure_exposed_stan_functions()
+  ka <- 0.1
+  ke <- 100
+  t <- 10
+  A0 <- 1
+  expected <- A0 * ka / (ka - ke) * (exp(-ke * t) - exp(-ka * t))
+  actual <- two_cpt_central(t, ka, ke, 1, A0, 0)
+
+  expect_true(is.finite(actual))
+  expect_equal(actual, expected, tolerance = 1e-12)
+})
+
+test_that("Stan PK equal-rate guards match the first-order continuous limits", {
+  ensure_exposed_stan_functions()
+  ke <- 0.2
+  t <- 6
+  tau <- 12
+  dose <- 30
+  A0 <- 2
+  C0 <- 0.4
+
+  for (sign in c(-1, 1)) {
+    ka <- ke * (1 + sign * 5e-9)
+    k <- 0.5 * (ka + ke)
+
+    expected_transient <-
+      A0 * ka * t * exp(-k * t) + C0 * exp(-ke * t)
+    expect_equal(
+      two_cpt_central(t, ka, ke, 1, A0, C0),
+      expected_transient,
+      tolerance = 1e-13
+    )
+
+    tt <- t %% tau
+    f <- exp(-k * tt) / (-expm1(-k * tau))
+    expected_ss <- dose * ka * f * (tt + tau / expm1(k * tau))
+    expect_equal(
+      two_cpt_central_ss(t, tau, dose, ka, ke, 1),
+      expected_ss,
+      tolerance = 1e-13
+    )
+  }
+})
+
+test_that("subject-level Stan hazard multipliers match direct calculation", {
+  ensure_exposed_stan_functions()
+  x_haz <- list(c(-1, 0.5, 2), c(0.25, -0.5, 1))
+  beta_oth <- list(c(0.4, -0.2), c(-0.1, 0.3))
+  x_xpsr <- list(c(1, 2, 4))
+  beta_xpsr <- list(c(0.2, -0.4))
+  loc <- 2
+  scale <- 0.5
+
+  expected <-
+    x_haz[[1]] %o% beta_oth[[1]] +
+    x_haz[[2]] %o% beta_oth[[2]] +
+    ((x_xpsr[[1]] - loc) / scale) %o% beta_xpsr[[1]]
+  actual <- compute_log_hazard_multiplier(
+    3,
+    2,
+    beta_oth,
+    beta_xpsr,
+    x_haz,
+    x_xpsr,
+    loc,
+    scale
+  )
+  expect_equal(actual, expected, tolerance = 1e-14)
+})
diff --git a/vignettes/math.Rmd b/vignettes/math.Rmd
index 1704e82..c5fe4cb 100644
--- a/vignettes/math.Rmd
+++ b/vignettes/math.Rmd
@@ -96,11 +96,23 @@ transitions.
 
 Currently the package uses a standard normal prior for the coefficients 
 $\beta$. The prior for the weights of the spline basis functions is
-set hierarchically so that transitions
-of same type have a shared mean. Also the average log hazard rates
-$b_0^{(h)}$ have a hierarchical prior so that transitions of the same type
-have a shared mean. This prior mean is estimated from the data using a standard
-Cox proportional hazards model fit.
+set hierarchically so that transitions of the same type have a shared mean.
+Because the B-spline basis partitions unity, the likelihood alone cannot
+distinguish a common shift of the spline weights from an opposite shift of the
+log baseline-hazard intercept. The proper priors regularize this decomposition,
+but inference about baseline hazard should be made from the complete hazard
+curve rather than either component in isolation. For transition $h$, the
+intercept prior is centered at $\log(\widehat r_h)$, where
+$\widehat r_h$ is calculated by pooling event counts and at-risk time within
+transition type. Every transition retained in the fitted transition matrix
+must have at least one observed occurrence and positive at-risk time; fitting
+stops otherwise. Consequently, $\widehat r_h$ is a natural-scale prior
+reference rate (a geometric location or conditional median), not an arithmetic
+prior mean.
+Indeed, under the current $\sigma_0 \sim \operatorname{HalfNormal}(0,3)$
+hyperprior, the marginal arithmetic mean of the baseline hazard is not finite;
+prior-predictive checks and sensitivity to this deliberately broad scale are
+therefore important.
 
 You can view the prior in the Stan code for example like this.
 
@@ -510,7 +522,3 @@ pd$paths$prop_matrix()
 
 
 # References
-
-
-
-