I just wanted to make two comments, then first on interval censoring in survival modeling and the second on posterior retrodictive checks.
Interval censoring is quite natural to incorporate into survival models. Recall that a survival model is defined by some probability density function \rho(x; \theta) for an event taking place at x \in [0, \infty). A typical survival observation is defined by not only by an event occurring at x but also having not occurred at any previous x; mathematically the corresponding observational model is given by
\begin{align*}
\pi(x; \theta)
&= \; \text{Probability of not occurring before }x
\\
&\quad \cdot \text {Probability density of occurring at }x
\\
&= \; (1 - \mathbb{P}_{\rho}[I(0, x)])
\\
&\quad \cdot \rho(x; \theta)
\\
&=
\left[ \int_{x}^{\infty} \mathrm{d} x' \, \rho(x'; \theta) \right] \cdot \rho(x; \theta),
\end{align*}
where I(0, x) is the interval from 0 to x.
A right censored observation occurs when x is only known to be larger than some threshold, x_{R}. Mathematically the observational model for these observations is given by the probability
\begin{align*}
P_{R}(\theta)
&=
\text{Probability of occurring after }x_{R}
\\
&= \mathbb{P}_{\rho}[I(x_{R}, \infty)])
\\
&=\int_{x_{R}}^{\infty} \mathrm{d} x' \, \rho(x'; \theta).
\end{align*}
Superficially we just drop the \rho(x; \theta) term compared to the non-censored observation and then x = x_{R}.
Mathematically an interval censored observation is defined similarly. Here observations are known to be above x_{R} (right-censored) but also below x_{L} (left censored) so that the observational model is given by the probability
\begin{align*}
P_{I}(\theta)
&=
\text{Probability of occurring after }x_{R}\text{ and before }x_{L}
\\
&= \mathbb{P}_{\rho}[I(x_{R}, x_{L})])
\\
&= \int_{x_{R}}^{x_{L}} \mathrm{d} x' \, \rho(x'; \theta).
\end{align*}
By the way most of Stan’s log CDF implementations are just naive logarithms of the CDF implementations so
rho_cdf(x_{L} | theta) - rho_cdf(x_{R} | theta)
will be more stable and efficient than
log_diff_exp(rho_lcdf(x_{L} | theta), rho_lcdf(x_{R} | theta),
as is done in the current Stan program.
One way to simulate all of these observations is to draw directly from the occurrence model
\tilde{x} \sim \rho(x; \theta),
and then return a censored observation if \tilde{x} falls into any of the censoring intervals, as is done in Felipe’s code. Another possibility is to partition the observational space into disjoint possibilities, compute the probability of each, sample one of those possibilities, and then sample any additional structure within that possibility.
For example let’s say that we observe survival events fully up until x_{1} at which point observations are censored between x_{1} and x_{2}, x_{2} and $x_{3}, and the right censored above $x_{3}. We now have four basic possibilities:
\begin{align*}
1: p_{1} = P[ 0 < x < x_{1} ] &= \int_{0}^{x_{1}} \mathrm{d} x' \, \rho(x'; \theta)
\\
2: p_{2} = P[ x_{1} < x < x_{2} ] &= \int_{x_{1}}^{x_{2}} \mathrm{d} x' \, \rho(x'; \theta)
\\
3: p_{3} = P[ x_{2} < x < x_{3} ] &= \int_{x_{2}}^{x_{3}} \mathrm{d} x' \, \rho(x'; \theta)
\\
4: p_{4} = P[ x_{3} < x ] &= \int_{x_{3}}^{\infty} \mathrm{d} x' \, \rho(x'; \theta)
\end{align*}
Once we’ve evaluated (p_{1}, p_{2}, p_{3}, p_{4}) we can draw from a categorical RNG using those weights to determine which kind of observation we make. If we draw 2, 3, or 4 then we’re done and if we draw 1 then we further have to simulate
\tilde{x} \sim \rho(x; \theta)
subject to \tilde{x} < x_{1}.
While this might seem like a mess compared to the more straightforward simulation approach it provides motivation for summary statistics that make for very useful retrodictive checks. In particular given the number of observations of each type we can construct empirical probabilities
\begin{align*}
\tilde{p}_{1} &= \frac{ n_{1} }{ N } \approx p_{1}
\\
\tilde{p}_{2} &= \frac{ n_{2} }{ N } \approx p_{2}
\\
\tilde{p}_{3} &= \frac{ n_{3} }{ N } \approx p_{3}
\\
\tilde{p}_{4} &= \frac{ n_{4} }{ N } \approx p_{4}.
\end{align*}
both from the observed data and posterior predictive simulations. This then defines four scalar posterior retrodictive checks between the observed empirical probabilities and the distribution of simulated posterior predictive empirical probabilities. Given enough observations the uncensored interval [0, x_{1}) can also be discretized into smaller intervals; the corresponding empirical probability retrodictive checks better probe the shape of the survival function there.
One can also combine probabilities, for example looking at p_{1}, p_{1} + p_{2}, p_{1} + p_{2} + p_{3}, and p_{1} + p_{2} + p_{3} + p_{4} which recreates a empirical CDF over the censoring intervals. This begins to approach the quantile residual summary statistics mentioned above only using the natural structure of the problem to define intervals that doesn’t require auxiliary simulations to work around the censoring. Also by directly comparing the empirical cumulative probabilities one doesn’t rely on assumed calibration behavior that isn’t guaranteed to hold.
