• The primary event falls somewhere in a window \([0, w_P]\)
• It need not be uniform, but a uniform primary is fair when the window is short
• Mix the delay over where in the window the event sat, a convolution on the CDF
• A uniform primary has closed forms; we are working on more analytical solutions

\[ F_{\mathrm{PEC}}(t) = \int_{0}^{w_P} g(\tau)\, F(t - \tau)\, d\tau \]

• The secondary event is recorded to the day as well
• The probability of a daily bin is a difference of the primary-censored CDF
• It doubles as the probability mass function for daily discrete-time models

\[ \Pr\!\big(T \in [n, n + w_S]\big) = F_{\mathrm{PEC}}(n + w_S) - F_{\mathrm{PEC}}(n) \]

• A pair enters the data only if the secondary event happened before the cutoff \(C\)
• Condition on being observed, so divide by the CDF at the horizon
• Short delays are over-represented early in an epidemic

\[ L^{\mathrm{RT}} = \frac{F_{\mathrm{PEC}}(n + w_S) - F_{\mathrm{PEC}}(n)} {F_{\mathrm{PEC}}(C - P)} \]

• At a cutoff seven days after the primary events, and 5.3 days at fourteen

• d, p, q and r prefixes, so it reads like dgamma. The delay family is whichever p* function you hand it
• S3 dispatch picks the closed form when the pair has one, and stats::integrate when it does not
• The primary event distribution defaults to dunif

• A distribution is the string "norm", with d and p pasted on the front and looked up by name
• Three copies of this censored objective in mledist.R, branching on log support and on weights
• To reach it, primarycensored invents a distribution called pcens_dist and hides two closures in an empty environment for the name lookup to find

• 1,215 lines of Stan next to 2,703 lines of R
• The R lookup table names 25 delay families. dist_lcdf in Stan has 18 branches, because Stan has no CDF for the other seven
• Every fix lands twice. test-stan-dist_lcdf.R is 246 lines whose job is checking that the two halves still agree

Error in function tanh_sinh<double>::integrate: The
tanh_sinh quadrature evaluated your function at a
singular point and got -inf. Please narrow the bounds
of integration or check your function for
singularities.

Warning: 243 chain(s) finished unexpectedly!

• Post-warmup, in sampling. 243 of 256 chains gone, and the ones that lived returned the right answer
• integrate_1d errors rather than rejecting, so the chain dies instead of the proposal

• xc, Stan’s high-precision distance argument, turned singularities into an integral that would not solve at all
• Then d from parameter to data, integrating over the delay instead of the window, a guard on tiny delays, and the tolerance from 1e-6 to 1e-2. Twelve per cent at best
• The recast to an ODE worked because the Stan devs had recently reworked AD support—not because it is a good or bad approach

// before
integrate_1d(primary_censored_integrand, ...)
// after
ode_rk45(primary_censored_ode, y0, ...)[1, 1]

• A downstream model cannot depend on the package. It has to physically contain the text of the functions
• So the copying code, R/pcd-stan-tools.R, is 453 lines of R that reads Stan source, builds a call graph and sorts it topologically. It counts braces by hand
• epinowcast carries a 521 line copy with its own licence file beside it. epidist splices the whole 1,103 line bundle in at build time

• Delays vary by district, by age and by wave. brms, by Paul-Christian Bürkner, is one way (via Stan) to get partial pooling and time-varying effects in R
• Every brms family has a mu. Stan’s gamma_lpdf takes the shape first, so the intercept we read as the mean was the shape
• To learn what order brms used, epidist compiles a dummy model, runs a regex over the generated Stan and reads the order back out
• That string is then gsub-ed into our own Stan template, with dist_id and three more. Five passes before anything compiles

• Each adjustment is a wrapper on a Distributions.jl object. primary_censored is the \(w_P\) integral, upper extends Distributions.jl’s truncated at \(C\), interval_censored is \(w_S\)
• The stack is still a distribution, so Turing takes it as it is
• A new distribution would just work. For example, cdf(primary_censored(Frechet(2, 1), Uniform(0, 1)), 3.0) returns 0.847
• In Stan Frechet needs a branch, an integer and a rebuild

# the Gamma closed form
MethodError: no method matching
  _gamma_inc(::Dual, ::Float64, ::Int64)

# and once a rule was written for it
Enzyme forward -> [-0.21736, -0.21296, 0.19054]
ForwardDiff    -> [-0.23692, -0.21296, 0.19054]

• SpecialFunctions leaves the shape-parameter partial of gamma_inc not implemented, so our Gamma closed form would not differentiate
• Adaptive quadrature also errored under auto-differentiation
• Lifting a correct ChainRulesCore rule into Enzyme returned a wrong shape partial and two right ones, with no error

• 63 scenarios and a CI workflow per backend. DifferentiationInterfaceTest, by Guillaume Dalle and Adrian Hill, is what made any of it testable, and it is good
• The 408 are for one Gamma CDF derivative, across gamma_ad.jl and five autodiff extensions
• The default solver is GaussLegendre(; n = 64), because fixed nodes trace through every backend. I got there by hacking around the integral in sad confusion

• Nine repositories on GitHub depend on primarycensored, two of them on CRAN. Seven are from my own orgs. Two are CDC’s, both as suggests
• This means multiple public-health agencies use it as well as researchers
• For example, the NEJM correspondence on Bundibugyo virus disease (DRC 2026) used epidist, and therefore primarycensored — doi:10.1056/NEJMc2608070
• Everything that depends on CensoredDistributions.jl is mine, Sam Brand’s or Sebastian Funk’s

• Somebody should build brms for Julia, building on the composability the language offers. Chat to me if you are interested — a few people are already working on it — or join the #bayes channel on the Julia Slack
• CensoredDistributions.jl, the Julia package
• primarycensored, the R and Stan one
• epiaware.org, the rest of the packages
• A workflow for infectious disease modelling
• The prompts that wrote these decks