1 Introduction
1.1 What are we going to do in this vignette
In this vignette, we’ll explore how to use the primarycensored package in your Stan modelling workflow. We’ll cover the following key points:
- Introduction to Stan and its relevance for our analysis
- Overview of the packages we’ll be using
- How to access and use Stan functions provided by
primarycensored - Methods for integrating these Stan functions into your workflow
If you are instead interested in fitting a delay distribution using primarycensored in R see the vignette("fitting-dists-with-stan") vignette or epidist package (which uses primarycensored under the hood).
1.2 What is Stan and why are we using it
Stan is a probabilistic programming language for statistical inference. It provides a flexible and efficient platform for Bayesian modeling and is widely used in various fields of data science and statistics. In this vignette, we’ll use Stan in conjunction with primarycensored to perform Bayesian inference on censored data.
For more information on Stan:
- Stan’s official website
- Stan documentation
- Stan forums for community support and discussions
2 Using Stan code in primarycensored
primarycensored includes a set of Stan functions that mirror the R functions in primarycensored. Documentation for these functions can be found here. We support a range of approaches to integrate this Stan code into your workflow.
2.1 Checking available Stan functions using pcd_stan_functions()
Aside from reading the documentation it is also possible to list the available Stan functions using a helper function directly in R.
## [1] "expgrowth_pdf"
## [2] "expgrowth_lpdf"
## [3] "expgrowth_cdf"
## [4] "expgrowth_lcdf"
## [5] "expgrowth_rng"
## [6] "pstep_lcdf"
## [7] "hazards_to_pmf"
## [8] "phazard_lcdf"
## [9] "primary_lcdf_vec"
## [10] "discretestep_lcdf"
## [11] "discretehazard_lcdf"
## [12] "check_for_uniform_terms"
## [13] "check_for_analytical"
## [14] "check_for_analytical_params"
## [15] "check_for_analytical_delay"
## [16] "primarycensored_uniform_lcdf_from_terms"
## [17] "primarycensored_gamma_uniform_terms"
## [18] "primarycensored_lognormal_uniform_terms"
## [19] "log_weibull_g"
## [20] "primarycensored_weibull_uniform_terms"
## [21] "primarycensored_gengamma_uniform_terms"
## [22] "primarycensored_loglogistic_lcdf"
## [23] "primarycensored_loglogistic_lccdf"
## [24] "loglogistic_ratio_small"
## [25] "loglogistic_partial_at_one"
## [26] "loglogistic_tau_integral"
## [27] "loglogistic_tail_term"
## [28] "loglogistic_ratio_large"
## [29] "loglogistic_moment_ratio"
## [30] "loglogistic_log_error_scale"
## [31] "loglogistic_window_ill_conditioned"
## [32] "loglogistic_numeric_ode"
## [33] "loglogistic_numeric_lcdf"
## [34] "primarycensored_loglogistic_uniform_terms"
## [35] "primarycensored_uniform_terms"
## [36] "primarycensored_gamma_uniform_lcdf"
## [37] "primarycensored_lognormal_uniform_lcdf"
## [38] "primarycensored_weibull_uniform_lcdf"
## [39] "primarycensored_gengamma_uniform_lcdf"
## [40] "primarycensored_loglogistic_uniform_lcdf"
## [41] "primarycensored_analytical_lcdf_raw"
## [42] "primarycensored_analytical_lcdf"
## [43] "primarycensored_analytical_cdf"
## [44] "check_for_analytical_vectorized"
## [45] "primarycensored_analytical_lcdf_vectorized"
## [46] "primarycensored_log_diff_exp"
## [47] "primarycensored_log_std_normal_cdf"
## [48] "primarycensored_stacy_series"
## [49] "check_for_tilt_transform"
## [50] "check_for_tilt_transform_at"
## [51] "primarycensored_gauss_legendre_nodes"
## [52] "primarycensored_gauss_legendre_log_weights"
## [53] "primarycensored_lognormal_tilt_one_panel"
## [54] "primarycensored_lognormal_tilt_panel"
## [55] "primarycensored_lognormal_tilt_bump"
## [56] "primarycensored_lognormal_tilt_quadrature"
## [57] "primarycensored_lognormal_tilt_series"
## [58] "primarycensored_lognormal_tilt_pair"
## [59] "log_tilt_transform_context"
## [60] "log_tilt_transform_pair_shared"
## [61] "log_tilt_transform_pair"
## [62] "primarycensored_moments_from_partial"
## [63] "primarycensored_gamma_tilt_moments"
## [64] "primarycensored_stacy_tilt_moments"
## [65] "primarycensored_tilt_moments"
## [66] "check_for_exptilt"
## [67] "primarycensored_exptilt_terms"
## [68] "primarycensored_exptilt_form"
## [69] "primarycensored_tail_diff"
## [70] "primarycensored_exptilt_lcdf_from_terms"
## [71] "primarycensored_exptilt_small_window_lcdf_from_terms"
## [72] "primarycensored_exptilt_small_delay_lcdf_from_terms"
## [73] "primarycensored_exptilt_lcdf"
## [74] "gumbel_n_terms"
## [75] "check_for_gumbel"
## [76] "check_for_gumbel_params"
## [77] "gumbel_series_accepted"
## [78] "gumbel_log_sum_exp"
## [79] "gumbel_finite_abs"
## [80] "primarycensored_gumbel_terms"
## [81] "primarycensored_gumbel_lcdf_from_terms"
## [82] "gumbel_numeric_z_limits"
## [83] "gumbel_delay_lcdf"
## [84] "gumbel_numeric_delay"
## [85] "gumbel_numeric_point"
## [86] "primarycensored_gumbel_ode"
## [87] "gumbel_numeric_piece_lcdf"
## [88] "gumbel_numeric_spike_lcdf"
## [89] "gumbel_z_log_integrand"
## [90] "gumbel_numeric_z_half"
## [91] "gumbel_numeric_z_lcdf"
## [92] "primarycensored_gumbel_numeric_lcdf"
## [93] "primarycensored_gumbel_numeric_cdf"
## [94] "primarycensored_gumbel_lcdf"
## [95] "gamma_log_lead_logx"
## [96] "gamma_log_pq_logx"
## [97] "gamma_lcdf_logx"
## [98] "gamma_lcdf_logx_pair"
## [99] "gamma_log_rate"
## [100] "gengamma_lcdf"
## [101] "dist_has_positive_support"
## [102] "lognormal_lcdf_underflows"
## [103] "dist_lcdf"
## [104] "primary_lcdf"
## [105] "primary_lpdf"
## [106] "primarycensored_ode"
## [107] "check_for_tlogis"
## [108] "tlogis_log_bound"
## [109] "tlogis_series_terms"
## [110] "tlogis_weights"
## [111] "tlogis_plan"
## [112] "tlogis_n_pos"
## [113] "tlogis_n_neg"
## [114] "tlogis_is_valid"
## [115] "tlogis_plan_applies"
## [116] "check_for_tlogis_params"
## [117] "primarycensored_tlogis_terms"
## [118] "tlogis_terms_diff"
## [119] "tlogis_weighted_sum"
## [120] "tlogis_part_before_window"
## [121] "tlogis_part_above_location"
## [122] "tlogis_part_below_location"
## [123] "primarycensored_tlogis_lcdf_from_terms"
## [124] "tlogis_is_small_delay"
## [125] "primarycensored_tlogis_small_delay_lcdf"
## [126] "primarycensored_tlogis_lcdf"
## [127] "check_for_tlogis_vectorized"
## [128] "primarycensored_tlogis_lcdf_vectorized"
## [129] "primarycensored_log_normalizer"
## [130] "primarycensored_apply_truncation"
## [131] "primarycensored_truncation_bounds"
## [132] "primarycensored_numeric_cdf"
## [133] "primarycensored_cdf"
## [134] "primarycensored_lcdf"
## [135] "primarycensored_lpmf"
## [136] "primarycensored_pmf"
## [137] "primarycensored_lcdf_vectorized"
## [138] "primarycensored_sone_lpmf_vectorized"
## [139] "primarycensored_sone_pmf_vectorized"
## [140] "tgumbel_log_delta_window"
## [141] "tgumbel_log1m_exp_neg_exp"
## [142] "tgumbel_log_norm"
## [143] "tgumbel_lpdf"
## [144] "tgumbel_lcdf"
## [145] "tgumbel_rng"
## [146] "tlogis_log_mass"
## [147] "tlogis_lpdf"
## [148] "tlogis_lcdf"
## [149] "tlogis_cdf"
## [150] "tlogis_rng"
## [151] "tlogis_spike_times"2.2 Accessing Stan functions
Stan functions are accessed using the pcd_load_stan_functions() function.
This function takes the name of the function as an argument and returns the function as a string.
It can additionally write the functions to a file and wrap them in a functions{} block.
pcd_load_stan_functions("primarycensored_lpmf")## [1] "// Stan functions from primarycensored version 1.6.0\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}"2.2.1 Including dependencies automatically
Many Stan functions in primarycensored depend on other functions.
For example, primarycensored_lpmf calls primarycensored_lcdf, which in turn may call analytical or ODE-based implementations.
When using these functions in your own Stan models, you need all the dependencies to be available.
The dependencies argument automatically resolves and includes all required functions in the correct order (dependencies before the functions that use them):
pcd_load_stan_functions("primarycensored_lpmf", dependencies = TRUE)## [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}\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 || dist_id == 29;\n}\nint check_for_analytical(int dist_id, int primary_id) {\n // Gamma, Lognormal, Weibull, generalised gamma and log-logistic with a\n // 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 || primary_id == 3;\n }\n return 0; // No analytical solution for other combinations\n}\nint check_for_tilt_transform(int dist_id, real xi, array[] real params) {\n if (dist_id == 4) return params[1] - xi > 0;\n if (dist_id == 2) return params[2] - xi > 0;\n if (dist_id == 18 || dist_id == 3 || dist_id == 5) return 1;\n if (dist_id == 1) {\n if (params[2] <= 0) return 0;\n if (xi >= 0) return 1;\n return log(-xi) + 2 * log(params[2]) + params[1] < 690;\n }\n return 0;\n}\nint check_for_exptilt(int dist_id, int primary_id) {\n return primary_id == 2\n && (dist_id == 1 || dist_id == 2 || dist_id == 3 || dist_id == 4\n || dist_id == 5 || dist_id == 18);\n}\nreal primarycensored_log_diff_exp(real a, real b) {\n if (b == negative_infinity()) return a;\n if (a <= b) return negative_infinity();\n return log_diff_exp(a, b);\n}\nint check_for_analytical_params(int dist_id, array[] real params,\n int primary_id, array[] real primary_params,\n data real pwindow) {\n if (check_for_exptilt(dist_id, primary_id)) {\n return check_for_tilt_transform(dist_id, -primary_params[1], params);\n }\n if (check_for_tlogis(dist_id, primary_id)) {\n return check_for_tlogis_params(dist_id, params, primary_params, pwindow);\n }\n if (check_for_gumbel(dist_id, primary_id)) {\n return check_for_gumbel_params(\n dist_id, params, primary_params[1], primary_params[2]\n );\n }\n // The log-logistic partial expectation needs a shape of at least 0.01\n if (dist_id == 29 && primary_id == 1) {\n return params[2] >= 0.01;\n }\n return check_for_analytical(dist_id, primary_id);\n}\nint check_for_tilt_transform_at(int dist_id, real xi, array[] real params,\n data real t, data real pwindow) {\n if (!check_for_tilt_transform(dist_id, xi, params)) return 0;\n if (dist_id == 1 && xi > 0 && t > 0) {\n if (xi * t + 9 * sqrt(xi * t) + 30 > 20000) return 0;\n return xi * t <= 60 || xi * pwindow > 2;\n }\n return 1;\n}\nint check_for_analytical_delay(int dist_id, array[] real params,\n int primary_id, array[] real primary_params,\n data real d, data real pwindow) {\n if (!check_for_analytical_params(dist_id, params, primary_id,\n primary_params, pwindow)) {\n return 0;\n }\n if (check_for_exptilt(dist_id, primary_id)) {\n return check_for_tilt_transform_at(\n dist_id, -primary_params[1], params, d, pwindow\n );\n }\n return 1;\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 real log_rate = gamma_log_rate(rate);\n // log E where E = k * theta = shape / rate is the mean of the delay\n real log_E = log(shape) - log_rate;\n // Both CDFs share one series or fraction evaluation at log(rate * t)\n real log_t = log(t);\n vector[2] log_F_T = gamma_lcdf_logx_pair(log_t + log_rate, shape);\n return [log_t + log_F_T[1], log_E + log_F_T[2]]';\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 // gamma_lcdf(x | a, 1) is log(gamma_p(a, x)), but reverse-mode gamma_p()\n // returns zero gradients for x / a > 10\n // (https://github.com/stan-dev/math/issues/2006). gamma_lcdf() has the\n // same value and computes its own gradients without that cutoff.\n return gamma_lcdf(x | a, 1) + 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 real log_t = log(t);\n real log_x = shape * (log_t - log(scale));\n return [\n log_t + gamma_lcdf_logx(log_x, k),\n log_E + gamma_lcdf_logx(log_x, 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 loglogistic_log_error_scale(real x) {\n if (x == negative_infinity()) return negative_infinity();\n return x + log1p(abs(x));\n}\nint loglogistic_window_ill_conditioned(real log_t_cdf_d, real log_t_cdf_q,\n real d, data real pwindow,\n real log_cdf) {\n if (is_nan(log_cdf)) return 1;\n if (log_t_cdf_d == negative_infinity()) return 0;\n real log_error = log(4e-16) - log(pwindow) + log_sum_exp(\n loglogistic_log_error_scale(log_t_cdf_d),\n loglogistic_log_error_scale(log_t_cdf_q)\n );\n real log_cdf_q = log_t_cdf_q == negative_infinity()\n ? negative_infinity() : log_t_cdf_q - log(d - pwindow);\n real log_density = primarycensored_log_diff_exp(\n log_t_cdf_d - log(d), log_cdf_q\n ) - log(pwindow);\n real log_capped = fmin(log_cdf, 0);\n real log_tail = fmin(log_capped, log1m_exp(log_capped));\n real log_limit = log(1e-8) + log_tail;\n if (log_capped > -log2()) {\n log_limit = fmax(log_limit, log(1e-15));\n }\n log_limit = fmin(\n log_limit, fmax(log(1e-6) + log_density, log(1e-13) + log_tail)\n );\n return log_error > log_limit;\n}\nreal primarycensored_loglogistic_lcdf(real y, real scale, real shape) {\n return -log1p_exp(-shape * (log(y) - log(scale)));\n}\nreal primarycensored_loglogistic_lccdf(real y, real scale, real shape) {\n return -log1p_exp(shape * (log(y) - log(scale)));\n}\nvector loglogistic_numeric_ode(real x, vector y, array[] real theta,\n array[] real x_r, array[] int x_i) {\n real u = x_r[1] - x;\n real log_delay;\n if (u <= 0) {\n log_delay = x_i[1] ? 0 : negative_infinity();\n } else if (x_i[1]) {\n log_delay = primarycensored_loglogistic_lccdf(u | theta[1], theta[2]);\n } else {\n log_delay = primarycensored_loglogistic_lcdf(u | theta[1], theta[2]);\n }\n return rep_vector(exp(log_delay - theta[3]), 1);\n}\nreal loglogistic_numeric_lcdf(data real d, array[] real params,\n data real pwindow) {\n if (d <= 0) return negative_infinity();\n real x_end = fmin(d, pwindow);\n real mid = d - 0.5 * pwindow;\n int use_survival = mid > 0\n && primarycensored_loglogistic_lcdf(mid | params[1], params[2])\n > -log2();\n real log_scale = use_survival\n ? 0\n : primarycensored_loglogistic_lcdf(d | params[1], params[2]);\n real integral = ode_rk45_tol(\n loglogistic_numeric_ode, rep_vector(0.0, 1), 0.0, {x_end}, 1e-13, 1e-20,\n 1000000, {params[1], params[2], log_scale}, {d}, {use_survival}\n )[1, 1];\n if (use_survival) {\n real survival = integral / pwindow;\n if (d < pwindow) {\n survival += 1 - d / pwindow;\n }\n return survival >= 1 ? negative_infinity() : log1m(survival);\n }\n return integral > 0 ? fmin(log(integral) + log_scale - log(pwindow), 0)\n : negative_infinity();\n}\nreal loglogistic_ratio_small(real a, real log_A, real tol) {\n real A = exp(log_A);\n real B = A / (1 + A);\n real weight = exp(-a * log1p(A));\n real total = weight / (a + 1);\n for (m in 0:5000) {\n weight *= B * (a + m) / (m + 1);\n real term = weight / (a + m + 2);\n total += term;\n if (term <= tol * total) break;\n }\n return total;\n}\nreal loglogistic_partial_at_one(real a) {\n if (a <= 0) return 0.5;\n return 0.5 * (a * (digamma((a + 1) / 2) - digamma(a / 2)) - 1);\n}\nreal loglogistic_tau_integral(real a, real tau, real tol) {\n real shift = 2 * a * tau;\n real tau_sq = square(tau);\n real previous = 0;\n real current = 1;\n real total = 1;\n for (j in 0:5000) {\n real following = (shift * current + (j - 1) * tau_sq * previous)\n / (j + 1);\n previous = current;\n current = following;\n total += current / (j + 2);\n if (current <= tol * total) break;\n }\n return 0.5 * tau * total;\n}\nreal loglogistic_tail_term(real a, int m, real log_A, real shrink) {\n real ell = log_A - log(3);\n real x = a - m - 1;\n real z = x * ell;\n if (abs(z) < 1) {\n // Series for (1 - exp(-z)) / z, to keep its derivative in a\n real exprel = abs(z) < 1e-2\n ? 1 - z / 2 + square(z) / 6 - z * square(z) / 24\n + square(square(z)) / 120 - z * square(square(z)) / 720\n : -expm1(-z) / z;\n return exp(-(m + 1) * log_A) * ell * exprel;\n }\n return (exp(-(m + 1) * log_A)\n - shrink * pow(3.0, -(m + 1))) / x;\n}\nreal loglogistic_ratio_large(real a, real log_A, real tol) {\n // M_a(3) / 3^a, the integral up to the start of the tail series\n real at_start = (loglogistic_partial_at_one(a)\n + loglogistic_tau_integral(a, 0.5, tol))\n * exp(-a * log(3));\n real shrink = exp(-a * (log_A - log(3)));\n real total = shrink * at_start;\n for (m in 0:200) {\n real term = (m + 1) * loglogistic_tail_term(a, m, log_A, shrink);\n total += m % 2 == 0 ? term : -term;\n if (abs(term) <= tol * abs(total)) break;\n }\n return (1 + exp(-log_A)) * total;\n}\nreal loglogistic_moment_ratio(real a, real log_A, real tol) {\n if (log_A <= 0) {\n return loglogistic_ratio_small(a, log_A, tol);\n }\n if (log_A > log(3)) {\n return loglogistic_ratio_large(a, log_A, tol);\n }\n // 1 < A <= 3\n real A = exp(log_A);\n real tau = (A - 1) / (A + 1);\n return (loglogistic_partial_at_one(a) + loglogistic_tau_integral(a, tau, tol))\n * (1 + A) / A * exp(-a * log_A);\n}\nvector primarycensored_loglogistic_uniform_terms(real t,\n array[] real params) {\n if (t <= 0) {\n return rep_vector(negative_infinity(), 2);\n }\n real shape = params[2];\n real log_ratio = shape * (log(t) - log(params[1]));\n real log_t_F_T = log(t) - log1p_exp(-log_ratio);\n return [\n log_t_F_T,\n log_t_F_T + log(loglogistic_moment_ratio(1 / shape, log_ratio, 1e-16))\n ]';\n}\nreal primarycensored_loglogistic_uniform_lcdf(data real d, real q,\n array[] real params,\n data real pwindow) {\n vector[2] terms_d = primarycensored_loglogistic_uniform_terms(d, params);\n vector[2] terms_q = primarycensored_loglogistic_uniform_terms(q, params);\n real log_cdf = primarycensored_uniform_lcdf_from_terms(\n terms_d, terms_q, pwindow\n );\n if (loglogistic_window_ill_conditioned(terms_d[1], terms_q[1], d, pwindow,\n log_cdf)) {\n return loglogistic_numeric_lcdf(d | params, pwindow);\n }\n return fmin(log_cdf, 0);\n}\nvector primarycensored_gauss_legendre_nodes() {\n return [\n 0.048307665687738324, 0.14447196158279649, 0.23928736225213706,\n 0.33186860228212767, 0.42135127613063533, 0.50689990893222936,\n 0.5877157572407623, 0.66304426693021523, 0.73218211874028971,\n 0.79448379596794239, 0.84936761373256997, 0.89632115576605209,\n 0.93490607593773967, 0.96476225558750639, 0.98561151154526838,\n 0.99726386184948157\n ]';\n}\nvector primarycensored_gauss_legendre_log_weights() {\n return [\n -2.3377969318791947, -2.3471775191742132, -2.3661171972104831,\n -2.3949868407368724, -2.4343797887852916, -2.4851635920714812,\n -2.5485636934665021, -2.626297960160425, -2.7207977655059796,\n -2.8355865687925124, -2.9759677006607701, -3.1503787891021764,\n -3.3733722294882642, -3.6733185431605069, -4.1181622816923937,\n -4.9591900848988502\n ]';\n}\nreal primarycensored_lognormal_tilt_one_panel(real a, real b, real mu,\n real sigma, real rho) {\n if (b <= a) return negative_infinity();\n real half = 0.5 * (b - a);\n real mid = 0.5 * (a + b);\n vector[16] nodes = primarycensored_gauss_legendre_nodes();\n vector[16] log_weights = primarycensored_gauss_legendre_log_weights();\n vector[32] terms;\n for (i in 1:16) {\n real z_up = mid + half * nodes[i];\n real z_down = mid - half * nodes[i];\n terms[i] = log_weights[i] - rho * exp(mu + sigma * z_up)\n - 0.5 * square(z_up);\n terms[16 + i] = log_weights[i] - rho * exp(mu + sigma * z_down)\n - 0.5 * square(z_down);\n }\n return log_sum_exp(terms) + log(half) - 0.5 * log(2 * pi());\n}\nreal primarycensored_lognormal_tilt_panel(real a, real b, real mu,\n real sigma, real rho) {\n if (b <= a) return negative_infinity();\n int n_panels = 1;\n while (n_panels * 1.8 < sigma) {\n n_panels += 1;\n }\n if (n_panels == 1) {\n return primarycensored_lognormal_tilt_one_panel(a, b, mu, sigma, rho);\n }\n real width = (b - a) / n_panels;\n vector[n_panels] panels;\n for (j in 1:n_panels) {\n panels[j] = primarycensored_lognormal_tilt_one_panel(\n a + (j - 1) * width, a + j * width, mu, sigma, rho\n );\n }\n return log_sum_exp(panels);\n}\nvector primarycensored_lognormal_tilt_bump(real mu, real sigma, real rho) {\n real w0 = lambert_w0(exp(log(rho) + 2 * log(sigma) + mu));\n real z0 = -w0 / sigma;\n real u0 = exp(mu - w0);\n real hi = fmin(\n z0 + sqrt(80 / (1 + w0)), fmax((log(u0 + 40 / rho) - mu) / sigma, -z0)\n );\n real lo = fmax(\n z0 - sqrt(80), -sqrt(square(z0) + 2 * w0 / square(sigma) + 80)\n );\n real log_total = log_sum_exp(\n primarycensored_lognormal_tilt_panel(lo, z0, mu, sigma, rho),\n primarycensored_lognormal_tilt_panel(z0, hi, mu, sigma, rho)\n );\n return [z0, lo, hi, log_total]';\n}\nvector log_tilt_transform_context(int dist_id, real xi, array[] real params) {\n if (dist_id == 1 && xi < 0) {\n return primarycensored_lognormal_tilt_bump(params[1], params[2], -xi);\n }\n return rep_vector(0, 4);\n}\nreal primarycensored_log_std_normal_cdf(real z) {\n if (z >= 0) return log(Phi(z));\n if (z > -37) return log(0.5 * erfc(-z / sqrt(2)));\n real inv_z2 = 1 / square(z);\n real series = 1 - inv_z2 * (1 - 3 * inv_z2 * (1 - 5 * inv_z2 * (\n 1 - 7 * inv_z2 * (1 - 9 * inv_z2 * (1 - 11 * inv_z2 * (\n 1 - 13 * inv_z2)))\n )));\n return std_normal_lpdf(z) - log(-z) + log(series);\n}\nvector primarycensored_moments_from_partial(real log_t, vector log_m) {\n real log_g1 = primarycensored_log_diff_exp(log_t + log_m[1], log_m[2]);\n real log_h = primarycensored_log_diff_exp(log_t + log_m[2], log_m[3]);\n real log_g2 = primarycensored_log_diff_exp(log_t + log_g1, log_h);\n real log_a = primarycensored_log_diff_exp(log_t + log_m[3], log_m[4]);\n real log_b = primarycensored_log_diff_exp(log_t + log_h, log_a);\n real log_g3 = primarycensored_log_diff_exp(log_t + log_g2, log_b);\n return [log_g1, log_g2, log_g3]';\n}\nvector primarycensored_gamma_tilt_moments(real t, real shape, real rate) {\n real log_t = log(t);\n real log_x = log_t + log(rate);\n // log_mk is the log of int_0^t u^k f(u) du\n vector[2] log_p01 = gamma_lcdf_logx_pair(log_x, shape);\n vector[2] log_p23 = gamma_lcdf_logx_pair(log_x, shape + 2);\n real log_m0 = log_p01[1];\n real log_m1 = log(shape) - log(rate) + log_p01[2];\n real log_m2 = log(shape) + log(shape + 1) - 2 * log(rate) + log_p23[1];\n real log_m3 = log(shape) + log(shape + 1) + log(shape + 2) - 3 * log(rate)\n + log_p23[2];\n return primarycensored_moments_from_partial(\n log_t, [log_m0, log_m1, log_m2, log_m3]'\n );\n}\nvector primarycensored_stacy_tilt_moments(real t, real shape, real scale,\n real k) {\n real log_x = shape * (log(t) - log(scale));\n real k1 = k + 1 / shape;\n real k2 = k + 2 / shape;\n real k3 = k + 3 / shape;\n real log_k = lgamma(k);\n real log_t = log(t);\n real log_m0 = gamma_lcdf_logx(log_x, k);\n real log_m1 = log(scale) + lgamma(k1) - log_k\n + gamma_lcdf_logx(log_x, k1);\n real log_m2 = 2 * log(scale) + lgamma(k2) - log_k\n + gamma_lcdf_logx(log_x, k2);\n real log_m3 = 3 * log(scale) + lgamma(k3) - log_k\n + gamma_lcdf_logx(log_x, k3);\n return primarycensored_moments_from_partial(\n log_t, [log_m0, log_m1, log_m2, log_m3]'\n );\n}\nvector primarycensored_tilt_moments(real t, int dist_id,\n array[] real params) {\n if (dist_id == 2) {\n if (t <= 0) return rep_vector(negative_infinity(), 3);\n return primarycensored_gamma_tilt_moments(t, params[1], params[2]);\n } else if (dist_id == 4) {\n if (t <= 0) return rep_vector(negative_infinity(), 3);\n return primarycensored_gamma_tilt_moments(t, 1, params[1]);\n } else if (dist_id == 18) {\n real mu = params[1];\n real sigma = params[2];\n real z = (t - mu) / sigma;\n real log_phi = std_normal_lpdf(z);\n real log_Phi = primarycensored_log_std_normal_cdf(z);\n real log_g1;\n real log_g2;\n real log_g3;\n if (z > -1) {\n // The direct form has a derivative at z = 0, unlike the log form\n log_g1 = log(exp(log_phi) + z * exp(log_Phi));\n log_g2 = log((square(z) + 1) * exp(log_Phi) + z * exp(log_phi));\n log_g3 = log(\n (square(z) + 2) * exp(log_phi) + z * (square(z) + 3) * exp(log_Phi)\n );\n } else {\n // The terms have opposite signs, so the log form keeps the tail\n log_g1 = primarycensored_log_diff_exp(log_phi, log(-z) + log_Phi);\n log_g2 = primarycensored_log_diff_exp(\n log1p(square(z)) + log_Phi, log(-z) + log_phi\n );\n log_g3 = primarycensored_log_diff_exp(\n log(square(z) + 2) + log_phi, log(-z) + log(square(z) + 3) + log_Phi\n );\n }\n return [\n log(sigma) + log_g1, 2 * log(sigma) + log_g2, 3 * log(sigma) + log_g3\n ]';\n } else if (dist_id == 3 || dist_id == 5) {\n if (t <= 0) return rep_vector(negative_infinity(), 3);\n return primarycensored_stacy_tilt_moments(\n t, params[1], params[2], dist_id == 3 ? 1.0 : params[3]\n );\n } else if (dist_id == 1) {\n if (t <= 0) return rep_vector(negative_infinity(), 3);\n real mu = params[1];\n real sigma = params[2];\n real log_t = log(t);\n real z = (log_t - mu) / sigma;\n vector[4] log_m;\n for (k in 0:3) {\n log_m[k + 1] = k * mu + 0.5 * square(k * sigma)\n + primarycensored_log_std_normal_cdf(z - k * sigma);\n }\n return primarycensored_moments_from_partial(log_t, log_m);\n }\n reject(\"Invalid distribution identifier: \", dist_id);\n}\nvector primarycensored_stacy_series(real t, real xi, real shape, real scale,\n real k) {\n real log_x = shape * (log(t) - log(scale));\n int alternating = xi < 0;\n real log_z = log(abs(xi)) + log(scale);\n real log_w = log(abs(xi)) + log(t);\n real w = exp(log_w);\n vector[2] unreliable = rep_vector(not_a_number(), 2);\n if (w >= 400) return unreliable;\n real log_k = lgamma(k);\n real log_pos = 0;\n real log_neg = negative_infinity();\n real log_first = 0;\n for (n in 0:400) {\n real s = k + n / shape;\n real log_term = n * log_z - lgamma(n + 1) + lgamma(s) - log_k\n + gamma_lcdf_logx(log_x, s);\n if (n == 0) {\n log_first = log_term;\n log_pos = log_term;\n } else if (alternating && n % 2 == 1) {\n log_neg = log_sum_exp(log_neg, log_term);\n } else {\n log_pos = log_sum_exp(log_pos, log_term);\n }\n real log_abs = log_sum_exp(log_pos, log_neg);\n // The transform is at most the CDF for negative xi\n if (alternating && log_abs - log_first > log(1e5)) return unreliable;\n if (n + 1 > w\n && log_term + log_w - log(n + 1 - w) < log_abs + log(1e-17)) {\n if (log_pos <= log_neg) return unreliable;\n real log_sum = primarycensored_log_diff_exp(log_pos, log_neg);\n if (log_abs - log_sum > log(1e5)) return unreliable;\n return [log_sum, log_abs - log_sum]';\n }\n }\n return unreliable;\n}\nvector primarycensored_lognormal_tilt_quadrature(real t, real mu,\n real sigma, real rho,\n vector bump) {\n real z = (log(t) - mu) / sigma;\n // The slope of the log integrand at t, positive left of the mode\n real slope = -rho * sigma * t - z;\n if (z < bump[1]) {\n real width = 80 / (slope + sqrt(square(slope) + 80));\n real log_lower = primarycensored_lognormal_tilt_panel(\n z - width, z, mu, sigma, rho\n );\n return [\n log_lower, primarycensored_log_diff_exp(bump[4], log_lower)\n ]';\n }\n real curvature = 1 + rho * square(sigma) * t;\n real width = 80 / (abs(slope) + sqrt(square(slope) + 80 * curvature));\n // The integrand is below exp(-40) of its value at t beyond the point where\n // the tilt alone has decayed by 40\n width = fmin(\n width, fmax(abs(z), (log(t + 40 / rho) - mu) / sigma) - z\n );\n real log_upper = primarycensored_lognormal_tilt_panel(\n z, z + width, mu, sigma, rho\n );\n return [primarycensored_log_diff_exp(bump[4], log_upper), log_upper]';\n}\nreal primarycensored_lognormal_tilt_series(real t, real mu, real sigma,\n real xi) {\n real z = (log(t) - mu) / sigma;\n real log_xi_mu = log(xi) + mu;\n real log_factorial = 0;\n real total = primarycensored_log_std_normal_cdf(z);\n for (k in 1:20000) {\n log_factorial += log(k);\n real log_term = k * log_xi_mu + 0.5 * square(k * sigma) - log_factorial\n + primarycensored_log_std_normal_cdf(z - k * sigma);\n total = log_sum_exp(total, log_term);\n if (k > xi * t && log_term < total - 41) return total;\n }\n reject(\n \"The lognormal tilt transform needs more than 20000 terms for t = \", t,\n \" and xi = \", xi, \". Use the numerical path.\"\n );\n}\nvector primarycensored_lognormal_tilt_pair(real t, real mu, real sigma,\n real xi, vector bump) {\n if (xi == 0) {\n if (t <= 0) return [negative_infinity(), 0]';\n real z = (log(t) - mu) / sigma;\n return [\n primarycensored_log_std_normal_cdf(z),\n primarycensored_log_std_normal_cdf(-z)\n ]';\n }\n if (xi > 0) {\n if (t <= 0) return [negative_infinity(), positive_infinity()]';\n return [\n primarycensored_lognormal_tilt_series(t, mu, sigma, xi),\n positive_infinity()\n ]';\n }\n if (t <= 0) return [negative_infinity(), bump[4]]';\n return primarycensored_lognormal_tilt_quadrature(t, mu, sigma, -xi, bump);\n}\nvector log_tilt_transform_pair_shared(real t, int dist_id, real xi,\n array[] real params, vector context) {\n if (dist_id == 2) {\n real shape = params[1];\n real rate = params[2];\n real log_total = -shape * log1m(xi / rate);\n if (t <= 0) return [negative_infinity(), log_total, 0]';\n return append_row(\n log_total + gamma_log_pq_logx(log(t) + log(rate - xi), shape), 0\n );\n } else if (dist_id == 4) {\n real rate = params[1];\n real tilted_rate = rate - xi;\n real log_total = -log1m(xi / rate);\n if (t <= 0) return [negative_infinity(), log_total, 0]';\n return [\n log_total + log1m_exp(-tilted_rate * t), log_total - tilted_rate * t, 0\n ]';\n } else if (dist_id == 18) {\n // The upper tail is the lower tail of the reflected normal\n real mu = params[1];\n real sigma = params[2];\n real z = (t - mu - xi * square(sigma)) / sigma;\n real log_total = xi * mu + 0.5 * square(xi * sigma);\n return [\n log_total + primarycensored_log_std_normal_cdf(z),\n log_total + primarycensored_log_std_normal_cdf(-z), 0\n ]';\n } else if (dist_id == 3 || dist_id == 5) {\n real shape = params[1];\n real scale = params[2];\n real k = dist_id == 3 ? 1.0 : params[3];\n if (xi == 0) {\n if (t <= 0) return [negative_infinity(), 0, 0]';\n return append_row(\n gamma_log_pq_logx(shape * (log(t) - log(scale)), k), 0\n );\n }\n if (t <= 0) return [negative_infinity(), not_a_number(), 0]';\n vector[2] series = primarycensored_stacy_series(t, xi, shape, scale, k);\n return [series[1], not_a_number(), series[2]]';\n } else if (dist_id == 1) {\n vector[2] pair = primarycensored_lognormal_tilt_pair(\n t, params[1], params[2], xi, context\n );\n return append_row(pair, 0);\n }\n reject(\"Invalid distribution identifier: \", dist_id);\n}\nvector log_tilt_transform_pair(real t, int dist_id, real xi,\n array[] real params) {\n return log_tilt_transform_pair_shared(\n t, dist_id, xi, params, log_tilt_transform_context(dist_id, xi, params)\n );\n}\nvector primarycensored_exptilt_terms(real t, int dist_id, real rho,\n array[] real params, vector context) {\n return append_row(\n log_tilt_transform_pair(t, dist_id, 0, params)[1:2],\n log_tilt_transform_pair_shared(t, dist_id, -rho, params, context)\n );\n}\nint primarycensored_exptilt_form(data real d, int dist_id, real rho,\n data real pwindow) {\n int positive = dist_has_positive_support(dist_id);\n if (abs(rho) * pwindow < (positive ? 1e-2 : 1e-3)\n && abs(rho) * (abs(d) + pwindow) < 0.1) {\n return 1;\n }\n if (positive && d < pwindow && abs(rho) * d < 5e-3) return 2;\n return 0;\n}\nreal primarycensored_tail_diff(real lower_d, real lower_q, real upper_d,\n real upper_q) {\n if (is_nan(upper_d) || is_nan(upper_q)) {\n return primarycensored_log_diff_exp(lower_d, lower_q);\n }\n // NaN from terms that both underflow is a zero difference\n if (is_nan(lower_q - lower_d) || is_inf(upper_q)\n || lower_q - lower_d <= upper_d - upper_q) {\n return primarycensored_log_diff_exp(lower_d, lower_q);\n }\n return primarycensored_log_diff_exp(upper_q, upper_d);\n}\nreal primarycensored_exptilt_lcdf_from_terms(vector terms_d, vector terms_q,\n data real d, real rho,\n data real pwindow) {\n if (is_nan(terms_d[3]) || is_nan(terms_q[3])) return not_a_number();\n real log_diff_f = primarycensored_tail_diff(\n terms_d[1], terms_q[1], terms_d[2], terms_q[2]\n );\n real log_diff_j = primarycensored_tail_diff(\n terms_d[3], terms_q[3], terms_d[4], terms_q[4]\n );\n real log_num;\n real log_den;\n if (rho > 0) {\n log_num = primarycensored_log_diff_exp(rho * d + log_diff_j, log_diff_f);\n log_den = log_diff_exp(rho * pwindow, 0);\n } else {\n log_num = primarycensored_log_diff_exp(log_diff_f, rho * d + log_diff_j);\n log_den = log1m_exp(rho * pwindow);\n }\n // log_sum_exp of two `-inf` has a NaN derivative\n if (terms_q[1] == negative_infinity() && log_num == negative_infinity()) {\n return negative_infinity();\n }\n // Rounding can put the log CDF above 0\n real result = fmin(log_sum_exp(terms_q[1], log_num - log_den), 0);\n int no_tail = is_nan(terms_d[4]) || is_nan(terms_q[4]);\n real log_size = rho * d + terms_d[3] + fmax(terms_d[5], terms_q[5]);\n if (no_tail && result > negative_infinity()\n && log_size - log_den > log(1e5)) {\n return not_a_number();\n }\n return result;\n}\nreal primarycensored_exptilt_small_window_lcdf_from_terms(\n vector moments_d, vector moments_q, real rho, data real pwindow\n) {\n real scale = moments_d[1];\n if (scale == negative_infinity()) return negative_infinity();\n real g1_q = exp(moments_q[1] - scale);\n real g2_d = exp(moments_d[2] - scale);\n real g2_q = exp(moments_q[2] - scale);\n real g3_d = exp(moments_d[3] - scale);\n real g3_q = exp(moments_q[3] - scale);\n real relative = (1 - g1_q)\n + 0.5 * rho * (g2_d - pwindow - g2_q - pwindow * g1_q)\n + square(rho) * ((g3_d - g3_q) / 6\n - pwindow * (g2_d + g2_q) / 4\n + square(pwindow) * (1 - g1_q) / 12);\n if (relative <= 0) return negative_infinity();\n return fmin(scale + log(relative) - log(pwindow), 0);\n}\nreal primarycensored_exptilt_small_delay_lcdf_from_terms(\n vector moments_d, real rho, data real pwindow\n) {\n if (moments_d[1] == negative_infinity()) return negative_infinity();\n real log_den = rho > 0 ? log_diff_exp(rho * pwindow, 0)\n : log1m_exp(rho * pwindow);\n return fmin(\n moments_d[1] + log(abs(rho)) - log_den\n + log1p(0.5 * rho * exp(moments_d[2] - moments_d[1])\n + square(rho) * exp(moments_d[3] - moments_d[1]) / 6),\n 0\n );\n}\nreal primarycensored_exptilt_lcdf(data real d, int dist_id,\n array[] real params, data real pwindow,\n real rho) {\n int positive = dist_has_positive_support(dist_id);\n if (positive && d <= 0) {\n return negative_infinity();\n }\n if (pwindow == 0) return dist_lcdf(d | params, dist_id);\n real q = d - pwindow;\n int form = primarycensored_exptilt_form(d, dist_id, rho, pwindow);\n if (form == 1) {\n return primarycensored_exptilt_small_window_lcdf_from_terms(\n primarycensored_tilt_moments(d, dist_id, params),\n primarycensored_tilt_moments(q, dist_id, params), rho, pwindow\n );\n }\n if (form == 2) {\n return primarycensored_exptilt_small_delay_lcdf_from_terms(\n primarycensored_tilt_moments(d, dist_id, params), rho, pwindow\n );\n }\n vector[4] context = log_tilt_transform_context(dist_id, -rho, params);\n real log_cdf = primarycensored_exptilt_lcdf_from_terms(\n primarycensored_exptilt_terms(d, dist_id, rho, params, context),\n primarycensored_exptilt_terms(q, dist_id, rho, params, context),\n d, rho, pwindow\n );\n if (is_nan(log_cdf)) {\n return log(\n primarycensored_numeric_cdf(d | dist_id, params, pwindow, 2, {rho})\n );\n }\n return log_cdf;\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 (check_for_exptilt(dist_id, primary_id)) {\n if (!check_for_tilt_transform(dist_id, -primary_params[1], params)) {\n reject(\n \"The tilted delay distribution does not exist for tilt \",\n primary_params[1], \". Use the numerical path, see \",\n \"check_for_analytical_params().\"\n );\n }\n return primarycensored_exptilt_lcdf(\n d | dist_id, params, pwindow, primary_params[1]\n );\n }\n if (check_for_tlogis(dist_id, primary_id)) {\n return primarycensored_tlogis_lcdf(\n d | dist_id, params, pwindow, primary_params[1], primary_params[2]\n );\n }\n if (check_for_gumbel(dist_id, primary_id)) {\n if (!check_for_gumbel_params(\n dist_id, params, primary_params[1], primary_params[2])) {\n reject(\n \"The truncated Gumbel solution does not apply for mu \",\n primary_params[1], \" and beta \", primary_params[2],\n \". Use the numerical path, see check_for_analytical_params().\"\n );\n }\n return primarycensored_gumbel_lcdf(\n d | dist_id, params, pwindow, primary_params[1], primary_params[2]\n );\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 == 29 && primary_id == 1) {\n return primarycensored_loglogistic_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, also for a finite negative L for delays\n // on the reals\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_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 check_for_gumbel(int dist_id, int primary_id) {\n return primary_id == 4 && dist_id == 18;\n}\nint gumbel_n_terms(real log_s0) {\n int n = 2;\n while (n < 200\n && n * log_s0 - lgamma(n + 1) - fmax(log_s0, 0) >= log(1e-18)) {\n n += 1;\n }\n return n;\n}\nint check_for_gumbel_params(int dist_id, array[] real params, real mu,\n real beta) {\n if (!check_for_gumbel(dist_id, 4)) return 0;\n if (!(beta > 0) || is_inf(beta) || is_nan(mu) || is_inf(mu)) return 0;\n real log_s0 = mu / beta;\n if (log_s0 > log(15)) return 0;\n real xi = gumbel_n_terms(log_s0) / beta;\n real size = xi * abs(params[1]) + 0.5 * square(xi * params[2]);\n return log(machine_precision()) + log1p(size) + exp(log_s0)\n <= log(1e-7);\n}\nreal gumbel_log_sum_exp(real a, real b) {\n if (a == negative_infinity()) return b;\n if (b == negative_infinity()) return a;\n return log_sum_exp(a, b);\n}\nreal gumbel_delay_lcdf(real x, array[] real params, int dist_id) {\n if (dist_has_positive_support(dist_id) && x <= 0) {\n return negative_infinity();\n }\n if (dist_id == 1) {\n if (lognormal_lcdf_underflows(x, params[1], params[2])) {\n return negative_infinity();\n }\n return primarycensored_log_std_normal_cdf(\n (log(x) - params[1]) / params[2]\n );\n } else if (dist_id == 4) {\n return log1m_exp(-params[1] * x);\n } else if (dist_id == 18) {\n return primarycensored_log_std_normal_cdf((x - params[1]) / params[2]);\n } else if (dist_id == 3) {\n if (params[1] * (log(x) - log(params[2])) < -700) {\n return negative_infinity();\n }\n } else if (dist_id == 2 || dist_id == 13) {\n real log_f = dist_id == 2\n ? gamma_lcdf_logx(log(x) + gamma_log_rate(params[2]), params[1])\n : gamma_lcdf_logx(log(x) - log2(), params[1] / 2);\n return log_f < -700 ? negative_infinity() : log_f;\n } else if (dist_id == 16 || dist_id == 19 || dist_id == 22) {\n // Inverse gamma, with the shape and scale of each family\n real shape = dist_id == 16 ? params[1] : params[1] / 2;\n real scale = dist_id == 16 ? params[2]\n : (dist_id == 19 ? 0.5 : params[1] * square(params[2]) / 2);\n real log_f = gamma_log_pq_logx(log(scale) - log(x), shape)[2];\n return log_f < -700 ? negative_infinity() : log_f;\n } else if (dist_id == 9) {\n if (x >= 1) {\n return 0;\n }\n if (x < params[1] / (params[1] + params[2])\n && params[1] * log(x) - log(params[1]) - lbeta(params[1], params[2])\n < -700) {\n return negative_infinity();\n }\n }\n return dist_lcdf(x | params, dist_id);\n}\nreal gumbel_numeric_delay(real u, data real d, data real pwindow, real mu,\n real beta) {\n return (d - pwindow) + beta * log1p(u * exp(-(mu - pwindow) / beta));\n}\nvector gumbel_numeric_point(real tau, real lo, real hi, real a,\n data real d, data real pwindow, real mu,\n real beta) {\n real x = lo + (hi - lo) * tau;\n if (mu >= pwindow) {\n return [gumbel_numeric_delay(a - log1m(x), d, pwindow, mu, beta), 0]';\n }\n return [d - x, tgumbel_lpdf(x | 0, pwindow, mu, beta)]';\n}\nvector primarycensored_gumbel_ode(real tau, vector y, data real d,\n data real pwindow, int dist_id,\n array[] real params, real mu, real beta,\n real lo, real hi, real a,\n real log_shift) {\n vector[2] point = gumbel_numeric_point(\n tau, lo, hi, a, d, pwindow, mu, beta\n );\n return rep_vector(\n (hi - lo)\n * exp(gumbel_delay_lcdf(point[1] | params, dist_id) + point[2] - log_shift),\n 1\n );\n}\nreal gumbel_numeric_piece_lcdf(data real d, int dist_id, array[] real params,\n data real pwindow, real mu, real beta, real a,\n real c, real log_shift) {\n if (log_shift == negative_infinity()) {\n return negative_infinity();\n }\n vector[1] y0 = rep_vector(0.0, 1);\n real integral = ode_rk45_tol(\n primarycensored_gumbel_ode, y0, 0.0, {1.0}, 1e-10, 1e-14, 200000,\n d, pwindow, dist_id, params, mu, beta, 0.0, c, a, log_shift\n )[1, 1];\n if (!(integral > 0)) {\n return negative_infinity();\n }\n return -a + log_shift + log(integral);\n}\nreal gumbel_numeric_spike_lcdf(data real d, int dist_id, array[] real params,\n data real pwindow, real mu, real beta) {\n real log_delta = (mu - pwindow) / beta + log_diff_exp(pwindow / beta, 0);\n real a = 0;\n if (dist_has_positive_support(dist_id) && d < pwindow) {\n // The delay CDF is 0 for u below the kink, where d - z is 0\n real log_uk = (mu - pwindow) / beta\n + log_diff_exp((pwindow - d) / beta, 0);\n // The CDF is below exp(-1000) beyond u of 1000\n if (log_uk >= log_delta || log_uk > log(1000)) {\n return negative_infinity();\n }\n a = exp(log_uk);\n }\n real top = exp(log_delta);\n // A piece of at most 36 ends within 2^-53 of 1 in the upper quantile.\n // The lengths are literals so that they are not differentiated.\n real log_total = negative_infinity();\n int finished = 0;\n for (k in 1:300) {\n int big = top - a > 48;\n int j = 0;\n int last = top - a <= 12;\n real len = last ? top - a : (big ? 36.0 : 12.0);\n real log_f_end = gumbel_delay_lcdf(\n gumbel_numeric_delay(a + len, d, pwindow, mu, beta) | params, dist_id\n );\n real log_f_mid = gumbel_delay_lcdf(\n gumbel_numeric_delay(a + len / 2, d, pwindow, mu, beta)\n | params, dist_id\n );\n for (i in 1:30) {\n if (!(log_f_end - log_f_mid > 3)) {\n break;\n }\n if (last) {\n // Continue with pieces that are at most half of what is left\n last = 0;\n while (12.0 / pow(2, j) > len / 2) {\n j += 1;\n }\n } else {\n j += 1;\n }\n len = (big ? 36.0 : 12.0) / pow(2, j);\n log_f_end = gumbel_delay_lcdf(\n gumbel_numeric_delay(a + len, d, pwindow, mu, beta) | params, dist_id\n );\n log_f_mid = gumbel_delay_lcdf(\n gumbel_numeric_delay(a + len / 2, d, pwindow, mu, beta)\n | params, dist_id\n );\n }\n // The lengths of the pieces before the last are passed as literals\n if (last) {\n log_total = gumbel_log_sum_exp(\n log_total,\n gumbel_numeric_piece_lcdf(\n d | dist_id, params, pwindow, mu, beta, a, -expm1(-(top - a)),\n log_f_end\n )\n );\n } else {\n log_total = gumbel_log_sum_exp(\n log_total,\n gumbel_numeric_piece_lcdf(\n d | dist_id, params, pwindow, mu, beta, a,\n -expm1(-((big ? 36.0 : 12.0) / pow(2, j))), log_f_end\n )\n );\n }\n a += len;\n // The CDF is below exp(-1000) beyond u of 1000\n if (log_total == negative_infinity() && a > 1000) {\n return negative_infinity();\n }\n real slope = 2 * (log_f_end - log_f_mid) / len;\n if (last || a >= top\n || (slope < 0.99\n && -a + log_f_end - log1m(slope) < log_total + log(1e-15))) {\n finished = 1;\n break;\n }\n }\n if (!finished) {\n reject(\n \"The truncated Gumbel numerical integral did not finish for d \", d,\n \", pwindow \", pwindow, \", mu \", mu, \" and beta \", beta,\n \". The log CDF is extremely small, below about -1000, so the delay \",\n \"distribution and the window are far apart.\"\n );\n }\n return log_total - tgumbel_log1m_exp_neg_exp(log_delta);\n}\nvector gumbel_numeric_z_limits(data real d, int dist_id, data real pwindow,\n real mu, real beta) {\n real log_sw = (mu - pwindow) / beta;\n real z_low = fmax(0, mu - beta * log(745 + exp(log_sw)));\n real z_high = fmin(pwindow, fmax(mu, 0) + 46 * beta);\n if (dist_has_positive_support(dist_id)) {\n z_high = fmin(z_high, d);\n }\n return [z_low, z_high]';\n}\nreal gumbel_z_log_integrand(real z, data real d, int dist_id,\n array[] real params, data real pwindow, real mu,\n real beta) {\n return gumbel_delay_lcdf(d - z | params, dist_id)\n + tgumbel_lpdf(z | 0, pwindow, mu, beta);\n}\nreal gumbel_numeric_z_half(data real d, int dist_id, array[] real params,\n data real pwindow, real mu, real beta,\n real z_from, real z_to, real log_shift) {\n vector[1] y0 = rep_vector(0.0, 1);\n return ode_rk45_tol(\n primarycensored_gumbel_ode, y0, 0.0, {1.0}, 1e-10, 1e-14, 200000,\n d, pwindow, dist_id, params, mu, beta, z_from, z_to, 0.0, log_shift\n )[1, 1];\n}\nreal gumbel_numeric_z_lcdf(data real d, int dist_id, array[] real params,\n data real pwindow, real mu, real beta) {\n vector[2] z_limits = gumbel_numeric_z_limits(d, dist_id, pwindow, mu, beta);\n if (!(z_limits[2] > z_limits[1])) {\n return negative_infinity();\n }\n real z_low = z_limits[1];\n real z_high = z_limits[2];\n real step = (z_high - z_low) / 16;\n real z_peak = z_low;\n real log_shift = gumbel_z_log_integrand(\n z_low, d, dist_id, params, pwindow, mu, beta\n );\n for (i in 1:16) {\n real z = z_low + step * i;\n real g = gumbel_z_log_integrand(z, d, dist_id, params, pwindow, mu, beta);\n if (g > log_shift) {\n log_shift = g;\n z_peak = z;\n }\n }\n if (log_shift == negative_infinity()) {\n return negative_infinity();\n }\n // Golden section search around the best grid point\n real inv_phi = 0.6180339887498949;\n real lo = fmax(z_low, z_peak - step);\n real hi = fmin(z_high, z_peak + step);\n real z1 = hi - inv_phi * (hi - lo);\n real z2 = lo + inv_phi * (hi - lo);\n real g1 = gumbel_z_log_integrand(z1, d, dist_id, params, pwindow, mu, beta);\n real g2 = gumbel_z_log_integrand(z2, d, dist_id, params, pwindow, mu, beta);\n for (i in 1:40) {\n if (g1 > g2) {\n hi = z2;\n z2 = z1;\n g2 = g1;\n z1 = hi - inv_phi * (hi - lo);\n g1 = gumbel_z_log_integrand(z1, d, dist_id, params, pwindow, mu, beta);\n } else {\n lo = z1;\n z1 = z2;\n g1 = g2;\n z2 = lo + inv_phi * (hi - lo);\n g2 = gumbel_z_log_integrand(z2, d, dist_id, params, pwindow, mu, beta);\n }\n }\n if (g1 > log_shift && g1 >= g2) {\n log_shift = g1;\n z_peak = z1;\n } else if (g2 > log_shift) {\n log_shift = g2;\n z_peak = z2;\n }\n real integral = 0;\n if (z_peak > z_low) {\n integral += -gumbel_numeric_z_half(\n d, dist_id, params, pwindow, mu, beta, z_peak, z_low, log_shift\n );\n }\n if (z_high > z_peak) {\n // An upper limit at the window end or d is passed as data, to keep it\n // out of the sensitivities, which are singular at d for a delay CDF\n // with an infinite slope at 0. Local variables would be parameters.\n if (z_high >= (dist_has_positive_support(dist_id) ? fmin(pwindow, d)\n : pwindow)) {\n integral += gumbel_numeric_z_half(\n d, dist_id, params, pwindow, mu, beta, z_peak,\n dist_has_positive_support(dist_id) ? fmin(pwindow, d) : pwindow,\n log_shift\n );\n } else {\n integral += gumbel_numeric_z_half(\n d, dist_id, params, pwindow, mu, beta, z_peak, z_high, log_shift\n );\n }\n }\n if (!(integral > 0)) {\n return negative_infinity();\n }\n return log_shift + log(integral);\n}\nreal primarycensored_gumbel_numeric_lcdf(data real d, int dist_id,\n array[] real params,\n data real pwindow, real mu,\n real beta) {\n if (dist_has_positive_support(dist_id) && d <= 0) {\n return negative_infinity();\n }\n real log_cdf = mu >= pwindow\n ? gumbel_numeric_spike_lcdf(d | dist_id, params, pwindow, mu, beta)\n : gumbel_numeric_z_lcdf(d | dist_id, params, pwindow, mu, beta);\n if (is_nan(log_cdf)) {\n return negative_infinity();\n }\n return fmin(log_cdf, 0);\n}\nreal primarycensored_gumbel_numeric_cdf(data real d, int dist_id,\n array[] real params,\n data real pwindow, real mu,\n real beta) {\n return exp(primarycensored_gumbel_numeric_lcdf(\n d | dist_id, params, pwindow, mu, beta\n ));\n}\nint gumbel_series_accepted(vector fit) {\n return fit[2] <= 1e-8;\n}\nvector primarycensored_gumbel_terms(real t, int dist_id, real beta,\n int n_terms, array[] real params) {\n vector[2 * (n_terms + 1)] terms;\n real inv_beta = inv(beta);\n for (n in 0:n_terms) {\n terms[(2 * n + 1):(2 * n + 2)] = log_tilt_transform_pair(\n t, dist_id, n * inv_beta, params\n )[1:2];\n }\n return terms;\n}\nreal gumbel_finite_abs(real x) {\n return is_inf(x) || is_nan(x) ? 0 : abs(x);\n}\nvector primarycensored_gumbel_lcdf_from_terms(vector terms_d, vector terms_q,\n data real d,\n data real pwindow, real mu,\n real beta, int n_terms) {\n real log_s0 = mu / beta;\n real log_c = -(d - mu) / beta;\n real log_f_q = terms_q[1];\n real log_delta0 = primarycensored_tail_diff(\n terms_d[1], terms_q[1], terms_d[2], terms_q[2]\n );\n // No delay mass in the window, so every term vanishes together\n if (log_delta0 == negative_infinity()) {\n return [log_f_q, 0]';\n }\n real log_pos = log1m_exp(-exp(log_s0)) + log_delta0;\n real log_neg = negative_infinity();\n real scale = 0;\n for (n in 1:n_terms) {\n int i = 2 * n + 1;\n real log_delta = primarycensored_tail_diff(\n terms_d[i], terms_q[i], terms_d[i + 1], terms_q[i + 1]\n );\n real log_a = n * log_c + log_delta - lgamma(n + 1);\n if (n % 2 == 1) {\n log_neg = gumbel_log_sum_exp(log_neg, log_a);\n } else {\n log_pos = gumbel_log_sum_exp(log_pos, log_a);\n }\n scale = fmax(\n scale,\n n * abs(log_c) + fmax(\n fmax(gumbel_finite_abs(terms_d[i]), gumbel_finite_abs(terms_d[i + 1])),\n fmax(gumbel_finite_abs(terms_q[i]), gumbel_finite_abs(terms_q[i + 1]))\n )\n );\n }\n real log_bracket = primarycensored_log_diff_exp(log_pos, log_neg);\n // log D = -s_w + log(1 - exp(-delta)), delta = s_w (exp(w / beta) - 1)\n real log_s_w = -(pwindow - mu) / beta;\n real log_d = -exp(log_s_w)\n + log1m_exp(-exp(log_s_w + log_diff_exp(pwindow / beta, 0)));\n real error;\n if (log_bracket == negative_infinity()) {\n error = positive_infinity();\n } else {\n int n_last = n_terms + 1;\n real log_rounding = log(machine_precision()) + log1p(scale)\n + gumbel_log_sum_exp(log_pos, log_neg) - log_bracket;\n real log_trunc = log_delta0 + n_last * log_s0 - lgamma(n_last + 1)\n - log1m(fmin(exp(log_s0) / (n_last + 1), 0.5))\n - log_bracket;\n error = exp(log_rounding) + exp(log_trunc);\n }\n real log_cdf = gumbel_log_sum_exp(log_f_q, log_bracket - log_d);\n return [fmin(log_cdf, 0), error]';\n}\nreal primarycensored_gumbel_lcdf(data real d, int dist_id,\n array[] real params, data real pwindow,\n real mu, real beta) {\n int n_terms = gumbel_n_terms(mu / beta);\n vector[2] fit = primarycensored_gumbel_lcdf_from_terms(\n primarycensored_gumbel_terms(d, dist_id, beta, n_terms, params),\n primarycensored_gumbel_terms(d - pwindow, dist_id, beta, n_terms, params),\n d, pwindow, mu, beta, n_terms\n );\n if (gumbel_series_accepted(fit)) {\n return fit[1];\n }\n return primarycensored_gumbel_numeric_lcdf(\n d | dist_id, params, pwindow, mu, beta\n );\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 if (dist_id == 29) return 1; // Log-logistic\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 } else if (primary_id == 3) {\n return tlogis_lcdf(p | 0, pwindow, primary_params[1], primary_params[2]);\n } else if (primary_id == 4) {\n return tgumbel_lcdf(p | 0, pwindow, primary_params[1], primary_params[2]);\n }\n reject(\"primary_lcdf: unsupported primary_id \", primary_id);\n}\nreal gamma_log_lead_logx(real log_x, real a) {\n real x = exp(log_x);\n if (a < 100) {\n return a * log_x - x - lgamma(a + 1);\n }\n real inv_a = inv(a);\n real inv_a2 = square(inv_a);\n real delta = inv_a * (1.0 / 12\n - inv_a2 * (1.0 / 360\n - inv_a2 * (1.0 / 1260\n - inv_a2 / 1680)));\n real r = x / a;\n // r - 1 is exact near 1, so log1p avoids cancellation\n real log_r = r > 0.5 && r < 2 ? log1p(r - 1) : log_x - log(a);\n return a * (log_r - (r - 1))\n - 0.5 * (log(a) + log(2 * pi())) - delta;\n}\nvector gamma_log_pq_logx(real log_x, real a) {\n if (!(a > 0) || is_inf(a)) {\n reject(\"gamma_log_pq_logx: shape must be positive finite, found a = \", a);\n }\n real x = exp(log_x);\n if (log_x == negative_infinity()) {\n return [negative_infinity(), 0]';\n }\n if (is_inf(x)) {\n return [0, negative_infinity()]';\n }\n real log_lead = gamma_log_lead_logx(log_x, a);\n // Both need about 9 sqrt(a) terms at worst, near x = a + 1\n real max_terms = 20 * sqrt(a) + 200;\n if (x < a + 1) {\n real term = 1;\n real total = 1;\n int n = 0;\n while (term >= 1e-17 * total) {\n n += 1;\n if (n > max_terms) {\n reject(\"gamma_log_pq_logx: series did not converge for log_x = \",\n log_x, \" and a = \", a);\n }\n term *= x / (a + n);\n total += term;\n }\n real log_p = log_lead + log(total);\n return [log_p, log1m_exp(log_p)]';\n }\n real b = x + 1 - a;\n real c = 1e300;\n real d = 1 / b;\n real h = d;\n real delta = 1;\n // err tracks |delta - 1|, which each step multiplies by\n // |an * d * d_prev / delta_prev|. Flooring |i - a| in an at 1 keeps it\n // from reaching 0 at an integer shape, where the fraction terminates but\n // its shape gradient does not.\n real err = 1;\n int i = 0;\n while (err >= 1e-15) {\n i += 1;\n if (i > max_terms) {\n reject(\"gamma_log_pq_logx: continued fraction did not converge for \",\n \"log_x = \", log_x, \" and a = \", a);\n }\n real an = -i * (i - a);\n real d_prev = d;\n b += 2;\n d = an * d + b;\n if (abs(d) < 1e-300) d = 1e-300;\n c = b + an / c;\n if (abs(c) < 1e-300) c = 1e-300;\n d = 1 / d;\n err *= i * fmax(abs(i - a), 1) * abs(d * d_prev / delta);\n delta = d * c;\n h *= delta;\n }\n real log_q = log_lead + log(a) + log(h);\n return [log1m_exp(log_q), log_q]';\n}\nreal gamma_lcdf_logx(real log_x, real a) {\n return gamma_log_pq_logx(log_x, a)[1];\n}\nvector gamma_lcdf_logx_pair(real log_x, real a) {\n if (is_inf(log_x)) {\n real log_p = gamma_lcdf_logx(log_x, a);\n return [log_p, log_p]';\n }\n real log_lead = gamma_log_lead_logx(log_x, a);\n if (exp(log_x) < a + 1) {\n real log_p_ap1 = gamma_lcdf_logx(log_x, a + 1);\n return [log_sum_exp(log_p_ap1, log_lead), log_p_ap1]';\n }\n real log_q = gamma_log_pq_logx(log_x, a)[2];\n return [log1m_exp(log_q), log1m_exp(log_sum_exp(log_q, log_lead))]';\n}\nreal gamma_log_rate(real rate) {\n if (!(rate > 0) || is_inf(rate)) {\n reject(\"Gamma rate must be positive finite, found \", rate);\n }\n return log(rate);\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 if (y < 0) {\n reject(\"gengamma_lcdf: y must be non-negative, found y = \", y);\n }\n if (y == 0) {\n return negative_infinity();\n }\n return gamma_lcdf_logx(shape * (log(y) - log(scale)), k);\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) {\n return gamma_lcdf_logx(log(delay) + gamma_log_rate(params[2]), params[1]);\n }\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 == 29) {\n return primarycensored_loglogistic_lcdf(delay | params[1], params[2]);\n }\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 if (primary_id == 3) {\n return tlogis_lpdf(x | xmin, xmax, params[1], params[2]);\n }\n if (primary_id == 4) {\n return tgumbel_lpdf(x | xmin, xmax, params[1], params[2]);\n }\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}\nint check_for_tlogis(int dist_id, int primary_id) {\n return primary_id == 3 && (dist_id == 2 || dist_id == 4 || dist_id == 18);\n}\nreal tlogis_log_bound(int n0, int M, real log_lo, real log_hi) {\n if (M == 0) return n0 * log_hi;\n if (n0 == 0) return M * (log1m_exp(log_lo) - log2());\n real log_x = fmin(fmax(log(1.0 * n0 / (n0 + M)), log_lo), log_hi);\n return M * (log1m_exp(log_x) - log2()) + n0 * log_x;\n}\narray[] int tlogis_series_terms(real log_lo, real log_hi, real log_tol) {\n if (is_inf(log_tol) || is_nan(log_tol)) return {-1, -1};\n for (k in 1:64) {\n int taper = to_int(floor(0.35 * k));\n if (tlogis_log_bound(k, 0, log_lo, log_hi) <= log_tol) return {k, 0};\n if (tlogis_log_bound(taper, k - taper, log_lo, log_hi) <= log_tol) {\n return {taper, k - taper};\n }\n if (tlogis_log_bound(0, k, log_lo, log_hi) <= log_tol) return {0, k};\n }\n return {-1, -1};\n}\narray[] int tlogis_plan(real location, real scale, data real pwindow,\n real log_mass) {\n real log_tol = log(1e-10) + log_mass;\n array[2] int pos = {0, 0};\n array[2] int neg = {0, 0};\n int pos_form = 0;\n if (location < 0) {\n pos_form = 1;\n pos = tlogis_series_terms(\n (location - pwindow) / scale, location / scale, log_tol\n );\n if (pos[1] < 0) return {0, 0, 0, 0, 0, 0};\n } else {\n if (location < pwindow) {\n pos_form = 2;\n pos = tlogis_series_terms(-(pwindow - location) / scale, 0, log_tol);\n if (pos[1] < 0) return {0, 0, 0, 0, 0, 0};\n }\n if (location > 0) {\n neg = tlogis_series_terms(\n -location / scale, fmin(0, (pwindow - location) / scale), log_tol\n );\n if (neg[1] < 0) return {0, 0, 0, 0, 0, 0};\n }\n }\n return {1, pos_form, pos[1], pos[2], neg[1], neg[2]};\n}\nint tlogis_is_valid(real location, real scale, data real pwindow) {\n if (!(pwindow > 0) || is_inf(pwindow)) return 0;\n return scale > 0 && !is_inf(scale) && !is_nan(location)\n && !is_inf(location);\n}\nint tlogis_n_pos(array[] int plan) {\n if (plan[2] == 0) return 0;\n return plan[3] + plan[4] - (plan[2] == 2 ? 1 : 0);\n}\nint tlogis_n_neg(array[] int plan) {\n return plan[5] + plan[6];\n}\nint tlogis_plan_applies(int dist_id, array[] real params, real scale,\n array[] int plan) {\n if (plan[1] == 0) return 0;\n if (tlogis_n_pos(plan) > 0\n && !check_for_tilt_transform(\n dist_id, tlogis_n_pos(plan) / scale, params\n )) {\n return 0;\n }\n if (tlogis_n_neg(plan) > 0\n && !check_for_tilt_transform(\n dist_id, -tlogis_n_neg(plan) / scale, params\n )) {\n return 0;\n }\n return 1;\n}\nint check_for_tlogis_params(int dist_id, array[] real params,\n array[] real primary_params,\n data real pwindow) {\n real location = primary_params[1];\n real scale = primary_params[2];\n if (!tlogis_is_valid(location, scale, pwindow)) return 0;\n real log_mass = tlogis_log_mass(0, pwindow, location, scale);\n return tlogis_plan_applies(\n dist_id, params, scale, tlogis_plan(location, scale, pwindow, log_mass)\n );\n}\nvector tlogis_weights(int n0, int M) {\n vector[n0 + M] weights = rep_vector(1, n0 + M);\n real tail = 0;\n // The upper tail from the top, so that small weights keep their precision\n for (i in 1:M) {\n int top = M - i + 1;\n tail += exp(lchoose(M, top) - M * log2());\n weights[n0 + top] = tail;\n }\n return weights;\n}\nmatrix primarycensored_tlogis_terms(real t, int dist_id, array[] real params,\n real scale, int n_pos, int n_neg,\n int use_pos, int use_neg) {\n matrix[2, 1 + n_pos + n_neg] terms = rep_matrix(0, 2, 1 + n_pos + n_neg);\n terms[, 1] = log_tilt_transform_pair(t, dist_id, 0, params)[1:2];\n if (use_pos) {\n for (n in 1:n_pos) {\n terms[, 1 + n] = log_tilt_transform_pair(\n t, dist_id, n / scale, params\n )[1:2];\n }\n }\n if (use_neg) {\n for (n in 1:n_neg) {\n terms[, 1 + n_pos + n] = log_tilt_transform_pair(\n t, dist_id, -n / scale, params\n )[1:2];\n }\n }\n return terms;\n}\nreal tlogis_terms_diff(matrix terms_hi, matrix terms_lo, int column) {\n return primarycensored_tail_diff(\n terms_hi[1, column], terms_lo[1, column],\n terms_hi[2, column], terms_lo[2, column]\n );\n}\nreal tlogis_weighted_sum(vector log_c, vector weights, real log_scale) {\n int n = rows(weights);\n real total = 0;\n for (i in 1:n) {\n total += (i % 2 == 1 ? 1 : -1) * weights[i] * exp(log_c[i] - log_scale);\n }\n return total;\n}\nvector tlogis_part_before_window(matrix terms_a, matrix terms_d,\n data real d, real location, real scale,\n vector weights) {\n int n = rows(weights);\n real log_df = tlogis_terms_diff(terms_d, terms_a, 1);\n vector[n] log_c;\n for (k in 1:n) {\n log_c[k] = k * location / scale + primarycensored_log_diff_exp(\n log_df, -k * d / scale + tlogis_terms_diff(terms_d, terms_a, 1 + k)\n );\n }\n real log_scale = max(log_c);\n if (is_inf(log_scale)) return [0, 0]';\n return [log_scale, tlogis_weighted_sum(log_c, weights, log_scale)]';\n}\nvector tlogis_part_above_location(matrix terms_a, matrix terms_u,\n data real d, real location, real scale,\n vector weights) {\n int n = rows(weights);\n real log_df = tlogis_terms_diff(terms_u, terms_a, 1);\n vector[n] log_c;\n for (k in 0:(n - 1)) {\n log_c[k + 1] = (k == 0 ? log_df\n : tlogis_terms_diff(terms_u, terms_a, 1 + k))\n - k * (d - location) / scale;\n }\n // The constant term of the series, L(0) times the mass of the delay\n real log_constant = log_inv_logit(-location / scale) + log_df;\n real log_scale = fmax(max(log_c), log_constant);\n if (is_inf(log_scale)) return [0, 0]';\n return [\n log_scale,\n tlogis_weighted_sum(log_c, weights, log_scale)\n - exp(log_constant - log_scale)\n ]';\n}\nvector tlogis_part_below_location(matrix terms_u, matrix terms_d, int n_pos,\n data real d, real location, real scale,\n vector weights) {\n int n = rows(weights);\n real log_df = tlogis_terms_diff(terms_d, terms_u, 1);\n vector[n] log_c;\n for (k in 1:n) {\n log_c[k] = -k * location / scale + primarycensored_log_diff_exp(\n k * d / scale + tlogis_terms_diff(terms_d, terms_u, 1 + n_pos + k),\n log_df\n );\n }\n real log_scale = max(log_c);\n if (is_inf(log_scale)) return [0, 0]';\n return [log_scale, tlogis_weighted_sum(log_c, weights, log_scale)]';\n}\nreal primarycensored_tlogis_lcdf_from_terms(\n matrix terms_a, matrix terms_u, matrix terms_d, data real d, real location,\n real scale, real log_mass, array[] int plan, vector w_pos, vector w_neg\n) {\n real log_phi = negative_infinity();\n if (plan[2] == 1) {\n vector[2] part = tlogis_part_before_window(\n terms_a, terms_d, d, location, scale, w_pos\n );\n if (part[2] > 0) log_phi = part[1] + log(part[2]);\n } else if (plan[2] == 2) {\n vector[2] part = tlogis_part_above_location(\n terms_a, terms_u, d, location, scale, w_pos\n );\n if (part[2] > 0) log_phi = part[1] + log(part[2]);\n }\n if (tlogis_n_neg(plan) > 0) {\n vector[2] part = tlogis_part_below_location(\n terms_u, terms_d, tlogis_n_pos(plan), d, location, scale, w_neg\n );\n if (part[2] > 0) log_phi = log_sum_exp(log_phi, part[1] + log(part[2]));\n }\n real log_f_a = terms_a[1, 1];\n // Deep enough into the lower tail every term underflows together. Both\n // are then `-inf` and log_sum_exp would differentiate to NaN.\n if (log_f_a == negative_infinity() && log_phi == negative_infinity()) {\n return negative_infinity();\n }\n // The CDF is at most 1. Rounding can put the log CDF a little above 0.\n return fmin(log_sum_exp(log_f_a, log_phi - log_mass), 0);\n}\nint tlogis_is_small_delay(int dist_id, data real d, data real pwindow,\n real scale) {\n return dist_has_positive_support(dist_id) && d > 0 && d < pwindow\n && d < 1e-4 * scale;\n}\nreal primarycensored_tlogis_small_delay_lcdf(data real d, int dist_id,\n array[] real params,\n real location, real scale,\n real log_mass) {\n vector[3] moments = primarycensored_tilt_moments(d, dist_id, params);\n if (moments[1] == negative_infinity()) return negative_infinity();\n real log_l1 = log_inv_logit(-location / scale)\n + log_inv_logit(location / scale) - log(scale);\n real ratio = (1 - 2 * inv_logit(-location / scale)) / scale;\n return fmin(\n moments[1] + log_l1 - log_mass\n + log1p(0.5 * ratio * exp(moments[2] - moments[1])),\n 0\n );\n}\nreal primarycensored_tlogis_lcdf(data real d, int dist_id,\n array[] real params, data real pwindow,\n real location, real scale) {\n if (dist_has_positive_support(dist_id) && d <= 0) {\n return negative_infinity();\n }\n if (!tlogis_is_valid(location, scale, pwindow)) {\n reject(\"The truncated logistic solution needs a finite location, \",\n \"and a positive finite scale and window.\");\n }\n real log_mass = tlogis_log_mass(0, pwindow, location, scale);\n array[6] int plan = tlogis_plan(location, scale, pwindow, log_mass);\n if (!tlogis_plan_applies(dist_id, params, scale, plan)) {\n reject(\n \"The truncated logistic solution does not apply for location \",\n location, \", scale \", scale, \" and window \", pwindow,\n \". Use the numerical path, see check_for_analytical_params().\"\n );\n }\n if (tlogis_is_small_delay(dist_id, d, pwindow, scale)) {\n return primarycensored_tlogis_small_delay_lcdf(\n d | dist_id, params, location, scale, log_mass\n );\n }\n int n_pos = tlogis_n_pos(plan);\n int n_neg = tlogis_n_neg(plan);\n int has_pos = plan[2] > 0;\n int has_neg = n_neg > 0;\n vector[plan[3] + plan[4]] w_pos = tlogis_weights(plan[3], plan[4]);\n vector[plan[5] + plan[6]] w_neg = tlogis_weights(plan[5], plan[6]);\n real a = d - pwindow;\n // The split point is d for a location at most 0 and a at least pwindow\n int u_is_d = location <= 0;\n int u_is_a = location >= pwindow;\n matrix[2, 1 + n_pos + n_neg] terms_a = primarycensored_tlogis_terms(\n a, dist_id, params, scale, n_pos, n_neg, has_pos, has_neg && u_is_a\n );\n matrix[2, 1 + n_pos + n_neg] terms_d = primarycensored_tlogis_terms(\n d, dist_id, params, scale, n_pos, n_neg, has_pos && u_is_d, has_neg\n );\n if (u_is_d) {\n return primarycensored_tlogis_lcdf_from_terms(\n terms_a, terms_d, terms_d, d, location, scale, log_mass, plan, w_pos,\n w_neg\n );\n } else if (u_is_a) {\n return primarycensored_tlogis_lcdf_from_terms(\n terms_a, terms_a, terms_d, d, location, scale, log_mass, plan, w_pos,\n w_neg\n );\n }\n return primarycensored_tlogis_lcdf_from_terms(\n terms_a,\n primarycensored_tlogis_terms(\n d - location, dist_id, params, scale, n_pos, n_neg, has_pos, has_neg\n ),\n terms_d, d, location, scale, log_mass, plan, w_pos, w_neg\n );\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_numeric_cdf(data real d, data int dist_id,\n array[] real params, data real pwindow,\n data int primary_id,\n array[] real primary_params) {\n if (primary_id == 4) {\n return primarycensored_gumbel_numeric_cdf(\n d | dist_id, params, pwindow, primary_params[1], primary_params[2]\n );\n }\n real lower_bound = dist_has_positive_support(dist_id)\n ? fmax(d - pwindow, 0) : 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(\n params, primary_params\n );\n array[4] int ids = {dist_id, primary_id, n_params, n_primary_params};\n\n vector[1] y = rep_vector(0.0, 1);\n if (primary_id != 3) {\n return ode_rk45(\n primarycensored_ode, y, lower_bound, {d}, theta, {d, pwindow}, ids\n )[1, 1];\n }\n // Solve in parts around a narrow truncated logistic density, with tighter\n // tolerances, so that the solver does not step over it\n vector[12] stops = append_row(\n tlogis_spike_times(d, pwindow, primary_params[1], primary_params[2]), d\n );\n real start = lower_bound;\n for (i in 1:12) {\n if (i < 12 && (stops[i] <= start || stops[i] >= d)) continue;\n y = ode_rk45_tol(\n primarycensored_ode, y, start, {stops[i]}, 1e-9, 1e-10, 1000000,\n theta, {d, pwindow}, ids\n )[1];\n start = stops[i];\n }\n return y[1];\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 and applies for these parameters\n if (check_for_analytical_delay(dist_id, params, primary_id,\n primary_params, d, pwindow)) {\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. For positive-support\n // delays the integral is 0 for d <= 0.\n if (dist_has_positive_support(dist_id) && d <= 0) {\n return 0;\n }\n result = primarycensored_numeric_cdf(\n d | dist_id, params, pwindow, primary_id, primary_params\n );\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_delay(dist_id, params, primary_id,\n primary_params, d, pwindow)) {\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 if (primary_id == 4) {\n // The truncated Gumbel integral is on the log scale, so a CDF below\n // 1e-300 is not lost\n result = primarycensored_gumbel_numeric_lcdf(\n d | dist_id, params, pwindow, primary_params[1], primary_params[2]\n );\n } else {\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\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}\nreal tgumbel_log_delta_window(real xmin, real xmax, real mu, real beta) {\n return -(xmax - mu) / beta + log_diff_exp((xmax - xmin) / beta, 0);\n}\nreal tgumbel_log1m_exp_neg_exp(real x) {\n if (x < -37) {\n return x;\n }\n if (x > 700) {\n return 0;\n }\n return log1m_exp(-exp(x));\n}\nreal tgumbel_lcdf(real x, real xmin, real xmax, real mu, real beta) {\n if (x <= xmin) {\n return negative_infinity();\n }\n if (x >= xmax) {\n return 0;\n }\n real log_delta_window = tgumbel_log_delta_window(xmin, xmax, mu, beta);\n real log_delta_lower = -(x - mu) / beta\n + log_diff_exp((x - xmin) / beta, 0);\n real log_delta_upper = -(xmax - mu) / beta\n + log_diff_exp((xmax - x) / beta, 0);\n real log_norm = tgumbel_log1m_exp_neg_exp(log_delta_window);\n return tgumbel_log1m_exp_neg_exp(log_delta_lower) - exp(log_delta_upper)\n - log_norm;\n}\nreal tgumbel_log_norm(real xmin, real xmax, real mu, real beta) {\n return tgumbel_log1m_exp_neg_exp(\n tgumbel_log_delta_window(xmin, xmax, mu, beta)\n );\n}\nreal tgumbel_lpdf(real x, real xmin, real xmax, real mu, real beta) {\n if (x < xmin || x > xmax) {\n return negative_infinity();\n }\n real log_s = -(x - mu) / beta;\n real exponent = 0;\n if (x < xmax) {\n exponent = exp(-(xmax - mu) / beta + log_diff_exp((xmax - x) / beta, 0));\n }\n return -log(beta) + log_s - exponent - tgumbel_log_norm(xmin, xmax, mu, beta);\n}\nreal tlogis_log_mass(real lo, real hi, real location, real scale) {\n if (lo >= location) {\n return primarycensored_log_diff_exp(\n log_inv_logit(-(lo - location) / scale),\n log_inv_logit(-(hi - location) / scale)\n );\n }\n return primarycensored_log_diff_exp(\n log_inv_logit((hi - location) / scale),\n log_inv_logit((lo - location) / scale)\n );\n}\nreal tlogis_lcdf(real x, real xmin, real xmax, real location, real scale) {\n if (x <= xmin) {\n return negative_infinity();\n }\n if (x >= xmax) {\n return 0;\n }\n return tlogis_log_mass(xmin, x, location, scale)\n - tlogis_log_mass(xmin, xmax, location, scale);\n}\nreal tlogis_lpdf(real x, real xmin, real xmax, real location, real scale) {\n if (x < xmin || x > xmax) {\n return negative_infinity();\n }\n return logistic_lpdf(x | location, scale)\n - tlogis_log_mass(xmin, xmax, location, scale);\n}\nvector tlogis_spike_times(data real d, data real pwindow, real location,\n real scale) {\n vector[11] multiples = [30, 12, 6, 3, 1.5, 0, -1.5, -3, -6, -12, -30]';\n return d - (fmin(fmax(location, 0), pwindow) + scale * multiples);\n}"This is useful when you want to extract a self-contained set of Stan functions for use in another project.
2.2.2 Exploring function dependencies
To understand which functions a particular Stan function depends on, use pcd_stan_function_deps():
pcd_stan_function_deps("primarycensored_lpmf")## [1] "dist_has_positive_support"
## [2] "primarycensored_log_normalizer"
## [3] "check_for_uniform_terms"
## [4] "check_for_analytical"
## [5] "check_for_tilt_transform"
## [6] "check_for_exptilt"
## [7] "check_for_gumbel"
## [8] "gumbel_n_terms"
## [9] "check_for_gumbel_params"
## [10] "check_for_tlogis"
## [11] "tlogis_log_bound"
## [12] "tlogis_series_terms"
## [13] "tlogis_plan"
## [14] "tlogis_is_valid"
## [15] "tlogis_n_pos"
## [16] "tlogis_n_neg"
## [17] "tlogis_plan_applies"
## [18] "primarycensored_log_diff_exp"
## [19] "tlogis_log_mass"
## [20] "check_for_tlogis_params"
## [21] "check_for_analytical_params"
## [22] "check_for_tilt_transform_at"
## [23] "check_for_analytical_delay"
## [24] "expgrowth_cdf"
## [25] "expgrowth_lcdf"
## [26] "tgumbel_log_delta_window"
## [27] "tgumbel_log1m_exp_neg_exp"
## [28] "tgumbel_lcdf"
## [29] "tlogis_lcdf"
## [30] "primary_lcdf"
## [31] "primary_lcdf_vec"
## [32] "discretestep_lcdf"
## [33] "hazards_to_pmf"
## [34] "discretehazard_lcdf"
## [35] "primarycensored_uniform_lcdf_from_terms"
## [36] "gamma_log_lead_logx"
## [37] "gamma_log_pq_logx"
## [38] "gamma_lcdf_logx"
## [39] "gamma_lcdf_logx_pair"
## [40] "gamma_log_rate"
## [41] "primarycensored_gamma_uniform_terms"
## [42] "primarycensored_gamma_uniform_lcdf"
## [43] "lognormal_lcdf_underflows"
## [44] "primarycensored_lognormal_uniform_terms"
## [45] "primarycensored_lognormal_uniform_lcdf"
## [46] "log_weibull_g"
## [47] "primarycensored_weibull_uniform_terms"
## [48] "primarycensored_weibull_uniform_lcdf"
## [49] "primarycensored_gengamma_uniform_terms"
## [50] "primarycensored_gengamma_uniform_lcdf"
## [51] "loglogistic_log_error_scale"
## [52] "loglogistic_window_ill_conditioned"
## [53] "primarycensored_loglogistic_lcdf"
## [54] "primarycensored_loglogistic_lccdf"
## [55] "loglogistic_numeric_ode"
## [56] "loglogistic_numeric_lcdf"
## [57] "loglogistic_ratio_small"
## [58] "loglogistic_partial_at_one"
## [59] "loglogistic_tau_integral"
## [60] "loglogistic_tail_term"
## [61] "loglogistic_ratio_large"
## [62] "loglogistic_moment_ratio"
## [63] "primarycensored_loglogistic_uniform_terms"
## [64] "primarycensored_loglogistic_uniform_lcdf"
## [65] "primarycensored_gauss_legendre_nodes"
## [66] "primarycensored_gauss_legendre_log_weights"
## [67] "primarycensored_lognormal_tilt_one_panel"
## [68] "primarycensored_lognormal_tilt_panel"
## [69] "primarycensored_lognormal_tilt_bump"
## [70] "log_tilt_transform_context"
## [71] "primarycensored_log_std_normal_cdf"
## [72] "primarycensored_moments_from_partial"
## [73] "primarycensored_gamma_tilt_moments"
## [74] "primarycensored_stacy_tilt_moments"
## [75] "primarycensored_tilt_moments"
## [76] "primarycensored_stacy_series"
## [77] "primarycensored_lognormal_tilt_quadrature"
## [78] "primarycensored_lognormal_tilt_series"
## [79] "primarycensored_lognormal_tilt_pair"
## [80] "log_tilt_transform_pair_shared"
## [81] "log_tilt_transform_pair"
## [82] "primarycensored_exptilt_terms"
## [83] "primarycensored_exptilt_form"
## [84] "primarycensored_tail_diff"
## [85] "primarycensored_exptilt_lcdf_from_terms"
## [86] "primarycensored_exptilt_small_window_lcdf_from_terms"
## [87] "primarycensored_exptilt_small_delay_lcdf_from_terms"
## [88] "pstep_lcdf"
## [89] "phazard_lcdf"
## [90] "gengamma_lcdf"
## [91] "dist_lcdf"
## [92] "gumbel_log_sum_exp"
## [93] "gumbel_delay_lcdf"
## [94] "gumbel_numeric_delay"
## [95] "tgumbel_log_norm"
## [96] "tgumbel_lpdf"
## [97] "gumbel_numeric_point"
## [98] "primarycensored_gumbel_ode"
## [99] "gumbel_numeric_piece_lcdf"
## [100] "gumbel_numeric_spike_lcdf"
## [101] "gumbel_numeric_z_limits"
## [102] "gumbel_z_log_integrand"
## [103] "gumbel_numeric_z_half"
## [104] "gumbel_numeric_z_lcdf"
## [105] "primarycensored_gumbel_numeric_lcdf"
## [106] "primarycensored_gumbel_numeric_cdf"
## [107] "expgrowth_lpdf"
## [108] "tlogis_lpdf"
## [109] "primary_lpdf"
## [110] "primarycensored_ode"
## [111] "tlogis_spike_times"
## [112] "primarycensored_numeric_cdf"
## [113] "primarycensored_exptilt_lcdf"
## [114] "gumbel_series_accepted"
## [115] "primarycensored_gumbel_terms"
## [116] "gumbel_finite_abs"
## [117] "primarycensored_gumbel_lcdf_from_terms"
## [118] "primarycensored_gumbel_lcdf"
## [119] "tlogis_weights"
## [120] "primarycensored_tlogis_terms"
## [121] "tlogis_terms_diff"
## [122] "tlogis_weighted_sum"
## [123] "tlogis_part_before_window"
## [124] "tlogis_part_above_location"
## [125] "tlogis_part_below_location"
## [126] "primarycensored_tlogis_lcdf_from_terms"
## [127] "tlogis_is_small_delay"
## [128] "primarycensored_tlogis_small_delay_lcdf"
## [129] "primarycensored_tlogis_lcdf"
## [130] "primarycensored_analytical_lcdf_raw"
## [131] "primarycensored_apply_truncation"
## [132] "primarycensored_truncation_bounds"
## [133] "primarycensored_analytical_lcdf"
## [134] "primarycensored_analytical_cdf"
## [135] "primarycensored_cdf"
## [136] "primarycensored_lcdf"
## [137] "primarycensored_lpmf"The result is ordered so that dependencies come before the functions that use them, with the requested function last. This can help you understand the structure of the package’s Stan code or decide which functions you need to include.
2.3 Linking the Stan functions to your workflow
2.3.1 Writing functions to a file
One option for using Stan functions is to write them to a file and then compile them using cmdstanr. This is a good approach as it means that once the functions are written they can be used in the same way as any other stan functions you might use. The downside is that it may mean more work keeping up to date with changes to the functions. We can do this using the pcd_load_stan_functions() function.
expgrowth_rng_file <- file.path(tempdir(), "expgrowth_rng.stan")
exp_model <- pcd_load_stan_functions(
"expgrowth_rng",
write_to_file = TRUE,
output_file = expgrowth_rng_file,
wrap_in_block = TRUE
)## Stan functions written to: /tmp/RtmpggKq0Z/expgrowth_rng.stanThis can now be compiled and used in the same way as any other cmdstanr model.
model <- cmdstan_model(expgrowth_rng_file)
model## functions {
## // Stan functions from primarycensored version 1.6.0
## real expgrowth_rng(real xmin, real xmax, real r) {
## real u = uniform_rng(0, 1);
## if (abs(r) < 1e-10) {
## return xmin + u * (xmax - xmin);
## }
## return log(u * exp(r * xmax) + (1 - u) * exp(r * xmin)) / r;
## }
## }Alternatively, you could use #include expgrowth_rng.stan in a stan file functions block to include the function along with the path to that file as with any other stan file (see here).
2.3.2 Including the functions directly via include_paths
Rather than writing the functions to a file it is also possible to include them directly in the stan file using the include_paths argument to cmdstan_model(). This is useful if you don’t to clutter your model with the stan code from primarycensored and want automatic updating of the functions. To demonstrate we will first write a small model has has expgrowth.stan in its include paths (rather than writing it to a file and then including it). The first step is find the file and path for the expgrowth_rng function.
pcd_stan_files("expgrowth_rng")## [1] "expgrowth.stan"With that done we now write stan wrapper model.
expgrowth_stan_file <- file.path(tempdir(), "expgrowth.stan")
writeLines(
text = c(
"functions {",
"#include expgrowth.stan",
"}",
"generated quantities {",
" real y = expgrowth_rng(0, 1, 0.4);",
"}"
),
con = expgrowth_stan_file
)We can now use this file to compile a model. Note that we need to include the path to the primarycensored Stan functions using the include_paths argument to cmdstan_model().
model <- cmdstan_model(expgrowth_stan_file, include_paths = pcd_stan_path())
model## functions {
## #include expgrowth.stan
## }
## generated quantities {
## real y = expgrowth_rng(0, 1, 0.4);
## }We can then sample from the model (we set fixed_param = TRUE here as our toy example doesn’t require MCMC sampling).
samples <- model$sample(chains = 1, fixed_param = TRUE)## Running MCMC with 1 chain...
##
## Chain 1 Iteration: 1 / 1000 [ 0%] (Sampling)
## Chain 1 Iteration: 100 / 1000 [ 10%] (Sampling)
## Chain 1 Iteration: 200 / 1000 [ 20%] (Sampling)
## Chain 1 Iteration: 300 / 1000 [ 30%] (Sampling)
## Chain 1 Iteration: 400 / 1000 [ 40%] (Sampling)
## Chain 1 Iteration: 500 / 1000 [ 50%] (Sampling)
## Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 1 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
samples## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## y 0.52 0.54 0.29 0.37 0.06 0.95 1.00 979 8022.4 Using Stan functions directly in R
Whilst it is possible to use Stan functions directly in R it is not recommended for most use cases (use the R functions in primarycensored instead). However, it can be useful to understand what is going on under the hood or for exploration (indeed we use this internally in primarycensored to check our functions against the R implementations). To do this we use the expose_functions() method on our already compiled model. This can take some time (~30 seconds) to compile all of the functions.
model$expose_functions(global = TRUE)We can now use the function in R. Note that this may get slightly more complicated if our stan function depends on other stan functions (i.e. you need to have those included in your compiled model as well).
expgrowth_rng(0, 1, 0.4)## [1] 0.1772395