This function is used within epidist() to create any custom Stan code which
is injected into brms via the stanvars argument. It is unlikely that
as a user you will need this function, but we export it nonetheless to be
transparent about what exactly is happening inside of a call to epidist().
See also
Other stan:
epidist_stancode.default()
Examples
data <- sierra_leone_ebola_data |>
as_epidist_linelist_data(
pdate_lwr = "date_of_symptom_onset",
sdate_lwr = "date_of_sample_tested"
) |>
as_epidist_aggregate_data() |>
as_epidist_marginal_model()
#> ℹ No primary event upper bound provided, using the primary event lower bound + 1 day as the assumed upper bound.
#> ℹ No secondary event upper bound provided, using the secondary event lower bound + 1 day as the assumed upper bound.
#> ℹ No observation time column provided, using 2015-09-14 as the observation date (the maximum of the secondary event upper bound).
#> ! Setting 2394 relative observation times (`relative_obs_time`) greater than 98
#> (2x the maximum delay) to Inf.
#> ℹ This improves model efficiency by reducing the number of unique observation
#> times in the data.
#> ℹ The impact on model accuracy should be negligible because these relative
#> observation times are high enough to cause very limited right truncation.
#> ℹ The original relative observation times are available in
#> `orig_relative_obs_time`.
#> ℹ Raise `obs_time_threshold` to avoid this behaviour.
family <- epidist_family(data, family = lognormal())
formula <- epidist_formula(data, family = family, formula = mu ~ 1)
#> Warning: Found infinite values in the data, which may cause issues for Stan.
epidist_stancode(data, family = family, formula = formula)
#> [[1]]
#> [[1]]$name
#> [1] ""
#>
#> [[1]]$sdata
#> NULL
#>
#> [[1]]$scode
#> [1] "// code chunks used from epidist 1.0.0\n"
#>
#> [[1]]$block
#> [1] "functions"
#>
#> [[1]]$position
#> [1] "start"
#>
#> [[1]]$pll_args
#> character(0)
#>
#>
#> [[2]]
#> [[2]]$name
#> [1] ""
#>
#> [[2]]$sdata
#> NULL
#>
#> [[2]]$scode
#> [1] "/**\n * Compute the log probability mass function for a marginal model with\n * censoring\n *\n * This function is designed to be read into R where:\n * - 'lognormal' is replaced with the target distribution (e.g., 'lognormal')\n * - 'real mu, real sigma' is replaced with multiple distribution parameters in the format\n * \"real paramname1, real paramname2, ...\".\n * - 'mu, sigma' is replaced with the same parameters as real mu, real sigma but\n * reparameterised according to the brms parameterisation for Stan.\n *\n * @param y Real value of observed delay\n * @param real mu, real sigma Distribution parameters (replaced via regex)\n * @param relative_obs_t Observation time relative to primary window start.\n * Passed as D (the relative observation time) to primarycensored.\n * @param pwindow_width Primary window width (actual time scale)\n * @param swindow_width Secondary window width (actual time scale)\n * @param y_upper Upper bound of delay interval\n * @param delay_min Minimum delay. Passed as L (the left truncation point)\n * to primarycensored.\n * @param primary_params Array of parameters for primary distribution\n *\n * @return Log probability mass with censoring adjustment for marginal model\n */\n real marginal_lognormal_lpmf(data int y, real mu, real sigma, data real relative_obs_t,\n data real pwindow_width, data real swindow_width,\n data real y_upper, data real delay_min,\n array[] real primary_params) {\n\n return primarycensored_lpmf(\n y | 1, {mu, sigma}, pwindow_width, y_upper,\n delay_min, relative_obs_t, 1, primary_params\n );\n}"
#>
#> [[2]]$block
#> [1] "functions"
#>
#> [[2]]$position
#> [1] "start"
#>
#> [[2]]$pll_args
#> character(0)
#>
#>
#> [[3]]
#> [[3]]$name
#> [1] ""
#>
#> [[3]]$sdata
#> NULL
#>
#> [[3]]$scode
#> [1] "// Stan functions from primarycensored version 1.6.0\nreal expgrowth_cdf(real x, real xmin, real xmax, real r) {\n if (x < xmin) {\n return 0;\n }\n if (x > xmax) {\n return 1;\n }\n if (abs(r) < 1e-10) {\n return (x - xmin) / (xmax - xmin);\n }\n return (exp(r * x) - exp(r * xmin)) / (exp(r * xmax) - exp(r * xmin));\n}\nreal expgrowth_lcdf(real x, real xmin, real xmax, real r) {\n if (x < xmin) {\n return negative_infinity();\n }\n if (x > xmax) {\n return 0;\n }\n return log(expgrowth_cdf(x | xmin, xmax, r));\n}\nreal expgrowth_lpdf(real x, real xmin, real xmax, real r) {\n if (x < xmin || x > xmax) {\n return negative_infinity();\n }\n if (abs(r) < 1e-10) {\n return -log(xmax - xmin);\n }\n return log(abs(r)) + r * x -\n log(abs(exp(r * xmax) - exp(r * xmin)));\n}\nvector primary_lcdf_vec(vector p, int primary_id,\n array[] real primary_params, data real pwindow) {\n int N = num_elements(p);\n vector[N] out;\n for (i in 1:N) {\n out[i] = primary_lcdf(p[i] | primary_id, primary_params, pwindow);\n }\n return out;\n}\nreal discretestep_lcdf(\n data real d, vector boundaries, vector pmf,\n int primary_id, array[] real primary_params, data real pwindow\n) {\n int K = num_elements(pmf);\n // Integration support in u = d - p for p in [0, pwindow]. It is not\n // clipped at 0 so boundaries that start below zero (delays with negative\n // support) are handled; for non-negative boundaries the per-bin clip to\n // [boundaries[k], boundaries[k + 1]] below gives the same result.\n real u_min = d - pwindow;\n real u_max = d;\n\n // Structural-zero short-circuit. Below `boundaries[2]` F_step is zero\n // and the bin-1 contribution carries `cum_before = 0`, so the integral\n // collapses to 0. Returning `negative_infinity()` directly keeps\n // `log(0)` off the autodiff tape so downstream `log_diff_exp(a, -inf)`\n // evaluates cleanly with a zero gradient w.r.t. `pmf`.\n if (u_max <= boundaries[2]) return negative_infinity();\n\n // Sub-interval endpoints in u-space, clipped to [u_min, u_max].\n vector[K] lo = fmax(u_min, head(boundaries, K));\n vector[K] hi = fmin(u_max, tail(boundaries, K));\n\n // F_step is right-continuous and on [b_k, b_{k+1}) takes the value\n // sum_{j < k} pmf[j] (mass before bin k). cumulative_sum(pmf) gives\n // the mass through and including bin k, so we shift right by one.\n vector[K] cum_before;\n cum_before[1] = 0;\n if (K > 1) cum_before[2:K] = head(cumulative_sum(pmf), K - 1);\n\n // 0/1 mask drops bins with `hi <= lo` from the reduction without a\n // branch in the inner expression. Built on `data`-level inputs.\n vector[K] active;\n for (k in 1:K) active[k] = hi[k] > lo[k] ? 1 : 0;\n\n // F_primary at lo/hi via two vectorised calls; one masked subtraction\n // gives the per-bin difference for the dot product.\n vector[K] f_lo = primary_lcdf_vec(d - lo, primary_id, primary_params,\n pwindow);\n vector[K] f_hi = primary_lcdf_vec(d - hi, primary_id, primary_params,\n pwindow);\n vector[K] f_diff = (exp(f_lo) - exp(f_hi)) .* active;\n\n real integral = dot_product(cum_before, f_diff);\n\n // Tail region [boundaries[K+1], u_max]: F_step = 1, contributing\n // F_primary(d - tail_start) - F_primary(d - u_max).\n real tail_start = fmax(boundaries[K + 1], u_min);\n if (tail_start < u_max) {\n real fp_tail = exp(primary_lcdf(d - tail_start | primary_id,\n primary_params, pwindow));\n real fp_end = exp(primary_lcdf(d - u_max | primary_id,\n primary_params, pwindow));\n integral += fp_tail - fp_end;\n }\n\n return log(integral);\n}\nvector hazards_to_pmf(vector hazards) {\n int K = num_elements(hazards);\n vector[K] log_surv;\n log_surv[1] = 0;\n if (K > 1) {\n log_surv[2:K] = cumulative_sum(log1m(hazards[1:(K - 1)]));\n }\n return hazards .* exp(log_surv);\n}\nreal discretehazard_lcdf(\n data real d, vector boundaries, vector hazards,\n int primary_id, array[] real primary_params, data real pwindow\n) {\n return discretestep_lcdf(\n d | boundaries, hazards_to_pmf(hazards), primary_id, primary_params,\n pwindow\n );\n}\nreal pstep_lcdf(real t, vector boundaries, vector pmf) {\n int K = num_elements(pmf);\n if (t < boundaries[2]) return negative_infinity();\n if (t >= boundaries[K + 1]) return 0;\n // Right-continuous CDF with jumps at the right edges\n // boundaries[2], ..., boundaries[K + 1]. F(t) = cum_pmf[k] for\n // t in [boundaries[k + 1], boundaries[k + 2]); equivalently the\n // largest k with boundaries[k + 1] <= t. Boundary-on-jump cases\n // (t == boundaries[k + 1]) advance k, matching R's\n // `findInterval(left.open = FALSE)`.\n int k = 1;\n while (k < K && boundaries[k + 2] <= t) k += 1;\n return log(cumulative_sum(pmf)[k]);\n}\nreal phazard_lcdf(real t, vector boundaries, vector hazards) {\n return pstep_lcdf(t | boundaries, hazards_to_pmf(hazards));\n}\nreal primarycensored_log_normalizer(real log_cdf_D, real log_cdf_L, real L) {\n if (!is_inf(L)) {\n return log_diff_exp(log_cdf_D, log_cdf_L);\n } else {\n return log_cdf_D;\n }\n}\nreal primarycensored_apply_truncation(real log_cdf, real log_cdf_L,\n real log_normalizer, real L) {\n if (!is_inf(L)) {\n return log_diff_exp(log_cdf, log_cdf_L) - log_normalizer;\n } else {\n return log_cdf - log_normalizer;\n }\n}\nvector primarycensored_truncation_bounds(\n data real L, data real D,\n data int dist_id, array[] real params, data real pwindow,\n data int primary_id, array[] real primary_params\n) {\n vector[2] result;\n // Internal lower bound for the un-truncated distribution: 0 lets the\n // `d <= L` early-exit in primarycensored_lcdf return -inf for delays below\n // the natural support of positive-support distributions; -inf disables that\n // short-circuit so distributions with support on the reals are integrated.\n // Expression is inlined (rather than bound to a local) so Stan's data-flow\n // checker recognises it as data-only.\n\n // Get CDF at lower truncation point L\n if (is_inf(L)) {\n result[1] = negative_infinity();\n } else {\n result[1] = primarycensored_lcdf(\n L | dist_id, params, pwindow,\n dist_has_positive_support(dist_id) ? 0.0 : negative_infinity(),\n positive_infinity(), primary_id, primary_params\n );\n }\n\n // Get CDF at upper truncation point D\n if (is_inf(D)) {\n result[2] = 0;\n } else {\n result[2] = primarycensored_lcdf(\n D | dist_id, params, pwindow,\n dist_has_positive_support(dist_id) ? 0.0 : negative_infinity(),\n positive_infinity(), primary_id, primary_params\n );\n }\n\n return result;\n}\nreal primarycensored_cdf(data real d, data int dist_id, array[] real params,\n data real pwindow, data real L, data real D,\n data int primary_id,\n array[] real primary_params) {\n real result;\n if (d <= L) {\n return 0;\n }\n\n if (d >= D) {\n return 1;\n }\n\n // Check if an analytical solution exists\n if (check_for_analytical(dist_id, primary_id)) {\n // Use analytical solution\n result = primarycensored_analytical_cdf(\n d | dist_id, params, pwindow, L, D, primary_id, primary_params\n );\n } else {\n // Use numerical integration for other cases. The integration variable\n // ranges over the primary-event time, so the natural lower bound is\n // d - pwindow. For positive-support delays the integrand `F_delay(t)` is\n // 0 for t <= 0, so an unclipped lower bound just adds a flat zero region\n // for negative t. Distributions with support on the reals also accept the\n // unclipped lower bound directly.\n real lower_bound = d - pwindow;\n int n_params = num_elements(params);\n int n_primary_params = num_elements(primary_params);\n array[n_params + n_primary_params] real theta = append_array(params, primary_params);\n array[4] int ids = {dist_id, primary_id, n_params, n_primary_params};\n\n vector[1] y0 = rep_vector(0.0, 1);\n result = ode_rk45(primarycensored_ode, y0, lower_bound, {d}, theta, {d, pwindow}, ids)[1, 1];\n\n // Apply truncation normalization on log scale for numerical stability.\n // Skip when F(L) = 0 makes it a no-op (positive support, L <= 0).\n if (!is_inf(D) || L > 0 ||\n (!is_inf(L) && !dist_has_positive_support(dist_id))) {\n real log_result = log(result);\n vector[2] bounds = primarycensored_truncation_bounds(\n L, D, dist_id, params, pwindow, primary_id, primary_params\n );\n real log_cdf_L = bounds[1];\n real log_cdf_D = bounds[2];\n\n real log_normalizer = primarycensored_log_normalizer(log_cdf_D, log_cdf_L, L);\n log_result = primarycensored_apply_truncation(\n log_result, log_cdf_L, log_normalizer, L\n );\n result = exp(log_result);\n }\n }\n\n return result;\n}\nreal primarycensored_lcdf(data real d, data int dist_id, array[] real params,\n data real pwindow, data real L, data real D,\n data int primary_id,\n array[] real primary_params) {\n real result;\n\n if (d <= L) {\n return negative_infinity();\n }\n\n if (d >= D) {\n return 0;\n }\n\n // Check if an analytical solution exists. The internal lower bound is 0 for\n // positive-support delays (lets the d <= L early-exit return -inf for d <= 0)\n // and -inf for distributions with support on the reals.\n if (check_for_analytical(dist_id, primary_id)) {\n result = primarycensored_analytical_lcdf(\n d | dist_id, params, pwindow,\n dist_has_positive_support(dist_id) ? 0.0 : negative_infinity(),\n positive_infinity(), primary_id, primary_params\n );\n } else {\n // Use numerical integration\n result = log(primarycensored_cdf(\n d | dist_id, params, pwindow,\n dist_has_positive_support(dist_id) ? 0.0 : negative_infinity(),\n positive_infinity(), primary_id, primary_params\n ));\n }\n\n // Handle truncation normalization. Skip when F(L) = 0 makes it a no-op\n // (positive support, L <= 0) to avoid the cancelling log_diff_exp.\n if (!is_inf(D) || L > 0 ||\n (!is_inf(L) && !dist_has_positive_support(dist_id))) {\n vector[2] bounds = primarycensored_truncation_bounds(\n L, D, dist_id, params, pwindow, primary_id, primary_params\n );\n real log_cdf_L = bounds[1];\n real log_cdf_D = bounds[2];\n\n real log_normalizer = primarycensored_log_normalizer(log_cdf_D, log_cdf_L, L);\n result = primarycensored_apply_truncation(result, log_cdf_L, log_normalizer, L);\n }\n\n return result;\n}\nreal primarycensored_lpmf(data int d, data int dist_id, array[] real params,\n data real pwindow, data real d_upper,\n data real L, data real D, data int primary_id,\n array[] real primary_params) {\n if (d_upper > D) {\n reject(\"Upper truncation point is greater than D. It is \", d_upper,\n \" and D is \", D, \". Resolve this by increasing D to be greater or equal to d + swindow or decreasing swindow.\");\n }\n if (d_upper <= d) {\n reject(\"Upper truncation point is less than or equal to d. It is \", d_upper,\n \" and d is \", d, \". Resolve this by increasing d to be less than d_upper.\");\n }\n if (d < L) {\n return negative_infinity();\n }\n real log_cdf_upper = primarycensored_lcdf(\n d_upper | dist_id, params, pwindow,\n dist_has_positive_support(dist_id) ? 0.0 : negative_infinity(),\n positive_infinity(), primary_id, primary_params\n );\n real log_cdf_lower = primarycensored_lcdf(\n d | dist_id, params, pwindow,\n dist_has_positive_support(dist_id) ? 0.0 : negative_infinity(),\n positive_infinity(), primary_id, primary_params\n );\n\n // Apply truncation normalization: log((F(d_upper) - F(d)) / (F(D) - F(L))).\n // Skip when F(L) = 0 makes it a no-op (positive support, L <= 0).\n if (!is_inf(D) || L > 0 ||\n (!is_inf(L) && !dist_has_positive_support(dist_id))) {\n real log_cdf_D;\n real log_cdf_L;\n\n // Get CDF at lower truncation point L\n if (is_inf(L)) {\n // No left truncation (L = -inf sentinel)\n log_cdf_L = negative_infinity();\n } else if (d == L) {\n // Reuse already computed CDF at d\n log_cdf_L = log_cdf_lower;\n } else {\n // Compute CDF at L directly\n log_cdf_L = primarycensored_lcdf(\n L | dist_id, params, pwindow,\n dist_has_positive_support(dist_id) ? 0.0 : negative_infinity(),\n positive_infinity(), primary_id, primary_params\n );\n }\n\n // Get CDF at upper truncation point D\n if (d_upper == D) {\n log_cdf_D = log_cdf_upper;\n } else if (is_inf(D)) {\n log_cdf_D = 0;\n } else {\n log_cdf_D = primarycensored_lcdf(\n D | dist_id, params, pwindow,\n dist_has_positive_support(dist_id) ? 0.0 : negative_infinity(),\n positive_infinity(), primary_id, primary_params\n );\n }\n\n real log_normalizer = primarycensored_log_normalizer(log_cdf_D, log_cdf_L, L);\n return log_diff_exp(log_cdf_upper, log_cdf_lower) - log_normalizer;\n } else {\n return log_diff_exp(log_cdf_upper, log_cdf_lower);\n }\n}\nint check_for_uniform_terms(int dist_id, int primary_id) {\n if (primary_id != 1) return 0;\n return dist_id == 2 || dist_id == 1 || dist_id == 3 || dist_id == 5;\n}\nint check_for_analytical(int dist_id, int primary_id) {\n // Gamma, Lognormal, Weibull and generalised gamma with a Uniform primary\n if (check_for_uniform_terms(dist_id, primary_id)) return 1;\n // Keep this primary list in sync with `primary_lcdf`; see the note above.\n if (dist_id == 26 || dist_id == 27 || dist_id == 28) {\n return primary_id == 1 || primary_id == 2;\n }\n return 0; // No analytical solution for other combinations\n}\nreal primarycensored_uniform_lcdf_from_terms(vector terms_d, vector terms_q,\n data real pwindow) {\n real log_A = log_sum_exp(terms_d[1], terms_q[2]);\n real log_B = log_sum_exp(terms_q[1], terms_d[2]);\n // Deep enough into the lower tail every term underflows together. Both\n // are then `-inf` and `log_diff_exp` would give NaN, so return the limit\n // directly.\n if (log_A == negative_infinity() && log_B == negative_infinity()) {\n return negative_infinity();\n }\n return log_diff_exp(log_A, log_B) - log(pwindow);\n}\nvector primarycensored_gamma_uniform_terms(real t,\n array[] real params) {\n if (t <= 0) {\n return rep_vector(negative_infinity(), 2);\n }\n real shape = params[1];\n real rate = params[2];\n // log E where E = k * theta = shape / rate is the mean of the delay\n real log_E = log(shape) - log(rate);\n // F_T(t; k) and the recursion to F_T(t; k+1):\n // P(k+1, y) = P(k, y) - y^k e^{-y} / Gamma(k+1), with y = rate * t\n real log_F_T_k = gamma_lcdf(t | shape, rate);\n real gamma_kp1_pdf_log = shape * log(rate * t) - rate * t\n - lgamma(shape + 1);\n real log_F_T_kp1 = log_diff_exp(log_F_T_k, gamma_kp1_pdf_log);\n return [log(t) + log_F_T_k, log_E + log_F_T_kp1]';\n}\nreal primarycensored_gamma_uniform_lcdf(data real d, real q,\n array[] real params,\n data real pwindow) {\n return primarycensored_uniform_lcdf_from_terms(\n primarycensored_gamma_uniform_terms(d, params),\n primarycensored_gamma_uniform_terms(q, params), pwindow\n );\n}\nvector primarycensored_lognormal_uniform_terms(real t,\n array[] real params) {\n real mu = params[1];\n real sigma = params[2];\n real mu_sigma2 = mu + square(sigma);\n // log E where E = exp(mu + sigma^2/2) is the mean of the delay\n real log_E = mu + 0.5 * square(sigma);\n real log_t_F_T = lognormal_lcdf_underflows(t, mu, sigma)\n ? negative_infinity()\n : log(t) + lognormal_lcdf(t | mu, sigma);\n real log_E_tF_T = lognormal_lcdf_underflows(t, mu_sigma2, sigma)\n ? negative_infinity()\n : log_E + lognormal_lcdf(t | mu_sigma2, sigma);\n return [log_t_F_T, log_E_tF_T]';\n}\nreal primarycensored_lognormal_uniform_lcdf(data real d, real q,\n array[] real params,\n data real pwindow) {\n return primarycensored_uniform_lcdf_from_terms(\n primarycensored_lognormal_uniform_terms(d, params),\n primarycensored_lognormal_uniform_terms(q, params), pwindow\n );\n}\nreal log_weibull_g(real t, real shape, real scale) {\n real x = pow(t * inv(scale), shape);\n real a = 1 + inv(shape);\n return log(gamma_p(a, x)) + lgamma(a);\n}\nvector primarycensored_weibull_uniform_terms(real t,\n array[] real params) {\n if (t <= 0) {\n return rep_vector(negative_infinity(), 2);\n }\n real shape = params[1];\n real scale = params[2];\n return [\n log(t) + weibull_lcdf(t | shape, scale),\n log(scale) + log_weibull_g(t, shape, scale)\n ]';\n}\nreal primarycensored_weibull_uniform_lcdf(data real d, real q,\n array[] real params,\n data real pwindow) {\n return primarycensored_uniform_lcdf_from_terms(\n primarycensored_weibull_uniform_terms(d, params),\n primarycensored_weibull_uniform_terms(q, params), pwindow\n );\n}\nvector primarycensored_gengamma_uniform_terms(real t,\n array[] real params) {\n if (t <= 0) {\n return rep_vector(negative_infinity(), 2);\n }\n real shape = params[1];\n real scale = params[2];\n real k = params[3];\n real k_shift = k + inv(shape);\n real log_E = log(scale) + lgamma(k_shift) - lgamma(k);\n return [\n log(t) + gengamma_lcdf(t | shape, scale, k),\n log_E + gengamma_lcdf(t | shape, scale, k_shift)\n ]';\n}\nreal primarycensored_gengamma_uniform_lcdf(data real d, real q,\n array[] real params,\n data real pwindow) {\n return primarycensored_uniform_lcdf_from_terms(\n primarycensored_gengamma_uniform_terms(d, params),\n primarycensored_gengamma_uniform_terms(q, params), pwindow\n );\n}\nreal primarycensored_analytical_lcdf_raw(data real d, int dist_id,\n array[] real params,\n data real pwindow,\n int primary_id,\n array[] real primary_params) {\n real q = max({d - pwindow, 0});\n\n if (dist_id == 2 && primary_id == 1) {\n return primarycensored_gamma_uniform_lcdf(d | q, params, pwindow);\n } else if (dist_id == 1 && primary_id == 1) {\n return primarycensored_lognormal_uniform_lcdf(d | q, params, pwindow);\n } else if (dist_id == 3 && primary_id == 1) {\n return primarycensored_weibull_uniform_lcdf(d | q, params, pwindow);\n } else if (dist_id == 5 && primary_id == 1) {\n return primarycensored_gengamma_uniform_lcdf(d | q, params, pwindow);\n } else if (dist_id == 26) {\n // params = [boundaries (K+1), pmf (K)]; length 2*K + 1.\n int K = (size(params) - 1) %/% 2;\n return discretestep_lcdf(\n d | to_vector(segment(params, 1, K + 1)),\n to_vector(segment(params, K + 2, K)),\n primary_id, primary_params, pwindow\n );\n } else if (dist_id == 27 || dist_id == 28) {\n // params = [boundaries (K+1), hazards (K)]; length 2*K + 1. The last\n // hazard must equal 1. RW (27) and RE (28) only differ in their\n // prior so they share this likelihood dispatch.\n int K = (size(params) - 1) %/% 2;\n return discretehazard_lcdf(\n d | to_vector(segment(params, 1, K + 1)),\n to_vector(segment(params, K + 2, K)),\n primary_id, primary_params, pwindow\n );\n }\n return negative_infinity();\n}\nreal primarycensored_analytical_lcdf(data real d, int dist_id,\n array[] real params,\n data real pwindow, data real L,\n data real D, int primary_id,\n array[] real primary_params) {\n if (d <= L) return negative_infinity();\n if (d >= D) return 0;\n\n real result = primarycensored_analytical_lcdf_raw(\n d, dist_id, params, pwindow, primary_id, primary_params\n );\n\n // Apply truncation normalization\n if (!is_inf(D) || L > 0) {\n vector[2] bounds = primarycensored_truncation_bounds(\n L, D, dist_id, params, pwindow, primary_id, primary_params\n );\n real log_cdf_L = bounds[1];\n real log_cdf_D = bounds[2];\n\n real log_normalizer = primarycensored_log_normalizer(log_cdf_D, log_cdf_L, L);\n result = primarycensored_apply_truncation(result, log_cdf_L, log_normalizer, L);\n }\n\n return result;\n}\nreal primarycensored_analytical_cdf(data real d, int dist_id,\n array[] real params,\n data real pwindow, data real L,\n data real D, int primary_id,\n array[] real primary_params) {\n return exp(primarycensored_analytical_lcdf(d | dist_id, params, pwindow, L, D, primary_id, primary_params));\n}\nint dist_has_positive_support(data int dist_id) {\n if (dist_id == 1) return 1; // Lognormal\n if (dist_id == 2) return 1; // Gamma\n if (dist_id == 3) return 1; // Weibull\n if (dist_id == 4) return 1; // Exponential\n if (dist_id == 5) return 1; // Generalised gamma\n if (dist_id == 9) return 1; // Beta (support on [0, 1])\n if (dist_id == 13) return 1; // Chi-square\n if (dist_id == 16) return 1; // Inverse Gamma\n if (dist_id == 19) return 1; // Inverse Chi-square\n if (dist_id == 21) return 1; // Pareto\n if (dist_id == 22) return 1; // Scaled inverse Chi-square\n return 0;\n}\nreal primary_lcdf(real p, int primary_id, array[] real primary_params,\n data real pwindow) {\n if (primary_id == 1) {\n // Uniform on [0, pwindow]: built-in uniform_lcdf matches the package\n // primary semantics over [0, pwindow].\n if (p <= 0) return negative_infinity();\n if (p >= pwindow) return 0;\n return uniform_lcdf(p | 0, pwindow);\n } else if (primary_id == 2) {\n return expgrowth_lcdf(p | 0, pwindow, primary_params[1]);\n }\n reject(\"primary_lcdf: unsupported primary_id \", primary_id);\n}\nint lognormal_lcdf_underflows(real y, real mu, real sigma) {\n if (y <= 0) {\n return 1;\n }\n return (log(y) - mu) / sigma < -38 ? 1 : 0;\n}\nreal gengamma_lcdf(real y, real shape, real scale, real k) {\n return gamma_lcdf(pow(y / scale, shape) | k, 1);\n}\nreal dist_lcdf(real delay, array[] real params, int dist_id) {\n if (dist_has_positive_support(dist_id) && delay <= 0) {\n return negative_infinity();\n }\n\n // IDs match pcd_distributions$stan_id in R\n // Guarded so a lower-tail underflow cannot put a NaN partial on the tape.\n // The downstream `exp(-inf)` differentiates to 0.\n if (dist_id == 1) {\n return lognormal_lcdf_underflows(delay, params[1], params[2])\n ? negative_infinity()\n : lognormal_lcdf(delay | params[1], params[2]);\n }\n else if (dist_id == 2) return gamma_lcdf(delay | params[1], params[2]);\n else if (dist_id == 3) return weibull_lcdf(delay | params[1], params[2]);\n else if (dist_id == 4) return exponential_lcdf(delay | params[1]);\n else if (dist_id == 5) return gengamma_lcdf(delay | params[1], params[2], params[3]);\n else if (dist_id == 9) return beta_lcdf(delay | params[1], params[2]);\n else if (dist_id == 12) return cauchy_lcdf(delay | params[1], params[2]);\n else if (dist_id == 13) return chi_square_lcdf(delay | params[1]);\n else if (dist_id == 15) return gumbel_lcdf(delay | params[1], params[2]);\n else if (dist_id == 16) return inv_gamma_lcdf(delay | params[1], params[2]);\n else if (dist_id == 17) return logistic_lcdf(delay | params[1], params[2]);\n else if (dist_id == 18) return normal_lcdf(delay | params[1], params[2]);\n else if (dist_id == 19) return inv_chi_square_lcdf(delay | params[1]);\n else if (dist_id == 20) return double_exponential_lcdf(delay | params[1], params[2]);\n else if (dist_id == 21) return pareto_lcdf(delay | params[1], params[2]);\n else if (dist_id == 22) return scaled_inv_chi_square_lcdf(delay | params[1], params[2]);\n else if (dist_id == 23) return student_t_lcdf(delay | params[1], params[2], params[3]);\n else if (dist_id == 24) return uniform_lcdf(delay | params[1], params[2]);\n else if (dist_id == 25) return von_mises_lcdf(delay | params[1], params[2]);\n else if (dist_id == 26) {\n // Non-parametric step: params = [boundaries (K+1), pmf (K)].\n int K = (size(params) - 1) %/% 2;\n return pstep_lcdf(\n delay | to_vector(segment(params, 1, K + 1)),\n to_vector(segment(params, K + 2, K))\n );\n }\n else if (dist_id == 27 || dist_id == 28) {\n // Non-parametric discrete hazard: params = [boundaries (K+1),\n // hazards (K)] with hazards[K] = 1. RW (27) and RE (28) share the\n // same likelihood; they only differ in the prior.\n int K = (size(params) - 1) %/% 2;\n return phazard_lcdf(\n delay | to_vector(segment(params, 1, K + 1)),\n to_vector(segment(params, K + 2, K))\n );\n }\n else reject(\"Invalid distribution identifier: \", dist_id);\n}\nreal primary_lpdf(real x, int primary_id, array[] real params, real xmin, real xmax) {\n // Implement switch for different primary distributions\n if (primary_id == 1) return uniform_lpdf(x | xmin, xmax);\n if (primary_id == 2) return expgrowth_lpdf(x | xmin, xmax, params[1]);\n // Add more primary distributions as needed\n reject(\"Invalid primary distribution identifier\");\n}\nvector primarycensored_ode(real t, vector y, array[] real theta,\n array[] real x_r, array[] int x_i) {\n real d = x_r[1];\n int dist_id = x_i[1];\n int primary_id = x_i[2];\n real pwindow = x_r[2];\n int dist_params_len = x_i[3];\n int primary_params_len = x_i[4];\n\n // Extract distribution parameters\n array[dist_params_len] real params;\n if (dist_params_len) {\n params = theta[1:dist_params_len];\n }\n array[primary_params_len] real primary_params;\n if (primary_params_len) {\n int primary_loc = num_elements(theta);\n primary_params = theta[primary_loc - primary_params_len + 1:primary_loc];\n }\n\n real log_cdf = dist_lcdf(t | params, dist_id);\n real log_primary_pdf = primary_lpdf(d - t | primary_id, primary_params, 0, pwindow);\n\n return rep_vector(exp(log_cdf + log_primary_pdf), 1);\n}"
#>
#> [[3]]$block
#> [1] "functions"
#>
#> [[3]]$position
#> [1] "start"
#>
#> [[3]]$pll_args
#> character(0)
#>
#>
#> [[4]]
#> [[4]]$name
#> [1] ""
#>
#> [[4]]$sdata
#> NULL
#>
#> [[4]]$scode
#> [1] "array[0] real primary_params;"
#>
#> [[4]]$block
#> [1] "parameters"
#>
#> [[4]]$position
#> [1] "start"
#>
#> [[4]]$pll_args
#> character(0)
#>
#>
#> attr(,"class")
#> [1] "stanvars"