2.1. Definition and notation

Let $(\Omega, \mathcal{A}, {P})$ be a probability space. Consider the two-dimensional stochastic process $(S_i, X_i)$ , $i \in \mathbb{N}_0$ , on $(\Omega, \mathcal{A}, {P})$ with values in $\{ 0, 1 \} \times \mathbb{R}^+$ . Here, $S_i$ denotes the ith state that the process is in, and $0 = X_0 \leq X_1 \leq \cdots$ are the jump times. Call a time interval that the process spends in state 0 a Z-phase and in state 1 a Y-phase, in analogy to [Reference Van Lieshout and Markwitz25]. We set $S_0 = 1$ .

The non-homogeneous alternating renewal process can simply be seen as a special case of the semi-Markov process with only two states. In the semi-Markov paradigm, the tuple $(S_n, X_n)_{n=1}^{\infty}$ defines a non-homogeneous alternating renewal process if

(2.1) \begin{align} \mathbb{P}(S_{n+1} = j&, X_{n+1} \leq x\mid (S_0, X_0) = (s_0, x_0), \ldots, (S_n, X_n) = (s_n, x_n)) \nonumber \\ &= \mathbb{P}(S_{n+1} = j, X_{n+1} \leq x\mid (S_n, X_n) = (s_n, x_n)) \nonumber \\ &= \mathbb{P}(S_{n+1} = j, X_{n+1} - X_n \leq x - x_n \mid (S_n, X_n) = (s_n, x_n)) \nonumber \\ &= \mathbb{P}(S_{1} = j, X_{1} - X_0 \leq x - x_n \mid (S_0, X_0) = (s_n, x_n)), \end{align}

i.e. the joint conditional probability distribution of the sojourn time

(2.2) \begin{align} T_{n+1} = X_{n+1} - X_{n} \end{align}

and the next state $S_{n+1}$ depends only on the nth state $S_n$ and its jump time $X_n$ , not on the entire history of the process [Reference Çinlar6, Reference Janssen and Manca12, Reference Korolyuk and Swishchuk16, Reference Ross23] nor on the index n. This process is only Markov at the jump times, hence the name semi-Markov.

It follows from (2.1) that the distribution of a non-homogeneous alternating renewal process is completely specified by the starting state (or its probability distribution) and a semi-Markov kernel G that describes the transition rates from state i to state j. Formally, for $\tau \geq 0$ , $x\geq 0$ , and $i,j \in \{ 0,1\}$ ,

(2.3) \begin{equation} G_{ij}(x, \tau) = \mathbb{P}(S_{n+1} = j,\,T_{n+1} \leq \tau \mid S_{n} = i, X_{n} = x), \end{equation}

regardless of $n = 0, 1, \ldots$ Since only two transitions are possible, we proceed by using the notation $G_{Y}(x, \tau) =G_{10}(x, \tau)$ and $G_{Z}(x, \tau) = G_{01}(x, \tau)$ , where the subscript denotes the state that the process is in after jump time x. In the remainder of this paper, we shall assume that, for all $x\geq 0$ , $G_Y(x, \cdot)$ and $G_Z(x, \cdot)$ are absolutely continuous with respect to Lebesgue measure and write $g_Y(x,\cdot)$ and $g_Z(x, \cdot)$ respectively for their Radon–Nikodym derivatives.

2.2. Conditional intensities and non-explosion conditions

The conditional intensity, also known as the stochastic intensity, of a temporal point process describes the infinitesimal conditional probability of occurrence given the history of the process [Reference Karr14]. More precisely, for $n=0, 1, \ldots$ and $0 = x_0 \leq x_1 \leq \ldots \leq x_n\leq x$ ,

\begin{equation*} \lambda_{n+1}(x;\ x_1, \ldots, x_n) \, {\mathrm{d}} x = \mathbb{P}(X_{n+1} \leq x+{\mathrm{d}} x\mid X_{n+1} \geq x, X_0 = 0, X_1 = x_1, \ldots, X_n = x_n). \end{equation*}

For a non-homogeneous alternating renewal process, the $\lambda_{n+1}(\cdot\,; \cdot)$ are closely related to the hazard rates of the sojourn times. To see this, recall that $S_0 = 1$ and assume that $n+1$ is odd. Then, using (2.2), the conditional intensity of the jump process at time x given jumps at times $0 \leq x_1 \leq x_2 \leq \cdots \leq x_n$ can be simplified as

\begin{align*} \lambda_{n+1}(x;\ x_1, \ldots, x_n)\,{\mathrm{d}} x & = \mathbb{P}(X_{n+1} \leq x+ {\mathrm{d}} x\mid X_{n+1} \geq x, X_n = x_n, S_n = 1) \\ & = \frac{\mathbb{P}(x - x_n \leq T_{n+1} \leq x - x_n + {\mathrm{d}} x\mid X_n = x_n, S_n = 1)}{ \mathbb{P}(T_{n+1} \geq x-x_n\mid X_n = x_n, S_n = 1)} \\ & = \frac{g_Y(x_n, x - x_n) \, {\mathrm{d}} x }{1-G_Y(x_n, x - x_n)}\end{align*}

whenever well-defined, using the absolute continuity of the semi-Markov kernel. Note that this is exactly the hazard rate of the sojourn times. When $G_Y(x_n, x-x_n) = 1$ , the conditional intensity is set to zero. For even $n+1$ , a similar argument holds with $g_Z$ and $G_Z$ instead of $g_Y$ and $G_Y$ .

The conditional intensity is a convenient tool to guarantee the existence of the process. Indeed, [Reference Haezendonck and De Vylder10] developed suitable comparison criteria under which explosion, the situation in which there are infinitely many jumps in a finite time span, can be prevented. More formally, the following proposition holds.

the probability of explosion at or before time x of the point process defined by $\lambda$ is at most as big as that of the point process defined by $\lambda^*$ . Under the same conditions, for all $n\in{\mathbb{N}}$ and $x\geq 0$ , $\mathbb{P}(X_n \leq x ) \leq \mathbb{P}(X_{n}^* \leq x)$ .

Below, we establish existence for two common families of sojourn time distributions, the Gamma and the Weibull. We remark that in the semi-Markov paradigm, due to the Markov property, $\lambda_{n+1}(x;\ x_n) = \lambda_{n+1}(x;\ x_1, \ldots, x_n)$ , as there is no dependence on previous values $x_1, \ldots, x_{n-1}$ of $X_n$ . We will therefore use this notation in what follows.

The proof is deferred to Appendix A.1. Important special cases include $k_T(x) = 1$ for exponential distributions or, more generally, $k_T(x) \in \mathbb{N}$ corresponding to Erlang-distributed phases, where T can be either Y or Z.

For the proof, see Appendix A.2. The case where $k = 1$ corresponds to exponential sojourn times.

2.3. Renewal function: Existence and boundedness

The process counting the number of cycles having occurred by time $t \geq 0$ can be written as $N(t) = \sup\{n \in \mathbb{N}_0\colon X_{2n} \leq t\}$ , where a cycle is an interval of time within which each state occurs once. Thus, N is a counting measure. The distribution of $X_{2n}$ , the jump time after completing the nth cycle, is, for $n\in\mathbb{N}_0$ and $t\geq 0$ ,

\begin{equation*} F_{2n}(t) = \mathbb{P}\Bigg( \sum_{i=1}^{2n} T_i \leq t \Bigg) = \mathbb{P}(X_{2n} \leq t) = \mathbb{P}(N(t) \geq n).\end{equation*}

The renewal function is defined, analogously to that of the classic alternating renewal process, as $M(t) = \mathbb{E}N(t)$ , $t \geq 0$ [Reference Janssen and Manca12]. In our case,

\begin{equation*} M(t) = \mathbb{E}N(t) = \sum_{n=0}^{\infty} \mathbb{P}(N(t) \gt n) = \sum_{n=1}^{\infty} \mathbb{P}(X_{2n} \leq t) = \sum_{n=1}^{\infty} F_{2n}(t)\end{equation*}

is a 2n-fold convolution.

The following corollaries to Propositions 2.2 and 2.3 hold.

For the proof, see Appendix A.3.

For the proof, see Appendix A.4.

2.4. Age and excess distributions

Now that the theoretical groundwork for the censoring mechanism has been laid, we proceed by determining the joint distribution of age and excess. The age A(t) is the time elapsed since the last phase change, and B(t), the excess, is the time remaining until the next phase change. For all t where the process is in state 0, or the Z-phase, we assume that the occurrence time can be observed perfectly. Therefore we only consider age and excess with respect to state 1, or the Y-phase. This also provides us with a concrete censoring mechanism: depending on the phase within which a point falls, we note either the exact occurrence time, or the interval corresponding to the age and excess functions. Obtaining their joint distribution allows us to specify the likelihood of intervals based on their starting point and length in terms of the semi-Markov kernel $G_Y$ . See Figure 1 for a visualisation of the age and excess functions.

Figure 1. Visualisation of a non-homogeneous alternating renewal process with initial values $S_0 = 1$ and $X_0 = 0$ . At the dotted line, one cycle has passed – i.e. the process has taken both possible state values. The jump times correspond to a change of state. Since t falls in state 1, a non-zero age A(t) and excess B(t) are recorded. If t were to fall in state 0, an exact time would be recorded.

Proposition 2.4. (Age and excess.) Consider a non-homogeneous alternating renewal process $(S_n, X_n)_{n=1}^{\infty}$ with values in $\{0,1\} \times\mathbb{R}^+$ with $S_0=1$ , $X_0=0$ , semi-Markov kernels $G_Y$ and $G_Z$ , and associated counting measure N(t), $t\geq 0$ . Let the age process with respect to the Y-phase be $A(t) = (t - X_{2N(t)})\,\mathbf{1}\{X_{2N(t) + 1} \gt t\}$ and define the excess with respect to the Y-phase as $B(t) = (X_{2N(t) + 1} - t)\,\mathbf{1}\{X_{2N(t) + 1} \gt t\}$ , where $X_{2N(t)}$ is the jump time immediately after N(t) cycles have been completed. Then, for $t\geq 0$ , $0\leq x \leq t$ , and $z\geq 0$ ,

(2.4) \begin{align} & \mathbb{P}(A(t) \leq x, B(t) \leq z) \nonumber \\ & = G_Y(0,t) - \int_{t-x}^{t} [1 - G_Y(s, t + z - s)]\,{\mathrm{d}} M(s) - \int_0^{t-x} [1 - G_Y(s, t-s)]\,{\mathrm{d}} M(s) \nonumber \\ & \quad + \mathbf{1}\{x=t\}[G_Y(0, t+z) - G_Y(0,t)]. \end{align}

For the proof, see Appendix A.5.

From Proposition 2.4 we conclude that the probability that time $t \geq 0$ falls in a Z-phase is given by

(2.5) \begin{equation} w_t = \mathbb{P}(A(t) \leq 0, B(t) \leq 0) = G_Y(0,t) - \int_0^{t} [1 - G_Y(s, t-s)]\,{\mathrm{d}} M(s).\end{equation}

This case constitutes the atomic part of (2.4). The singular component on the line $x=t$ has total mass $1 - G_Y(0,t)$ and represents the case that t falls before the first jump of the alternating renewal process.

The absolutely continuous component of (2.4) can be written as

\begin{align*}\int_0^x \int_0^z g_Y(t-u, u+v) m(t-u)\, {\mathrm{d}} u\,{\mathrm{d}} v\end{align*}

provided that the Radon–Nikodym derivatives m of M and $g_Y$ of $G_Y$ exist. Recall that in our proposed censoring mechanism, when t falls in a Y-phase the entire interval $[t-A(t), t+B(t)]$ is reported, which may be parametrised by the left-most point $t-A(t)$ and length $A(t) + B(t)$ . Suppose that $A(t) = u$ , $B(t) = v$ , and apply the change of variables $a = t- u$ and $l = u + v$ . We find that the joint probability density function of left-most point and length is

(2.6) \begin{equation} q_t(a, l) = \frac{m(a)g_Y(a, l)}{\int_0^t [1 - G_Y(s, t-s)]\,{\mathrm{d}} M(s)} \mathbf{1} \{ 0 \leq a \leq t \leq a+l;\ l\geq 0 \}\end{equation}

upon scaling.

Proposition 2.5. (Marginal and conditional distribution.) Let $g_Y$ and m be as before, and let (A,L) be distributed according to $q_t(a,l)$ given by (2.6). Then the marginal probability density function of A at $a \in [0, t]$ is

(2.7) \begin{equation} f_t(a) = \frac{m(a)[1-G_Y(a,t-a)]}{\int_{0}^t [1 - G_Y(s, t-s)]\, {\mathrm{d}} M(s) }\end{equation}

and the conditional probability density function of L given $A=a$ is, for $l \in [t-a, \infty)$ ,

(2.8) \begin{equation} f_{t,\,L\,|A=a}(l) = \frac{g_Y(a, l)}{1-G_Y(a,t-a)}.\end{equation}

For the proof, see Appendix A.6.

3.1. Model formulation

The ensemble of potentially censored occurrence times can be mathematically formalised as a marked point process [Reference Daley and Vere-Jones7]. The ground process of points represents the uncensored event occurrences, which we model by a Markov point process [Reference Van Lieshout24] defined by a probability density with respect to a unit-rate Poisson process. Temporal variations can be taken into account as well as interactions between the points. Each point is subsequently marked, independently of other points, either by an atom at the point when it is observed perfectly, or by the interval in which the point lies in case of censoring. The mark kernel that governs the random censoring is based on the distribution of age and excess in a non-homogeneous alternating renewal process.

Formally, let $\mathcal{X}$ be an open set on the real line. The state space $\mathcal{N}_{\mathcal{X}}$ of a point process X consists of finite sets $\{x_1, x_2, \ldots, x_n\} \subset\mathcal{X}$ , $n \in \mathbb{N}_0$ , which we equip with the Borel $\sigma$ -algebra of the weak topology [Reference Daley and Vere-Jones7, Appendix A2]. Let p be a measurable, non-negative function on $\mathcal{X}$ that integrates to unity, and $\sim$ be a symmetric, reflexive relation on $\mathcal{X}$ . A point process X on $\mathcal{X}$ having probability density p with respect to a unit-rate Poisson process is Markov with respect to $\sim$ if [Reference Ripley and Kelly22, Reference Van Lieshout24]:

1. (i) p is hereditary, i.e. $p(\mathbf{x}) \gt 0$ implies that $p(\mathbf{y}) \gt 0$ for all subsets $\mathbf{y}$ of $\mathbf{x}$ , and

2. (ii) the conditional intensity, defined as ${p(\mathbf{x} \cup \{ t \} )}/{p(\mathbf{x})}$ , with $a/0 = 0$ for $a \geq 0$ , depends only on the neighbourhood $\{ x \in \mathbf{x} \colon x \sim t \}$ of t in $\mathbf{x}$ for every $t\in \mathcal{X}\setminus \mathbf{x}$ and every $\mathbf{x} = \{ x_1, \ldots, x_n\} \subset \mathcal{X}$ for which $p(\mathbf{x}) \gt 0$ .

An interaction function is a family $\phi_0, \phi_1, \phi_2,\ldots$ of non-negative functions $\phi_i$ defined on configurations of i points that take the value one whenever the configuration contains a pair $\{ x_1, x_2 \}$ of points that are unrelated, i.e. $x_1 \not\sim x_2$ . By the Hammersley–Clifford theorem [Reference Ripley and Kelly22], a Markov density p can be factorised as

(3.1) \begin{equation}p(\mathbf{x}) = \prod_{\mathbf{y}\subset\mathbf{x}} \phi_{|\mathbf{y}|}(\mathbf{y})\end{equation}

for some interaction function $\phi_i$ , writing $|\cdot|$ for cardinality. The function $\phi_1(x)$ can be used to model temporal variations in the likelihood of events occurring. Higher-order terms $\phi_2, \phi_3, \ldots$ govern interactions between pairs, triples, or tuples of points.

The points x in a realisation $\mathbf{x}$ of X are marked independently according to a mark kernel $\nu(\cdot \mid x)$ on ${\mathbb{R}}\times{\mathbb{R}}^+$ . A mark (a, l) represents an interval $[a,\ a+l]$ that starts at a and has length l. The mark kernel $\nu$ , which describes the conditional distribution of the mark given a location in time [Reference Daley and Vere-Jones7], formalises the non-homogeneous alternating renewal process censoring discussed in Section 2. For demonstrative purposes, we assumed a starting time of 0, which we now set to $-\infty$ . Doing so also allows us to ignore the singular component. Hence, the appropriate time-dependent mark kernel $\nu(\cdot \mid x)$ , $x\in\mathcal{X}$ , for a Borel subset $A \subset \mathbb{R} \times \mathbb{R}^+$ is

(3.2) \begin{align} \nu(A \mid x) & = \bigg(1 - {\int_{-\infty}^x[1 - G_Y(s, x-s)]\,{\mathrm{d}} M(s)}\bigg) \delta(\{(x,0)\} \cap A) \nonumber \\ & \quad + \int_{-\infty}^x\int_{x-a}^{\infty} \mathbf{1}\{ (a,l) \in A \} \, G_Y(a,{\mathrm{d}} l)\,{\mathrm{d}} M(a). \end{align}

Write W for the marked point process defined by $p{(}{\cdot}{)}$ and $\nu(\cdot \mid \cdot)$ [Reference Daley and Vere-Jones7, Proposition 4.IV]. A realisation $\mathbf{w}$ is of the form

\begin{align*}\mathbf{w} = \{w_1, w_2, \ldots, w_n\} = \{(x_1, (a_1, l_1)), (x_2, (a_2, l_2)), \ldots, (x_n, (a_n, l_n))\}\end{align*}

for $a_i \leq x_i \leq a_i+l_i$ for all $i = 1, 2, \ldots, n$ . We denote the set of realisations by $\mathcal{N}_{ {\mathcal{X}}\times (\mathbb{R} \times \mathbb{R}^+) }$ .

The model description is complete by noting that the observable pattern of marks after censoring is $U = \bigcup_{(x_i, (a_i,l_i)) \in W} \{ (a_i, l_i) \}$ . To obtain the probability distribution of U, write, for F in the Borel $\sigma$ -algebra of the weak topology on $\mathcal{N}_{ {\mathcal{\mathbb{R} \times \mathbb{R}^+}}}$ ,

\begin{align*} \mathbb{P}(U\in F \mid X = \mathbf{x} ) = \int_{(\mathbb{R}\times\mathbb{R}^+)^n} \mathbf{1}(\{ (a_1, l_1), \ldots, (a_n, l_n) \} \in F) \prod_{i=1}^n {\mathrm{d}} \nu((a_i, l_i)\mid x_i),\end{align*}

where $\mathbf{x} = \{x_1, \ldots, x_n\}$ , and then take the expectation with respect to X.

3.2. Conditional distribution

Write $\mathbf{u}$ for a realisation of the interval set U. We are interested in the conditional distributions of X and W given U.

Theorem 3.1. (Conditional distribution.) Let W be a marked point process with ground process X on the open set ${\mathcal{X}} \subset {\mathbb{R}}$ defined by its probability density function p with respect to the distribution of a unit-rate Poisson process having independent marks distributed according to the mark kernel $\nu( \cdot \mid x) $ for $x \in {\mathcal{X}}$ given by (3.2). Let ${\mathbf{u}}$ be a realisation of U that consists of an atomic part $\{ (a_1, 0), \ldots, (a_m, 0) \}$ , $m\in{\mathbb{N}}_0$ , and a non-atomic part $\{ (a_{m+1}, l_{m+1}), \ldots, (a_n, l_n) \}$ , $n \geq m$ . Then the conditional distribution of X given $U = {\mathbf{u}}$ satisfies, for A in the Borel $\sigma$ -algebra of the weak topology on ${\mathcal{N}}_{\mathcal{X}}$ ,

\begin{align*} \mathbb{P}(X \in A\mid U=\mathbf{u}) = & \int_{\mathcal{X}^{n-m}} p(\{a_1, \ldots, a_m, x_{1}, \ldots, x_{n-m}\}) \mathbf{1}_A(\{a_1, \ldots, a_m, x_{1}, \ldots, x_{n-m}\}) \\ & \quad \sum_{\substack{D_1, \ldots, D_{n-m} \\ \cup_j \{ D_j \} = \{1, \ldots, n-m\}}} \prod_{i=1}^{n-m} \mathbf{1}_{[a_{m+i}, a_{m+i} + l_{m+i}]}(x_{D_i})\,{\mathrm{d}} x_i \times c(\mathbf{u}),\end{align*}

provided that the normalisation constant

\begin{equation*} c(\mathbf{u}) = \Bigg[\int_{\mathcal{X}^{n-m}} p(\mathbf{x} \cup \{a_1, \ldots, a_m\}) \sum_{\substack{D_1, \ldots, D_{n-m} \\ \cup_j \{ D_j \} = \{1, \ldots, n-m\}}} \prod_{i=1}^{n-m} \mathbf{1}_{[a_{m+i}, a_{m+i} + l_{m+i}]}(x_{D_i})\,{\mathrm{d}} x_i\Bigg]^{-1}\end{equation*}

exists in $(0,\infty)$ .

Proof. We must prove (see, e.g., [Reference Daley and Vere-Jones7, Appendix 3.1]), for each A in the Borel $\sigma$ -algebra of ${\mathcal{N}}_{\mathcal{X}}$ with respect to the weak topology and each F in the Borel $\sigma$ -algebra of the weak topology on ${\mathcal{N}}_{{\mathbb{R}}\times{\mathbb{R}}^+}$ , that

\begin{equation*}{\mathbb{E}}[ {\mathbf{1}}_F(U) {\mathbb{P}}( X \in A \mid U )] ={\mathbb{E}}[ {\mathbf{1}}_F(U) {\mathbf{1}}_A(X)],\end{equation*}

as in [Reference Van Lieshout and Markwitz25, (4)]. From the model description, writing $w_t$ for the modification of (2.5) over $({-}{\infty,} t)$ , $\ell$ for Lebesgue measure, and $| \cdot |$ for cardinality,

(3.3) \begin{align} & \mathbb{E}[\mathbf{1}_F(U)\mathbf{1}_A(X)] \nonumber \\ & = \sum_{n=0}^{\infty}\frac{{\mathrm{e}}^{-\ell(\mathcal{X})}}{n!} \nonumber \\ & \qquad \int_{\mathcal{X}^n}\mathbf{1}_A(\mathbf{x})p(\{x_1,\ldots,x_n\})\sum_{C_0\subset\{1,\ldots,n\}} \frac{1}{(n - |C_0|)!} \prod_{i\in C_0} w_{x_i} \nonumber \\ & \qquad\quad \int_{(\mathbb{R}\times\mathbb{R}_0^+)^{n-|C_0|}}\mathbf{1}_F(\{(a_1,l_1),\ldots,(a_{n-|C_0|},l_{n-|C_0|})\} \cup (\mathbf{x}_{C_0} \times \{{0}\})) \nonumber \\ & \qquad\qquad \sum_{\substack{C_1,\ldots,C_{n-|C_0|} \\ \cup_j \{C_j\} = \{1,\ldots,n\}\setminus C_0}} \prod_{j=1}^{n-|C_0|} m(a_j)g_Y(a_j,l_j)\mathbf{1}_{[a_j, a_j + l_j]}(x_{C_j})\,{\mathrm{d}} a_j\,{\mathrm{d}} l_j \prod_{i=1}^n{\mathrm{d}} x_i. \end{align}

Similarly,

\begin{align*} & {\mathbb{E}}[{\mathbf{1}}_F(U){\mathbb{P}}(X \in A \mid U )] \\ & = \sum_{n=0}^{\infty}\frac{{\mathrm{e}}^{-\ell(\mathcal{X})}}{n!} \\ & \qquad \int_{\mathcal{X}^n} p(\{x_1, \ldots, x_n\})\sum_{C_0 \subset\{1,\ldots,n\}} \frac{1}{(n - |C_0|)!}\prod_{i \in C_0} w_{x_i} \\ & \qquad\quad \int_{(\mathbb{R}\times\mathbb{R}_0^+)^{n-|C_0|}}\mathbf{1}_F(\{(a_1,l_1),\ldots,(a_{n-|C_0|},l_{n-|C_0|})\} \cup (\mathbf{x}_{C_0} \times \{{0}\})) \\ & \qquad\qquad \mathbb{P}(X \in A\mid U = \{(a_1,l_1),\ldots,(a_{n-|C_0|},l_{n-|C_0|})\} \cup (\mathbf{x}_{C_0} \times \{{0}\})) \\ & \qquad\qquad \sum_{\substack{C_1,\ldots,C_{n-|C_0|} \\ \cup_j \{ C_j\} = \{1,\ldots, n\} \setminus C_0}} \prod_{j=1}^{n-|C_0|} m(a_j)g_Y(a_j,l_j)\mathbf{1}_{[a_j, a_j + l_j]}(x_{C_j})\,{\mathrm{d}} a_j\,{\mathrm{d}} l_j\prod_{i=1}^n{\mathrm{d}} x_i.\end{align*}

Next, we plug in the expression for $\mathbb{P}(X \in A\mid U=\mathbf{u})$ proposed in the statement of the theorem. Upon substitution, changing integration order, and rearranging, we obtain

\begin{align*} & {\mathbb{E}}[{\mathbf{1}}_F(U){\mathbb{P}}(X \in A \mid U )] \\ & = \sum_{n=0}^{\infty}\frac{{\mathrm{e}}^{-\ell(\mathcal{X})}}{n!} \\ & \qquad \sum_{C_0 \subset\{1,\ldots,n\}}\frac{1}{(n - |C_0|)!} \\ & \qquad\quad \int_{\mathcal{X}^n} p(\mathbf{y} \cup \mathbf{x}_{C_0})\mathbf{1}_A(\mathbf{y}\cup\mathbf{x}_{C_0}) \prod_{i \in C_0} w_{x_i}\end{align*}

\begin{align*} & \qquad\qquad \int_{(\mathbb{R}\times\mathbb{R}_0^+)^{n-|C_0|}} \mathbf{1}_F(\{(a_1,l_1),\ldots,(a_{n-|C_0|},l_{n-|C_0|})\} \cup (\mathbf{x}_{C_0} \times \{{0}\})) \\ & \qquad\qquad\quad c(\{(a_1,l_1),\ldots,(a_{n-|C_0|},l_{n-|C_0|}) \cup (\mathbf{x}_{C_0}\times\{0\})) \\ & \qquad\qquad\quad \sum_{\substack{D_1,\ldots,D_{n-|C_0|} \\ \cup_j \{D_j\} = \{1,\ldots,n\} \setminus C_0}} \prod_{j=1}^{n-|C_0|} m(a_j)g_Y(a_j,l_j)\mathbf{1}_{[a_j,a_j + l_j]}(y_{D_j}) \\ & \qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad \Bigg(\int_{\mathcal{X}^{n-|C_0|}}p(\mathbf{x}) \sum_{\substack{C_1,\ldots,C_{n-|C_0|} \\ \cup_j \{C_j\} = \{1,\ldots,n\} \setminus C_0}} \\ & \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad \prod_{j=1}^{n-|C_0|}\mathbf{1}_{[a_j, a_j+l_j]}(x_{C_j}) \prod_{j\not\in C_0} {\mathrm{d}} x_j \Bigg) \\ & \qquad\qquad\quad \prod_{j=1}^{n-|C_0|}{\mathrm{d}} a_j\,{\mathrm{d}} l_j\prod_{i\in C_0}{\mathrm{d}} x_i\prod_{k=1}^{n-|C_0|}{\mathrm{d}} y_k = \mathbb{E}[\mathbf{1}_F(U)\mathbf{1}_A(X)],\end{align*}

since the term within brackets cancels out against the normalisation constant $c{(}{\cdot}{)}$ .

Strikingly, although the marking mechanism is more complicated than that in [Reference Van Lieshout and Markwitz25], the conditional distribution of X has the same form.

The conditional distribution of W can be obtained in the same vein, by considering $\mathbf{1}_{A}(W)$ instead of $\mathbf{1}_A(X)$ for A a Borel set in $\mathcal{N}_{\mathcal{X}\times(\mathbb{R}\times\mathbb{R}^+)}$ , the space of marked point configurations. It is given by

\begin{align*} & \mathbb{P}(W \in A \mid U=\mathbf{u}) \\ & \propto \int_{\mathcal{X}^{n-m}}p(\{a_1,\ldots,a_m,x_{1},\ldots,x_{n-m}\}) \\ & \qquad \mathbf{1}_A(\{(a_1,(a_1,0)),\ldots,(a_m,(a_m,0)),(x_{1},(a_{m+1},l_{m+1})),\ldots,(x_{n-m},(a_{n},l_{n}))\}) \\ & \qquad \prod_{i=1}^{n-m}\mathbf{1}_{[a_{m+i}, a_{m+i} + l_{m+i}]}(x_{i})\,{\mathrm{d}} x_i .\end{align*}

4.1. Non-homogeneous point process densities

We look first at inhomogeneity that manifests itself via the occurrence time distribution. In view of (3.1), it is natural to add inhomogeneity by means of the first-order interaction function $\phi_1$ , a procedure known as type I inhomogeneity [Reference Jensen and Nielsen13]. The idea is to let $\phi_1(\{ x \}) = \beta(x)$ vary over time according to a measurable function $\beta$ that maps $x \in \mathcal{X}$ to $[0,\infty)$ . In many applications, it may make sense to model $\beta$ as a step function. More specifically, given a measurable partition $B_1, \ldots, B_K$ of $\mathcal{X}$ , set $\beta(x) = \sum_{k=1}^{K}\beta_k\mathbf{1}_{B_k}(x)$ , $x \in \mathcal{X}$ , where $\beta_k \geq 0$ is the value that $\beta$ takes in the corresponding set $B_k$ .

The function $\phi_1$ can be combined with classic second- and higher-order interaction functions. For instance, the probability density of the non-homogeneous area-interaction point process [Reference Baddeley and van Lieshout3] becomes

(4.1) \begin{align} p(\mathbf{x}) = \alpha_p \Bigg( \prod_{x \in \mathbf{x}} \beta(x) \Bigg) \exp[{-}(\log \gamma)\, \ell( {\mathcal{X}} \cap U_r(\mathbf{x}) )] \end{align}

with respect to a unit-rate Poisson process on $\mathcal{X}$ . The parameter $\gamma$ quantifies the interaction strength, r the radius of interaction, and $\alpha_p = c(\beta({\cdot}), \gamma)$ is a normalisation constant [Reference Baddeley and van Lieshout3] that depends on the function $\beta$ as well as on $\gamma$ . Additionally, $U_r(\mathbf{x}) =\bigcup_{i=1}^n B(x_i, r)$ , where $B(x_i,r)$ is the closed interval $[ x_i-r, x_i + r]$ . We observe regularity for $\gamma \lt 1$ and clustering for $\gamma \gt 1$ , and $\gamma = 1$ corresponds to a non-homogeneous Poisson process with intensity function $\beta$ . For further examples, we refer to [Reference Van Lieshout24].

4.2. Parametric modelling of the mark kernel

To proceed, parametric forms for $G_Y$ and m must be developed. We begin by modelling $G_Y$ , the semi-Markov kernel that determines the length of time until the next jump. We may take one of the time-dependent probability density functions considered in Section 2.2. For instance, $g_Y(a, \cdot)$ could be the density function of an exponential distribution with rate

(4.2) \begin{align} \lambda(a;\ \alpha) = \alpha (b + \sin(c a)), \quad a \in \mathbb{R}, \end{align}

where c specifies the period and $b \geq 1$ the elevation away from 0. The parameter $\alpha$ determines the amplitude of the harmonic.

We could proceed in a similar fashion for $g_Z$ . However, there are two problems with such an approach. From a probabilistic point of view, tractable expressions for the renewal density m in terms of the semi-Markov kernels $G_Y$ and $G_Z$ do not seem to exist, and, statistically speaking, lengths of Z-phases cannot be observed. Therefore, we shall model m directly. The following proposition justifies this approach.

Proof. As $0 \leq {h(t)} / {\tilde{m}(t)} \leq 1$ , we may use a time-dependent thinning approach with retention probability $p(t) = {h(t)} / {\tilde{m}(t)}$ . Algorithmically, the sought-after process can be constructed as follows. Initialise $\widehat{S}_0 = 1$ , $\widehat{X}_0 = 0$ , and $\widehat{X}_1= X_1$ . Also, set $\widehat{S}_{2i} = 1$ , $\widehat{S}_{2i-1} = 0$ for $i\in{\mathbb{N}}$ , and $j=1$ . For each jump time $X_i$ , $i=1, 2, \ldots$ ,

• • with probability $p( X_{i})$ , if j is even, update $\widehat{X}_{j+1} = X_{i+1}$ and increment j by 1; for odd j update $\widehat{X}_{j+1} = \widehat{X}_j$ , $\widehat{X}_{j+2} = X_{i+1}$ , and increment j by 2;

• • otherwise, if j is odd, update $\widehat{X}_{j+1} = X_{i+1}$ and increment j by 1; for even j update $\widehat{X}_j = X_{i+1}$ leaving j unchanged.

Because complete cycles correspond to intervals in between accepted points $X_i$ , $i=0, 1,2,\ldots$ ,

\begin{align*}H(t) = \sum_{n=1}^\infty {\mathbb{P}}( \widehat{X}_{2n} \leq t ) =\sum_{n=1}^\infty {\mathbb{P}}( X_n \leq t, X_n\,{\mathrm{retained}}) =\int_0^t \frac{h(s)}{\tilde{m}(s)}\,{\mathrm{d}} \tilde M(s).\end{align*}

We can therefore conclude that the intensity of the thinned process is $[{h(t)}/{\tilde{m}(t)}] \tilde{m}(t) = h(t)$ (see, e.g., [Reference Daley and Vere-Jones7, pp. 78–79]) and hence $(\widehat S_n, \widehat X_n)_{n=1}^\infty$ is a non-homogeneous alternating renewal process that satisfies the proposed conditions.

As an illustration, suppose that m is a step function,

(4.3) \begin{align} m(t) = \sum_{j=1}^{J} \delta_j\mathbf{1}_{A_j}(t), \quad t\in \mathbb{R}^+, \end{align}

that takes J different values $\delta_j \gt 0$ on bounded Borel sets $A_j$ ( $j=1,\ldots, J$ , with $J\in{\mathbb{N}}$ ) and 0 elsewhere. The following corollary lays out conditions under which $h=m$ is the renewal density of an alternating renewal process whose Y-phases are governed by (4.2).

As noted before, in practice, the starting point 0 is moved back to $-\infty$ . Realisations $\mathbf{u}$ from the specified model may be obtained as follows. First, a set of points $\mathbf{x} \subset \mathcal{X}$ in time is chosen according to the probability density function $p({\cdot})$ by, for example, coupling from the past [Reference Kendall15] or the Metropolis–Hastings algorithm [Reference Geyer9]. Next, for each point $x\in\mathbf{x}$ , it is determined whether or not it is an atom based on $w_x$ , see (2.5). If this is not the case, we appeal to Proposition 2.5 and use rejection sampling with a proposal distribution that simulates a uniformly distributed point in $A_j \cap ({-}\infty, x]$ chosen with probability $\delta_j \ell(A_j\cap({-}\infty,x]) / \sum_{i =1}^J \delta_i \ell(A_i \cap({-}\infty, x])$ and acceptance probability $\exp{[}{-}\lambda(a;\ \alpha) (x - a) ]$ . The result is a sample a from $f_x(a)$ , cf. (2.7). The length is then sampled according to an exponential distribution with parameter $\lambda(a; \alpha)$ shifted by $x-a$ (see (2.8)). It is interesting to observe that, in contrast to the alternating renewal case studied in [Reference Van Lieshout and Markwitz25], using the marginal distribution with respect to A and then the conditional given A is computationally simpler than sampling L first.

4.3. Statistical aspects

In practical applications, both the family of probability density functions $g_Y(t, \tau;\ \theta)$ for the sojourn times in phase Y and the function $m(t;\ \xi)$ rely on unknown parameters $\eta = (\theta, \xi)$ that must be estimated. The log-likelihood $L(\eta;\ \mathbf{u})$ follows directly from (3.3). Upon observing $\mathbf{u} = \{ (a_1, 0), \ldots, (a_m, 0),(a_{m+1}, l_{m+1}),\ldots, (a_n, l_n) \}$ ,

(4.4) \begin{equation} L(\eta;\ \mathbf{u}) = \sum_{i=1}^{m} \log\bigg(1 \!-\! \int_{-\infty}^{a_i}[1\!-\!G_Y(s,a_i-s;\ \theta)]m(s;\ \xi)\,{\mathrm{d}} s\bigg)\! +\! \sum_{i=m+1}^n \log(m(a_i;\ \xi)g_Y(a_i,l_i;\ \theta)).\end{equation}

When the semi-Markov kernels $G_Y$ and $G_Z$ , as well as the renewal density $m \equiv ( \mathbb{E} Y +\mathbb{E} Z )^{-1}$ , do not vary in time, (4.4) reduces to the renewal likelihood in [Reference Van Lieshout and Markwitz25].

We illustrate the procedure by means of a specific example. For the sojourn times, we take an exponential model; for the function m, we use (4.3). Assume that $G_Y(t, \cdot)$ , $t\in \mathbb{R}$ , is the cumulative distribution function of an exponential distribution with rate parameter $\lambda(t)$ as in (4.2), so that $G_Y(t, l;\ \lambda(t)) = 1 - {\mathrm{e}}^{-l \lambda(t)}$ . In the homogeneous case where $\lambda(t) \equiv \alpha \gt 0$ (that is, $b=1$ and $c=0$ ), $J = 1$ , $A_1 = \mathbb{R}$ , and $0 \leq \delta_1 \leq \alpha$ , we obtain

(4.5) \begin{equation} \int_{-\infty}^{t}[1 - G_Y(s, t-s)]m(s)\,{\mathrm{d}} s = \int_{-\infty}^{t}{\mathrm{e}}^{-\alpha(t-s)}\,{\mathrm{d}} s = \frac{1}{\alpha} = \mathbb{E}[Y],\end{equation}

where $\mathbb{E}[Y]$ refers to the expected length of a Y-phase in a homogeneous alternating renewal process. Using this result, we obtain $f_Y(l;\ \alpha) = g_Y(t,l;\ \lambda(t))$ , where $f_Y$ is the probability density function of the length of such a Y-phase (see, e.g., [Reference Van Lieshout and Markwitz25, Theorem 2.2]). Following through, we additionally obtain

\begin{equation*} L(\alpha,\delta_1;\ \mathbf{u}) = m\log\bigg(1-\frac{\delta_1}{\alpha}\bigg) + (n-m)\log\delta_1 + (n-m)\log\alpha - \alpha \sum_{i=m+1}^n l_i,\end{equation*}

which corresponds exactly to the simplified version of [Reference Van Lieshout and Markwitz25, (7)], the log-likelihood in the homogeneous case.

Returning to the non-homogeneous case, the atom probability for a given time $x \in \mathcal{X}$ is

\begin{equation*} w_x = 1 - \sum_{j=1}^{J}\delta_j\int_{A_j \cap ({-}\infty, x]}{\mathrm{e}}^{-(x-s)\lambda(s)}\,{\mathrm{d}} s.\end{equation*}

The likelihood equation (4.4), after substitution and discarding of terms that do not depend on the parameters, becomes

\begin{align*} L(\delta,\alpha;\ \mathbf{u}) & = \sum_{i=1}^{m}\log\Bigg(1-\sum_{j=1}^{J}\delta_j\int_{(a_i-A_j)\cap[0,\infty)}{\mathrm{e}}^{-\alpha r(b+\sin(ca_i-cr))}\,{\mathrm{d}} r\Bigg) \\ & \quad + \sum_{i=m+1}^n\log\Bigg(\sum_{j=1}^{J}\delta_j\mathbf{1}_{A_j}(a_i)\Bigg) + (n-m)\log\alpha-\alpha\sum_{i=m+1}^{n}l_i(b + \sin(ca_i)).\end{align*}

The resulting equations can be solved numerically to find optimal values for $\delta_k$ , $k = 1, \ldots, J$ , and $\alpha$ under the inequality constraints $0 \leq \delta_j \leq \alpha(b-1)$ , $j=1, \ldots, J$ .

The distribution of unobserved occurrence times may also be considered as a parameter to be estimated using the reported intervals. To do so, since the form of the conditional distribution of W given U according to Theorem 3.1 is identical to that for alternating renewal process-based censoring, the simulation techniques developed in [Reference Van Lieshout and Markwitz25] to obtain realisations of the marked occurrence times given a sample $\mathbf{u}$ of U apply. A brief overview is provided below.

Let $\zeta$ be the parameter vector corresponding to the area-interaction process parameters $(\beta({\cdot}), \gamma)$ . Write $p(\mathbf{x)}$ , the probability density from (4.1), as $p(\mathbf{x};\ \zeta) = c(\zeta)h(\mathbf{x};\ \zeta)$ , where c refers to the normalisation constant as before and h is the unnormalised density. Reference parameter values $\zeta_0$ are picked using a Monte Carlo EM approach [Reference Geyer8], and the likelihood given the interval set $U = {\mathbf{u}}$ , for N samples, is estimated by

\begin{equation*} l_N(\zeta) = \log\Bigg(\frac{1}{N}\sum_{i=1}^N\frac{h(X_{{\mathbf{u}},i};\ \zeta)}{h(X_{{\mathbf{u}},i};\ \zeta_0)}\Bigg) - \log\Bigg(\frac{1}{N}\sum_{i=1}^N\frac{h(X_{i};\ \zeta)}{h(X_{i};\ \zeta_0)}\Bigg),\end{equation*}

where $X_i$ are samples from the area-interaction process, and $X_{\mathbf{u}, i}$ are samples from the conditional distribution of occurrence times given $U = {\mathbf{u}}$ [Reference Van Lieshout and Markwitz25]. The likelihood is then maximised to find the optimal value(s) of $\zeta$ [Reference Geyer8].

Once the area-interaction parameters have been estimated, a Metropolis–Hastings algorithm [Reference Brooks, Gelman, Jones and Meng5, Reference Meyn and Tweedie18, Reference Møller and Waagepetersen19] for a fixed number of points can be used to determine the most likely location of the point within each interval, given the parameter values. For further details and conditions under which the algorithm converges to the desired distribution, we refer to [Reference Van Lieshout and Markwitz25, Propositions 4.3–4.5].

5.1. Model misspecification

The first source of inhomogeneity in our model is the semi-Markov kernel $G_Y(a, \cdot)$ which determines the time until the next jump for starting point $a\in \mathbb{R}$ . In the homogeneous case, as in (4.5), these may be assumed to be independent and identically Gamma or Weibull distributed [Reference Van Lieshout and Markwitz25]. In the non-homogeneous case, the density function is dependent on a.

We have devised the following experiment to illustrate the effect of using a homogeneous model erroneously. We generate a realisation of intervals using the non-homogeneous alternating renewal process model, where the interval censoring mechanism is governed by a Weibull distribution with shape parameter $k=1$ and rate parameter $\lambda_Y(t;\ \alpha) = \alpha(1.6+\sin(2\pi t))$ (see (4.2)) for $\alpha = 1$ . Regarding the other model ingredients, we assume that $p({\cdot})$ is of the form given in (4.1) with $\beta = 400$ and $\gamma = 1$ , i.e. a homogeneous Poisson process on $\mathcal{X}$ with intensity 400. We additionally set $m(t) = 0.6\, \mathbf{1}_{ [-0.2, 1) }(t)$ .

After generating these intervals using the non-homogeneous model, we now assume, wrongly, that these intervals were instead generated by a homogeneous interval censoring scheme. Specifically, we fit a Weibull distribution with parameters $k \gt 0$ and constant rate $\lambda_Y(t;\ \alpha) \equiv \alpha$ for $\alpha \gt 0$ using maximum likelihood estimation (see, e.g., [Reference Van Lieshout and Markwitz25, Section 6.1.1]) and we obtain parameter estimates $\widehat{k} = 0.9$ and $\widehat{\alpha} = 2.0$ .

We plot the survival time densities for both sets of parameters. We must choose an exact point in time to obtain the Weibull parameters for the non-homogeneous model, so we assume that we are looking at time $t=0.6$ . This plot, for both models, can be seen in Figure 2. Compared to the actual model, the homogeneous model is able to roughly discern the shape of the distribution, but struggles with the scale. A homogeneous model would have generated more intervals shorter than about $0.5$ , and fewer of longer length.

Figure 2. The unbroken line corresponds to the actual probability density of interval length for $k=1$ and $\lambda(0.6;\ 1) = 1$ . The dotted line corresponds to the estimated survival time density.

5.2. Inhomogeneity in renewal density and survival time

In our second experiment, we add inhomogeneity in m to the the model and study the effect on $f_x$ , cf. (2.7). As in Section 5.1, consider an exponential semi-Markov kernel density $g_Y$ with rate parameter either constant, $\lambda(t;\ \alpha) = 1.3\alpha$ , or varying in time according to $\lambda(t;\ \alpha) = \alpha(1.3 + \sin(2\pi t))$ . Furthermore, set $m(t) = 0.4$ for $t \in [{-}0.2, 1)$ in the constant case, and

\begin{align*}m(t) = \begin{cases} 0.4 & t \in [{-}0.2, 0.4) \\ 0.1 & t \in [0.4, 1) \end{cases}\end{align*}

in the time-varying case. We set $\alpha = 1.6$ , so the largest value of $\delta_i$ which guarantees that $g_Y$ is the Radon-Nikodym derivative of a semi-Markov kernel is $1.6 \times ( 1.3 - 1 ) = 0.48$ . Figure 3 shows the graphs of $f_x({\cdot})$ for the four possible combinations of $g_Y$ and m obtained from 200,000 samples from $q_x$ for $x=1$ . In Figures 3a and 3b, we assume that $\lambda$ is constant. When m is also constant as in Figure 3a, the marginal distribution of the starting times, by Proposition 2.5, is a shifted exponential distribution. If m is allowed to vary in time, the exponential curve is broken at $t=0.4$ , the discontinuity point of m, resulting in a zigzag pattern. In both Figures 3a and 3c, m is constant, but in Figure 3c the rate parameter of $g_Y$ varies according to a harmonic. The resulting sinusoidal modulation is clearly visible. Finally, allowing m to vary too results in a break at its discontinuity point $t=0.4$ as seen in Figure 3d. We conclude that the non-homogeneous model is able to capture various types of inhomogeneity well and is far more flexible than the homogeneous one.

Figure 3. Probability density function of the starting time $f_x({\cdot})$ with $x=1$ for various choices of $g_Y$ and m.

5.3. Inhomogeneity in occurrence time distribution

In the previous subsections, we have assumed that the first-order interaction function $\beta$ of the point process X of occurrence times remains constant over the entire sampling window (0, 1). In our final example, we relax this assumption in that we consider a ‘peak time’ in which events are more likely to occur and investigate the effect on the conditional distribution of occurrences. More precisely, we take an area-interaction model (4.1) with

\begin{align*}\beta(y) = \begin{cases} 3, & y \in [0, c_1), \\ 5, & y \in [c_1, c_2), \\ 3, & y \in [c_2, 1], \end{cases} \end{align*}

and critical range $[c_1, c_2) = [0.81, 0.85)$ . The radius of interaction is set to $r=0.1$ and we consider both a regular $(\eta = -1.2)$ and a clustered $(\eta = 1.2)$ model.

As in [Reference Van Lieshout and Markwitz25], consider the set $\mathbf{u} = \{(0.45, 0.4),(0.51, 0), (0.58, 0)\}$ that contains one non-degenerate interval. Recall that the entries are parametrised as (a, l), where a is the starting point and l is the length. Figure 4 plots the conditional distribution of the occurrence time on the interval $[0.45, 0.85]$ given ${\mathbf{u}}$ for the regular and clustered models. To create this figure, a Metropolis–Hastings algorithm [Reference Van Lieshout and Markwitz25, Algorithm 4.2] was run for 600 000 time steps, with the first 100 000 iterations thrown out due to burn-in.

Figure 4. A comparison between a regular and clustered model with a ‘peak time’ added by changing the intensity function within a critical range.

The general shape of the graphs is similar to the corresponding plots for constant $\beta = 3$ in [Reference Van Lieshout and Markwitz25, Figures 2 and 3]. For the clustered model, the occurrence time is more likely to happen close to the atoms; for regular models the probability density is shifted away from the atoms. In the non-homogeneous case, the higher value of $\beta$ during the peak times causes a clear bump in the range $[c_1, c_2) =[0.81, 0.85)$ , demonstrating the ability of the non-homogeneous model to favour certain occurrence times over others.

The previous three experiments show the effects that the added non-homogeneity may have on the model. It is able to draw from a more complicated interval censoring scheme to generate intervals and provide a more realistic starting time distribution based on the semi-Markov kernel and Markov renewal function, while having the functionality of splitting the observation window into regions of differing likelihood, and deal with more complicated interaction structures in the underlying point process. For a concrete application of the non-homogeneous model to a crime dataset, we direct the reader to [Reference Markwitz17]. An application of the homogeneous model can be found in [Reference Van Lieshout and Markwitz25, Section 6].