1.1 Termination time \(T\)

Let the random variable \(T > 0\) represent the termination time at which an individual ceases its activity. Either by death or by capture in a trap.

We can describe the distribution of \(T\) using its:

• Probability Density Function (PDF): \(f(t) = \lim_{\d{t} \to 0} \frac{\p(t < T \leq t + \d{t})}{\d{t}} = F'(t)\)

• Cumulative Distribution Function (CDF): \(F(t) = \p(T \leq t) = \int_0^t f(u) \d{u}\),

• Survival function: \(S(t) = \p(T > t) = 1 - F(t)\)

• Hazard Rate function: \(h(t) = \lim_{\d{t} \to 0} \frac{\p(t < T \leq t + \d{t} \mid T > t)}{\d{t}}\)

The hazard function can be interpreted as the instantaneous rate of events (i.e. the expected number of terminations per unit of time). It is often more convenient to express hypotheses in terms of the hazard function and then derive the other expressions from it. The following relationships can be useful.

Since \(\p(t < T \leq t + \d{t} \mid T > t) = \frac{\p(t < T \leq t + \d{t}, T > t)}{\p(T > t)} = \frac{\p(t < T \leq t + \d{t})}{S(t)},\) it follows from the definition that,

\[\begin{equation} \tag{1.1} h(t) = \frac{f(t)}{S(t)} = -\frac{S'(t)}{S(t)} = - \frac{\d}{\d{t}} \log S(t) \end{equation}\]

Reciprocally, given a hazard function, survival can be obtained using

\[\begin{equation} \tag{1.2} S(t) = \exp \Big( - \int_0^t h(u) \d{u} \Big), \end{equation}\]

assuming that \(S(0) = 1\).

1.2 Termination cause \(C\)

Let the categorical random variable \(C \in \{0, 1, \ldots, I\}\) be an indicator of the termination cause. Either death (\(C = 0\)) or capture in trap \(i\) (\(C = i\)).

Its Probability Mass Function is:

\[\begin{equation} \tag{1.3} \p(C = i) = \omega_i, \end{equation}\]

with \(0 \leq \omega_i \leq 1\) and \(\sum \omega_i = 1\).

Due to the dispersal process, the relative risks of being caught in either trap will shift over time. Thus, it will be easier to specify these risks conditional on time and then average over the PDF of T:

\[\begin{equation} \tag{1.4} \omega_i = \int_0^\infty \p(C = i \mid T = t) f(t) \d{t}. \end{equation}\]

1.3 Specific termination times \(T_i\)

Let \(T_i\) be the termination time in a scenario where the only possible termination cause is \(i\).

\(T\) and \(T_i\) are different random variables, with different PDFs, CDFs, survivals and hazard rates, for all of which we define the corresponding specific versions \(f_i\), \(F_i\), \(S_i\) and \(h_i\).

The relationship between these variables is that \[\begin{equation} \tag{1.5} T = \min\{T_i:\; i = 0, \ldots, I\} \end{equation}\]

We can think of it as if at every instant, all specific termination causes were tested independently. The first that occurs determines both the termination time \(T\) and cause \(C\).

1.4 Position \(X\)

Finally, let \(X(t) \in \mathbb{R}^2\) be the random location of an individual at a given time \(0 \leq t \leq T\). Without loss of generality, we assume that \(X(0) = (0, 0)\).

2.1 Fixed and known position

2.2 Different constant conditional hazard rates imply non-constant overall hazard rate

Assume a fixed position (\(x_0\)) of an individual, and no mortality.

\[\begin{equation} \begin{aligned} \tag{2.1} \p(t < T \leq t + \d{t} \mid T > t) & = \frac{\p(t < T \leq t + \d{t})}{\p(T > t)} \\ & = \frac{ \sum \p(t < T \leq t + \d{t} \mid C = i) \p(C = i) }{ \p(T > t \mid C = i) \p(C = i) } \end{aligned} \end{equation}\]

Letting \(\omega_i = \p(C = i)\), assuming constant conditional hazard rates \(\lambda_i\) implies \(f_i(t \mid C = i) = \lambda_i e^{-\lambda_i t}\), thus:

\[\begin{equation} \begin{aligned} \tag{2.2} h(t) & = \lim_{\d{t} \to 0} \frac{\p(t < T \leq t + \d{t} \mid T > t)}{\d{t}} = \lim_{\d{t} \to 0} \frac{\p(t < T \leq t + \d{t})}{S(t) \d{t}}\\ & = \frac{\sum f_i(t \mid C = i) \omega_i}{\sum S(t \mid C = i) \omega_i} = \frac{\sum \omega_i \lambda_i e^{-\lambda_i t}}{\sum \omega_i e^{-\lambda_i t}} \\ & = \sum \lambda_i \theta_i(t), \end{aligned} \end{equation}\] with \(\theta_i(t) = \frac{\omega_i e^{-\lambda_i t}}{\sum \omega_j e^{-\lambda_j t}}\).

The overall hazard rate is thus a weighted average of conditional hazard rates, with coefficients that evolve in time.

The coefficients \(\theta_i(t)\) are a weighted version of the softmax transformation of \(\{\lambda_i\}\) with temprature \(t\). They verify \(\sum \theta_i(t) = 1,\; \forall t > 0\), \(\theta_i(0) = \omega_i / \sum w_j\) and \(\lim_{\d{t} \to \infty}\theta_i(t) = 0\) except for \(i\) such that \(\lambda_i\) is smallest.

Thus, the coefficients start by reproducing the distribution of the termination causes \(\omega_i\) and progressively concentrate on the trap with lower hazard rate.

2.3 The conditional hazard rates must be all equal to the overall hazard rate

The conditional hazard rates are the hazard rates of capture events in a specific trap, conditional to capture in that trap.

\[ \begin{aligned} h(t \mid C = i) & = \lim_{\d{t} \to 0} \frac{\p(t < T \leq t + \d{t} \mid C = i,\, T > t)}{\d{t}} \\ & = \lim_{\d{t} \to 0} \frac{\p(t < T \leq t + \d{t}, C = i \mid T > t)}{P(C = i \mid T > t) \d{t}} \\ \end{aligned} \]

The key observation is that since the position of the individual is fixed, the scenario is constantly the same and thus \(P(C = i \mid T > t) = P(C = i), \; \forall t\). So, the termination cause \(C\) is independent from the termination time \(T\), and we have, \[ \begin{aligned} \p(t < T \leq t + \d{t}, C = i \mid T > t) & = \p(t < T \leq t + \d{t} \mid T > t) \p(C = i \mid T > t) \\ & = \p(t < T \leq t + \d{t} \mid T > t) \p(C = i) \\ & = \omega_i\, \p(t < T \leq t + \d{t} \mid T > t) \end{aligned} \]

Thus, \[\begin{equation} \tag{2.3} \begin{aligned} h(t \mid C = i) & = \lim_{\d{t} \to 0} \frac{\omega_i \p(t < T \leq t + \d{t} \mid T > t)}{\omega_i \d{t}} \\ & = \lim_{\d{t} \to 0} \frac{\p(t < T \leq t + \d{t} \mid T > t)}{\d{t}} \\ & = h(t), \quad \forall i \end{aligned} \end{equation}\]

Alternatively, one could also argue that the overall hazard rate should be constant since the situation with fixed position remains unchanged.

Let \(h(t) = \lambda_0\). From (2.2), we have

\[\begin{equation} \sum \lambda_i \theta_i(t) = \lambda_0. \end{equation}\]

Which is only possible if \(\lambda_i = \lambda_0, \; \forall i\).

This is confirmed with an empirical simulation with 2 traps and fixed positions as displayed in figure 2.1).

Figure 2.1: Conditional and overall empirical hazard rates for a scenario with fixed position of individuals and 2 traps at different distances.

2.4 The overall hazard rate is the sum of specific hazard rates

Let \(\lambda_i = h_i(t), \; i = 0,\ldots, I\) the specific hazard rates associated to death (\(i = 0\)) and capture in trap \(i\). This is, the (constant) hazard functions of the corresponding events if the specific cause was the only one at play. In other words, if the individual was active in its fixed position until it was eventually terminated by \(i\).

Since \(T = \min\{T_i: i = 0, \ldots, I\}\), in an instant \(\d{t}\), \[\begin{equation} \begin{aligned} \p(t < T \leq t + \d{t} \mid T > t) & = 1 - \prod_{i = 0}^I (1 - \p(t < T_i \leq t + \d{t} \mid T > t)) \\ & = 1 - \prod_{i = 0}^I (1 - \lambda_i \d{t}) \\ & = 1 - \big[ 1 - \lambda_0 \d{t} - \cdots - \lambda_I \d{t} + O(\d{t}^k) \big], \quad k \geq 2 \\ & = \d{t} \sum_{i = 0}^I \lambda_i + O(\d{t}^k) \end{aligned} \end{equation}\]

Thus, the overall hazard rate is: \[\begin{equation} \tag{2.4} \begin{aligned} h(t) & = \lim_{\d{t} \to 0} \frac{\p(t < T \leq t + \d{t} \mid T > t)}{\d{t}} \\ & = \lim_{\d{t} \to 0} \sum_{i = 0}^I \lambda_i + O(\d{t}^{k-1}), \quad k \geq 2 \\ & = \sum_{i = 0}^I \lambda_i \end{aligned} \end{equation}\]

3.1 With only one trap

Only one trap means that \(I = 1\) and \(C\) can only take the value \(0\) in case of death prior to capture, or the value \(1\) in the converse.

We assume that, if an individual stays at the same location, the hazard function of the capture time \(T \mid C = 1, X(t) = x_0\) is constant in time, and decreases exponentially with the distance \(r = \norm{x_0} > 0\) from the trap:

\[ h_{1r}(t) = h(t \mid C = 1, X = x_0) = \lambda_r = \alpha e^{-\beta r}, \quad \alpha, \beta > 0. \]

From equation (1.2), it follows that:

\[ S(t \mid C = 1, X = x_0) = e^{-\lambda_r t} \]

Hence, \(f(t \mid C = 1, X = x_0) = -S'(t \mid C = 1, X = x_0) = \lambda_r e^{-\lambda_r t}\), which is the PDF of the exponential distribution of rate \(\lambda_r\).

\[\begin{equation} \tag{3.1} T \mid C = 1, \norm{X} = r \sim \mathrm{Ex}\big( \alpha e^{-\beta r} \big) \end{equation}\]

We can also compute the daily capture probability:

\[ \p(0\le T\le1 \mid C = 1, X = x_0) = \int_0^1\lambda_re^{-\lambda_r t} \d{t} = 1- e^{-\lambda_r} \]

3.2 With several traps

Consider now the case where we have \(I > 1\) traps in different locations. The aim is to compute the probability of capture where these traps are active together and to model the capture time of each trap.

Here, the difficulty comes from the dependence between these variables because if an individual is trapped by the trap \(i\), then the capture time of the other ones will be equal to infinity.

This is why we consider the global termination time \(T\) jointly with the termination cause \(C\) described initially.

4.1 With only one trap

All the previous section was about the computation of the probability of capture for a constant location. Here, we would like to extend the expression assuming the location is random.

Consider the hazard function, we have, by the definition that :

\[ h(t) = \lim_{\d{t} \rightarrow0} \frac{\p(t \le T \le \d{t} |T \ge t)}{\d{t}}= \lim_{\d{t} \rightarrow0} \int_{\mathbb{R}^2}\frac{\p(t \le T \le \d{t} |T \ge t, X_t=x)}{\d{t}}f_{X_t}(x) \d{x} \]

As everything is positive, we can switch between limit and integral :

\[ h(t) = \int_{\mathbb{R}^2} \lim_{\d{t} \rightarrow0}\frac{\p(t \le T \le t + \d{t} |T \ge t, X_t=x)}{\d{t}}f_{X_t}(x)\d{x} \]

More, the limit allows us to consider this quantity as the hazard function of the capture time with a constant location. So we can write this :

\[ h(t) =\int_{\mathbb{R}^2} \lambda(x) f_{X_t}(x)\d{x} \]

where \(\lambda(x) = \alpha e^{-\beta ||x-x_0||^2}\) where \(x_0\) is the location of the trap.

Computation of the integral

We would like to calculate :

\[ \int_{\mathbb{R}^2}\alpha e^{-\beta ||x-x_0||^2}f_{X_t}(x)\d{x} = \alpha \mathbb{E}[e^{-\beta||X_t-x_0||^2}] \]

We have :

\[\begin{align*} \mathbb{E}[e^{-\beta||X_t-x_0||^2}] & = \int_{\mathbb{R}^2} e^{-\beta ||x-x_0||^2} \frac{1}{8\pi Dt}e^{-\frac{||x||^2}{8Dt}}\d{x} \\ & = \frac{1}{8\pi Dt}\int_{\mathbb{R}^2} e^{-(\frac{||x||^2}{8Dt}+\beta ||x-x_0||^2)}\d{x} \end{align*}\]

But the term in the exponential can be written as :

\[\begin{align*} \frac{||x||^2}{8Dt}+\beta ||x-x_0||^2 & = \frac{||x||^2}{8Dt}+ \beta(||x||^2 - 2 <x|x_0> + ||x_0||^2) \\ & = \gamma_t(||x||^2 - 2 <x|\frac{\beta}{\gamma_t}x_0>)+ \beta||x_0||^2 \\ & = \gamma_t(||x-\frac{\beta}{\gamma_t}x_0||^2-||\frac{\beta}{\gamma_t}x_0||^2)+\beta||x_0||^2 \\ & = \gamma_t||x-\frac{\beta}{\gamma_t}x_0||^2 + \beta(1-\frac{\beta}{\gamma_t})||x_0||^2 \end{align*}\]

where \(\gamma_t = \frac{8D\beta t + 1}{8Dt}\).

We can use this equality to have :

\[ \mathbb{E}[e^{-\beta||X_t-x_0||^2}] = \frac{e^{-\beta(1-\frac{\beta}{\gamma_t})||x_0||^2}}{8\pi Dt}\int_{\mathbb{R}^2} e^{-\gamma_t||x-\frac{\beta}{\gamma_t}x_0||^2 }\d{x} \]

But, it is well known that : \[ \int_{\mathbb{R}^2} e^{-\gamma_t||x-\frac{\beta}{\gamma_t}x_0||^2 }\d{x} = \frac{\pi}{\gamma_t} \]

Thus we finally have:

\[ h(t) = \mathbb{E}[\alpha e^{-\beta||X_t-x_0||^2}] =\frac{\alpha}{1+8D\beta t}e^{-\frac{\beta||x_0||^2}{1+8D\beta t}} \]

These are the corresponding curve for different location of the trap :

Figure 4.1: Hazard function of the capture time for different location of the trap.

Integral modification for survival function approximation

We know that :

\[ S(t) = \exp(-\int_0^th(s)ds) \]

but this integral is not calculable for our hazard function.

Thus, we have to approximate this value. However, the fact that the bound of the interval goes to infinity makes the computation really difficult. But, we can rewrite the integral (by using the transformation \(u = \frac{1}{1+8\beta D s}\)) :

\[\begin{equation} \begin{aligned} \int_0^th(s)ds & = \int_0^t \frac{\alpha}{1+8D\beta s}e^{-\frac{\beta||x_0||^2}{1+8D\beta s}} ds \\ & = \int_1^{\frac{1}{1+ 8D \beta t}} \alpha u e^{-\beta||x_0||^2u} (-\frac{1}{8D \beta u^2})du \\ & = \frac{\alpha}{8D \beta} \int_{\frac{1}{1+ 8D \beta t}}^1 \frac{1}{u}e^{-\beta||x_0||^2u}du \end{aligned} \end{equation}\]

which will easier to compute for each terms (see function cumhazard_one in the file functions.R.

4.2 With several traps

Now, we want to compute this in the case that there are \(I\) trap located at position \(x_k\) during the experiment.

FM: What is “this”? where does this come from?

\[ h_k(t |x)= \frac{\omega_k(x)}{W(x)}\lambda_k(x) \]

Thus, the previous reasoning from the one trap case gives :

\[ h_k(t )= \int_{\mathbb{R}^2}\frac{\omega_k(x)}{W(x)}\lambda_k(x)f_{X_t}(x)\d{x} \]

We can rewrite this formula :

\[\begin{align*} h_k(t ) &= \int_{\mathbb{R}^2}\frac{\omega_k(x)}{W(x)}\lambda_k(x)f_{X_t}(x)\d{x} \\ & = \alpha\int_{\mathbb{R}^2}\frac{\omega_k(x)}{W(x)}e^{-\beta||x-x_k||^2}f_{X_t}(x)\d{x} \\ & = \frac{\alpha}{8 \pi Dt}\int_{\mathbb{R}^2}\frac{\omega_k(x)}{W(x)}e^{-(\frac{||x||^2}{8Dt}+\beta ||x-x_k||^2)}\d{x} \end{align*}\]

Using what we have done in the previous section, we can deduce that :

\[ h_k(t) = \frac{\alpha e^{-\frac{\beta||x_k||^2}{1+8D\beta t}}}{8 \pi Dt}\int_{\mathbb{R}^2}\frac{\omega_k(x)}{W(x)}e^{-\gamma_t||x-\frac{\beta}{\gamma_t}x_k||^2}\d{x} = \frac{\alpha e^{-\frac{\beta||x_k||^2}{1+8D\beta t}}}{1+8D\beta t}\int_{\mathbb{R}^2}\frac{\omega_k(y)}{W(y)}f_{Y_t}(y)dy \]

where \(Y_t \sim \mathcal{N}(\frac{\beta}{\gamma_t}x_k, \frac{\gamma_t}{2}I_2)\). Thus we can compute the hazard function seeing that :

\[ h_k(t)=\frac{\alpha e^{-\frac{\beta||x_k||^2}{1+8D\beta t}}}{1+8D\beta t} \mathbb{E}\Big[\frac{\omega_k(Y_t)}{W(Y_t)}\Big] \]

It’s funny because we can see \(Y_t\) as a path to the trap \(k\).

7.1 Type of censor

There is a lot of different type of censor. In our context, the capture times are censored randomly and at the end, we will get an interval-censored data.

Formally, the model of this type of data is defined as a triplet of random variable \((L, T, R)\) where \(R\) is the value of the right censored bound, \(L\) is the left censored bound and \(T\) the real capture time.

We also get another random variable, that we call \(\delta\) which can take three values :

• \(\delta =1\) if we know the true value of \(T\)

• \(\delta = 2\) if \(R = \infty\)

• \(\delta = 3\) else

The last case is the casual interval-censored case.

We can compute an estimator of the survival function with this type of data using the Turnbull’s estimator.

7.2 A different likelihood

With these paradigm, we have to change the structure of the likelihood and take in consideration this variation of the data. Consider a sample of censored data \((t_i,[l_i, r_i], \delta_i)_{i=1,...,n}\), we write :

\[ \mathcal L = \prod_{i: \delta_i=1}f(t_i) \prod_{i: \delta_i=2} S(l_i) \prod_{i: \delta_i=3}(S(l_i)-S(r_i)) \]

But, in our context, we will have \(\delta_i=3, \forall i\), so the likelihood can be written as :

\[ \mathcal L = \prod_{i=1}^n(S(l_i)-S(r_i)) \]

Now, we suppose that we have a sequence of time \((t_0=0,t_1,\dots, t_n)\) which correspond at the time when the trap are checked. For every check at \(t_j\), we get a list of numbers \((n_{1j},\dots, n_{Kj})\) which are the number of capture between \(t_{j_1}\) and \(t_j\) for every trap.

With this setting, we can write our likelihood as :

\[ \mathcal L = \prod_{j=1}^n \prod_{k=1}^K(S_k(t_{j})-S_k(t_{j-1}))^{n_{jk}} \]

So we get the following log-likelihood :

\[ l(N ;D, \beta, \alpha, \lambda) =\sum_{j=1}^n \sum_{k=1}^Kn_{jk} \log \Big(S_k(t_{j};D, \beta, \alpha, \lambda)-S_k(t_{j-1};D, \beta, \alpha, \lambda)\Big) \] where \(N = (n_{jk})\) is the matrix of capture count.

8.1 Maximum likelihood estimator

8.2 Bayesian estimator by MCMC

TODO : - end this section with survival times of the individuals - plot the histogram with the corresponding pdf (for one and several traps cases) - write the cumulative likelihood with specific data - introduce two methods : likelihood maximization, and Bayesian inference with priors defined by \(\alpha \sim \mathcal Beta(a,1)\) \(\beta \sim \mathcal Beta(1,b)\) + research a good candidate for D’s prior and \((a,b)\) parameters. - estimation of the parameters with these two methods and see results with non parametric estimator of the survival function. - application with real data.