A continuous-time adaptive particle filter for estimations under measurement time uncertainties with an application to a plasma-leucine mixed effects model

Abstract

Background

When mathematical modelling is applied to many different application areas, a common task is the estimation of states and parameters based on measurements. With this kind of inference making, uncertainties in the time when the measurements have been taken are often neglected, but especially in applications taken from the life sciences, this kind of errors can considerably influence the estimation results. As an example in the context of personalized medicine, the model-based assessment of the effectiveness of drugs is becoming to play an important role. Systems biology may help here by providing good pharmacokinetic and pharmacodynamic (PK/PD) models. Inference on these systems based on data gained from clinical studies with several patient groups becomes a major challenge. Particle filters are a promising approach to tackle these difficulties but are by itself not ready to handle uncertainties in measurement times.

Results

In this article, we describe a variant of the standard particle filter (PF) algorithm which allows state and parameter estimation with the inclusion of measurement time uncertainties (MTU). The modified particle filter, which we call MTU-PF, also allows the application of an adaptive stepsize choice in the time-continuous case to avoid degeneracy problems. The modification is based on the model assumption of uncertain measurement times. While the assumption of randomness in the measurements themselves is common, the corresponding measurement times are generally taken as deterministic and exactly known. Especially in cases where the data are gained from measurements on blood or tissue samples, a relatively high uncertainty in the true measurement time seems to be a natural assumption. Our method is appropriate in cases where relatively few data are used from a relatively large number of groups or individuals, which introduce mixed effects in the model. This is a typical setting of clinical studies. We demonstrate the method on a small artificial example and apply it to a mixed effects model of plasma-leucine kinetics with data from a clinical study which included 34 patients.

Conclusions

Comparisons of our MTU-PF with the standard PF and with an alternative Maximum Likelihood estimation method on the small artificial example clearly show that the MTU-PF obtains better estimations. Considering the application to the data from the clinical study, the MTU-PF shows a similar performance with respect to the quality of estimated parameters compared with the standard particle filter, but besides that, the MTU algorithm shows to be less prone to degeneration than the standard particle filter.

Background

Measurement time uncertainties

Uncertainty in the time at which a measurement is taken is an often neglected source of random error. While in many application areas, this kind of error is generally small and indeed neglectable (due to automated measurements and precise timings), in others it may be of real influence, especially in the life sciences. As a prominent example, one may consider pharmacokinetic and pharmacodynamic (PK/PD) models which are used to describe the metabolic interactions and the effects of a chemical agent (like a drug or a labelled substance) over time inside an organism, respectively.

A typical population experiment in the PK/PD context consists in the analysis of the contents of the blood plasma of several individuals with respect to concentrations of certain molecules of interest. For this purpose, blood probes have to be taken from each individual at certain (fixed) time points after a certain event has occurred (e.g. a drug or a labelled substance has been applied). It is clear from the setting of the experiments that there is some variation in the real point in time when the blood probe has been taken: the true time when the measurement value has been obtained might be shortly before or after the intended time, and this true measurement time is not known to us. Since the inclusion of those time uncertainties in the model usually makes the analysis more difficult, it is standard to lump the time uncertainties with the measurement error. But especially at early times when concentrations change quickly, this may easily lead to wrong estimations, even if one assumes very high variances of the measurement error (we will demonstrate this later on a simple example). On the other hand, the inclusion of measurement time uncertainties (MTU) in algorithms aiming at inference making in complex models is not straightforward. In this article, we will present a modification of the Particle Filter (PF) algorithm (which we call MTU-PF) which is able to fully include a statistical model of the time uncertainties.

Inference in complex systems

The assessment of the effectiveness of a drug in a clinical study has been done in the past by the direct computation of relatively simple statistical values. The enormous increase in complexity of the underlying models, due to present developments in medicine and biology, for instance in the areas of personalized medicine or systems biology, increases also the need for more sophisticated model-based inference methods.

The estimation of unobservable internal variables or model parameters from data which have been obtained from blood or tissue samples at several time points can reveal information on the concentrations and effectiveness of the substance under question. If these data come from individuals which belong to two different (or even more) groups, e.g. test and control group, mixed effects are introduced in the underlying models. The inherent non-linearity and high variability of biological processes adds considerably to the difficulties one faces during the inference step. Inference in connection with dynamic models plays a major role in many other application areas. State and parameter estimation as well as model discrimination and validation are most common, but also optimal control problems should be mentioned.

It is often not enough to consider (independent) measurement noise [1]. Correlations between residuals are not uncommon, and the violation of this statistical assumption may lead to wrong estimates. A natural way to include correlated noise is to model two different types of noise: the dynamic (process or system) noise which is present in the dynamics of the system states and originates either from true random fluctuations in the system or from unmodelled dynamics in the system, and the measurement noise which is introduced by the measurement procedure or equipment and modelled by independent residuals. One possible approach is to use state space models which consist of a time-continuous model for the system states, e.g. based on Stochastic Differential Equations (SDEs), and a separate model for the time-discrete measurements.

Parameter estimation with Maximum Likelihood approach

Parameter estimation in state space systems is a difficult problem. In a context where the system dynamics are modelled by Ordinary Differential Equations (ODEs) without correlated noise, the problem is most often considered as a (deterministic) optimization problem based on a Maximum Likelihood (ML) formulation. An overview of these approaches can be found in [2] and [3]; see also [4], which consider also other aspects like identifiability. A generalization of the ML approach including more flexible cost functions is given by the prediction error estimation method ([5]). The introduction of system noise in the state variables leads to optimization problems with SDE constraints. In this case, internal system states which cannot be directly observed need to be estimated jointly with the parameters, given the data. For this purpose, the parameter estimation methods must be augmented by appropriate state filtering methods. An overview of ML parameter estimation in these types of models is given in [6]. If the SDEs are non-linear, linearizations to the Kalman Filter, like the Extended Kalman Filter (EKF) or the Unscented Kalman Filter (UKF), are used to establish approximations to the means and covariances of the filter distributions over time. All those approximations suffer from the fact that they approximate the filtering distributions of the states by a Gaussian distribution at all time points and cannot adequately approximate skewed or multimodal distributions. Better approximations are provided by simulation based methods like Sequential Monte Carlo (SMC) algorithms where good convergence results have been established ([7]). Nevertheless, they suffer from several drawbacks when applied to the joint estimation of dynamic states and fixed parameters ([8–10], see also [11]).

Parameter estimation in a Bayesian context

In a Bayesian context, in contrast to the “classical” ML approach, a prior distribution is assigned to the parameter vector, hence the parameters can be treated as random variables. In this sense, parameter estimation is done by evaluating the so-called posterior distribution which can be computed (at least theoretically) by Bayes’ theorem given the observations (measurements) and the prior distribution. In the context of high-dimensional spaces, this requires the computation of high-dimensional integrals which is not possible to do analytically. For this purpose, Markov Chain Monte Carlo (MCMC) methods provide powerful tools for the computation of simulation-based approximations to the posterior distribution. Again, in the context of the joint estimation of dynamic states and fixed parameters, the design of good proposal densities is a very difficult problem which renders the use of standard MCMC methods like the Metropolis-Hastings sampler impractical for the purposes of parameter estimation in state space systems.

It has long been a wish to combine both (dynamic) SMC and (static) MCMC methods to provide a general tool for the joint estimation of dynamic states and static parameters. Only recently, Andrieu et al. [11] proposed a very promising combination of both types of Monte Carlo approaches called Particle Markov Chain Monte Carlo (PMCMC) which is generally applicable and where also convergence has been proved.

In the present article, even though the PMCMC approach might be the preferred method for parameter estimation in state space systems, we will concentrate solely on the SMC methods, since our modification affects only this part. However, to be able to do parameter estimation in a pure SMC context, we rely on an approach that is very often used to avoid problems with the estimation of constant parameters. This approach consists in the introduction of artificial dynamics in the parameters, that means the parameters are allowed to slightly change their values over time. In this way, and in a Bayesian context, the parameters can be treated exactly in the same way as the system states. After building an augmented system state by concatenating the parameter vector and the state vector, the joint estimation of states and parameters reduces to filtering of the augmented state vector which makes SMC methods directly applicable to the problem.

Particle filters for state and parameter estimation

Particle filters ([12–14]) belong to the class of SMC methods for state filtering in state space models. Using the state augmentation approach, the method is also capable of estimating system parameters. The standard particle filter is designed for discrete, non-linear, and non-Gaussian models and can routinely be adapted to the continuous case with measurements at discrete times. The idea of the particle filter is that, at each time point, there is a sample based representation (the weighted particles) of the current estimate of the inner states and parameters which is based on the measurements that have been obtained up to the current time point. The particle cloud is then propagated through time, and the particles and weights are updated accordingly at each time point where measurements are available.

Non-Linear Mixed Effects models

Estimation in a Non-linear Mixed Effects model (NLME) involves the estimation of both global and individual parameters. With classical maximum likelihood estimation, the individual parameters are random variables equipped with a distribution while the global parameters remain constants with a “true” but unknown value. If the underlying model equations are non-linear, this leads to likelihood functions which are not analytically accessible and one has to rely on approximations. In the context where the system dynamics are modelled by ODEs, the most popular algorithm for NLME parameter estimation in the PK/PD context is the tool NONMEM ([15]). In [1] an estimation algorithm for NLME models based on Stochastic Differential Equations (SDEs) was proposed that uses the First-Order Conditional Estimation (FOCE) method to approximate the likelihood in combination with the EKF estimation in the SDEs. This has been added to NONMEM ([16]). In [17], a comparison between ODE and SDE based parameter estimation has been performed which showed that the interindividual variabilities were in general estimated to be smaller for the SDE model. Donnet and Samson ([18]) proposed a stochastic version of the Expectation-Maximization (SAEM) algorithm (for the estimation of the global parameters) in combination with MCMC methods (for the estimation of states and individual parameters). However, since MCMC exhibits slow mixing properties in the context of the estimation of states and parameters in state space models, in [19] MCMC has been replaced by the more promising PMCMC approach of Andrieu et al. ([11]).

On the other hand, in a Bayesian context, also the global parameters are equipped with a (prior) probability distribution, and the conceptual difference between global and individual parameters vanishes. The mixed effects model can then be considered as a hierarchical model with dependent parameters ([20, 21], see also [22] for a more recent population-based Bayesian approach to PK/PD modelling). Simulation-based (Monte Carlo) methods can easily be adapted to this case. Nevertheless, the above mentioned challenges to both SMC and MCMC methods are even higher due to the increased number of states and parameters in NLME models (the number of states and individual parameters has to be multiplied by the number of individuals).

Aim of the article

Our goal is two-fold: Firstly, we want to show that the particle filter algorithm is applicable (with our modifications) also to more complex models when time uncertainties are formulated explicitly. Secondly, we want to show that the modification may even provide the possibility for further enhancement of the performance of the algorithm by presenting an adaptive time-stepping scheme which is only possible in the context of the new algorithm.

We do not claim that our MTU algorithm generally performs better or worse than the standard filter, nor that it should be the preferred method for estimation in non-linear mixed effects models. Rather, we provide a method which is usable for models where time uncertainties may play a major role. In these cases, it may indeed lead to better estimations. On the other hand, our method transfers the time-discrete particle filter approach, where updates based on the measurements very strictly depend on the measurement times, to a truly time-continuous approach, where updates to the filtering distributions can be performed at every point on the time-scale. Since we want to focus on the time uncertainties, we neglect discussing further issues like identifiability, model evaluation and model discrimination. In our application to the model of plasma-leucine kinetics, we try to avoid these issues by providing ad-hoc values to some of the parameters (especially to the variances of the system states).

Motivating example

Let us have a look at an example for illustrating the benefits of a separate modelling of measurement time uncertainties. Let us consider a state space system given by the ODE

$$\mathrm{d}q\left(t\right)=(-\mathrm{\alpha q}(t)+\beta )\mathrm{d}t$$

with parameters $\alpha ,\beta \in \mathbb{R}$. Here $q\left(t\right)\in \mathbb{R}$ denotes the state of the system at time t. We call the state trajectories obtained by this deterministic system the nominal evolutions of the states. We add noise to the system in a standard way by introducing an additional term $\sigma \mathrm{d}{\mathcal{W}}_t$, with a standard Wiener process ${\left({\mathcal{W}}_t\right)}_{t\in {\mathbb{R}}^{\ge 0}}$ and a diffusion parameter $\sigma \in \mathbb{R}$. This leads to the following SDE describing the evolution of the state q over time:

$$\mathrm{d}q\left(t\right)=(-\alpha q(t)+\beta )\mathrm{d}t+\sigma \mathrm{d}{\mathcal{W}}_t.$$

Furthermore, let the initial state q(0) be given by a log-normal distribution with parameters ${\mu }_{q_0}$ and ${\sigma }_{q_0}^2$ (mean and variance of the logarithm of q(0), respectively). The parameters chosen in our implementation of this example are shown in Table 1.

Table 1 Parameters for the motivating example

We assume that M measurements of the state q(t) will be taken at times T _j and that each measurement j, j=1,…,M is disturbed by normal noise with mean q(t _j ) and with a fixed variance ${\sigma }_y^2$, i.e. measurement j is distributed according to $y_j\sim \mathcal{N}\left(q\right(t_j),{\sigma }_y^2)$. Usually, the times T _j are assumed to be known. In contrast, we will assume that in addition to the measurement value error, there is some uncertainty about the exact times where the measurements have been taken. If we attempt to take the j th measurement at the intended (or nominal) time ${\widehat{t}}_j$, the measurement in fact takes place at time T _j which may be shortly before or after the intended time ${\widehat{t}}_j$. A natural way to model these uncertainties is to assume that the measurement time T _j is given as a realization of a random variable T _j . In our example, we assume that T _j follows a truncated normal distribution given by the density

$${\gamma }_j\left(t_j\right):=\frac{1}{{\Gamma }_j}\left\{\begin{array}{c}exp\left(\frac{{(t_j-{\widehat{t}}_j)}^2}{0.3^2}\right)\text{if}max\{0,{\widehat{t}}_j-1\}\\ \le t_j\le {\widehat{t}}_j+1,\\ 0\text{otherwise}\end{array}\right.$$

with normalizing constant

$${\Gamma }_j:={\int }_{max\{0,{\widehat{t}}_{j-1}\}}^{{\widehat{t}}_{j+1}}exp\left(\frac{{(t_j-{\widehat{t}}_j)}^2}{0.3^2}\right)$$

and given intended measurement times ${\widehat{t}}_j$. Figure 1 shows the different distributions for one measurement j for all possible intended times ${\widehat{t}}_j$. In each of the subfigures (a)-(d), the shaded green area gives an impression of the “density” of the distribution of the measurement value y _j in dependence of the intended measurement time ${\widehat{t}}_j$ on the x-axis, while the dark-green dashed line depicts the nominal evolution of the state q over time. Subfigure (a) shows the distribution of the measurement values with time uncertainties, while (b)-(d) depict the distribution of the measurement values with known measurement times (in this case ${\widehat{t}}_j=t_j$), for several standard deviations σ _y : in (b), the original standard deviation is used, while in (c) and (d), higher standard deviations are used which correspond to the cases with lumped value and time variations.

Figure 1

Assumed measurement distributions for the motivating example. Measurement distribution resulting from (a) separate modelling of measurement time uncertainties and measurement value uncertainties, with σ _y =0.005, and (b-d) lumped time-and-value uncertainties with several different assumed lumped measurement variances σ _y . The dashed dark-green line depicts the nominal evolution of the state q over time. The green shaded area depicts the region where the measurements are expected.

Comparing Figures 1(a) and 1(b-d), we observe that the distributions of the measurements exhibit clearly different shapes. For the “true” model depicted in Figure 1(a), if we consider a single point in time that lies in a time segment where the state values change quickly, the distribution of the measurement at this certain point in time is quite broad. The variance in the measured value is very high, whereas it is small in time segments where the state values change slowly. In contrast, for the standard particle filter, the measurement variance is constant and hence the assumed measurement distributions differ remarkably from the “true” distributions, howsoever the value of ${\sigma }_y^2$ is chosen. It must be expected that this leads to difficulties when inference on the states and parameters needs to be done based on these models. We will resume our example after having presented the MTU particle filter and will show that this is indeed the case.

Methods

We divide this section into three subsections. In the first subsection, we fix the state and observation model we want to consider. In the second subsection entitled “Standard case” we outline the standard particle filter algorithm in the context of time-continuous states with time-discrete measurements, and the various probability distributions involved. Although nothing is new in this subsection, it serves several purposes. Firstly, the time-continuous case is relatively rarely considered in the literature; secondly, the derivation of our modification needs a slightly more general formulation than it is standard for the discrete-time filter; and lastly, the comparison of our modified version with the standard case might more clearly reveal the differences between the two approaches. In the third subsection entitled “MTU particle filter”, we present our new modification of the particle filter. In the following section “Results and Discussion”, we compare the new MTU particle filter to the standard particle filter and to an alternative Maximum Likelihood estimation method on a simple artificial example. We also present an application of our MTU-PF method to a PK/PD study in a non-linear mixed-effects setting in direct comparison with the standard particle filter.

Note: a list of all used symbols with a short explanation can be found at the end of this paper.

The model

State process

Let $(\Omega ,\mathcal{A},\mathsf{\text{P}})$ be a probability space and for each t∈ [t₀,∞) with $t_0\in \mathbb{R}$ let $({\mathcal{X}}_t,{\mathcal{B}}_{{\mathcal{X}}_t})$ be an arbitrary measurable space. For each t∈ [t₀,∞) let further $X_t:\Omega \to {\mathcal{X}}_t$ be an $\mathcal{A}-{\mathcal{B}}_{{\mathcal{X}}_t}$ measurable random variable such that $X_{[t_0,\infty )}:={\left(X_t\right)}_{t\in [t_0,\infty )}$ is a continuous-time Markov process with general state space

$${\mathcal{X}}_{[t_0,\infty )}:=\prod _{t_0\le s}{\mathcal{X}}_s.$$

For each t∈ [t₀,∞), denote by ${\mathcal{L}}_{X_t}$ the pushforward measure of P under X _t , i.e. ${\mathcal{L}}_{X_t}\left(B\right):=\mathsf{\text{P}}\left(X_t^{-1}\right(B\left)\right)$ for all $B\in {\mathcal{B}}_{{\mathcal{X}}_t}$. Further, denote by ${\mathcal{L}}_{X_{[t_0,\infty )}}$ the pushforward measure of P under $X_{[t_0,\infty )}:={\left(X_s\right)}_{s\in [t_0,\infty )}$ (with the corresponding product algebra). Analogously, denote by

$${\mathcal{X}}_{[t_0,t]}:=\prod _{t_0\le s\le t}{\mathcal{X}}_s\text{for each}t\ge t_0$$

the state space restricted to the interval [t₀,t], and denote by ${\mathcal{L}}_{X_{[t_0,t]}}$ the corresponding pushforward measure. For each s and t with t>s≥t₀, let K_s,t(x _s , dx _t ) be the Markov kernel of the process $X_{[t_0,\infty )}$ from time s to time t.

An important special case for $X_{[t_0,\infty )}$ is given by a multidimensional Itô process on ${\mathcal{X}}_t={\mathbb{R}}^n$ (equipped with the corresponding Borel σ-algebra) defined through a stochastic differential equation (SDE)

$$\mathrm{d}{X}_t=a(X_t,t)\mathrm{d}t+B(X_t,t)\mathrm{d}{\mathcal{W}}_t$$

with drift a(x,t), diffusion matrix B(x,t), multidimensional standard Wiener process ${\mathcal{W}}_t$, and initial variable $X_{t_0}$. In this case, it is possible to sample directly (at least approximately) from the kernels K_s,t when a suitable discretization method is applied, for instance the Euler-Maruyama method.

Observations / measurements

Let the process $X_{[t_0,\infty )}$ be observed via M random variables Y_1:M with values in measurable spaces $({\mathcal{Y}}_j,{\mathcal{B}}_{{\mathcal{Y}}_j})$. Each single observation (measurement) y _j depends on the state variable $x_{t_j}$ at some time T _j and on the observation time (measurement time) T _j itself. We assume that, given the observation time T _j and the state $X_{t_j}=x_{t_j}$, the variable y _j is independent of all other variables, and the conditional measure can be expressed via some conditional probability density $g_j\left(y_j\right|x_{t_j},t_j)$ with respect to a reference measure ${\mu }_{{\mathcal{Y}}_j}$ on $({\mathcal{Y}}_j,{\mathcal{B}}_{{\mathcal{Y}}_j})$. We do not require any further restrictions on g such as linear dependence on the states or Gaussianity.

Observation / measurement times

The observation times (measurement times) T _j for j=1,…,M are usually assumed to be deterministically given and known. Our variant of the particle filter will be based on the assumption that the observation times T _j are themselves realizations of random variables T _j . These variables model the uncertainty about exact observation times. In contrast to the observation variables y _j , the observation times T _j are never observed (measured). We assume that all information available to us is their probability distribution on the half axis [t₀,∞), while in the case of the observations y _j , we know both the densities $g_j\left(y_j\right|x_{t_j},t_j)$and the observed values y _j .

In this article, we will only consider the simplest case where each variable T _j is independent of all others. Dependencies between the T _j ’s, especially concerning the order of the observation times, may be considered natural but would lead to more complicated algorithms. However, order dependencies can easily be introduced via restrictions on the support of the variables. In general, the probability distribution of every single variable T _j shall be given by a density γ _j (t _j ) with respect to the Lebesgue measure ${\lambda }_{[t_0,\infty )}$ on the interval [t₀,∞).

In the following, we will consider the two cases mentioned, where either all T _j are deterministic and known or all T _j are random and unknown. Note that the first case formally coincides with the case that T _j is random but observed. We will therefore stick to the notation $g_j\left(y_j\right|x_{t_j},t_j)$ for the observation densities in both cases.

Standard case: measurement times deterministic and known

We will first consider the standard case, where the observation times T _j are known. For simplicity, we assume here that the observation times t_1:M are strictly ordered increasingly, i.e. t₀<t₁⋯<t _M .

The standard case of the particle filter is usually formulated for discrete-time Markov processes $X_{t_{0:M}}:={\left(X_{t_j}\right)}_{j\in \{0,\dots ,M\}}$ with general state space where the state variables are only defined at the initial time t₀ and at the times t₁,…,t _M when measurements occur. Nevertheless, this case is included in our more general framework where X _t is defined for all t≥t₀. One just focuses on the state variables for those times only. In view of the later generalization to random observation times, we will consider the fixed values T _j as realizations of random variables T _j and condition all occurring densities on them. As mentioned above this assumption leads to the same results as if we assumed the values T _j to be given deterministically.

Full model and filter model

The full model is given by the joint density of the variables $X_{t_{0:M}}$ and Y_1:M (conditioned on the observation times T_1:M=t_1:M) with respect to the product measure ${\mathcal{L}}_{X_{t_{0:M}}}\prod _{j=1}^M{\mu }_{{\mathcal{Y}}_j}$:

$$f^{X_{t_{0:M}},Y_{1:M}|T_{1:M}}(x_{t_{0:M}},y_{1:M}|t_{1:M}):=\prod _{j=1}^M{g}_j(y_j|x_{t_j},t_j).$$

(1)

The filter at a given time t _k is based on a reduced model. This model is given by the joint density of the variables $X_{t_{0:k}}$ and Y_1:k (conditioned on T_1:M=t_1:M) with respect to the product measure ${\mathcal{L}}_{X_{t_{0:k}}}\prod _{j=1}^k{\mu }_{{\mathcal{Y}}_j}$:

$$f^{X_{t_{0:k}},Y_{1:k}|T_{1:M}}(x_{t_{0:k}},y_{1:k}|t_{1:M}):=\prod _{j=1}^k{g}_j(y_j|x_{t_j},t_j).$$

(2)

This density is based on the state sequence $X_{t_{0:k}}$. In contrast, we can focus on the single state $X_{t_k}$ by considering the joint density of the variables $X_{t_k}$ and Y_1:k (given T_1:M=t_1:M) with respect to ${\mathcal{L}}_{X_{t_k}}\prod _{j=1}^k{\mu }_{{\mathcal{Y}}_j}$. It can be computed by marginalization as follows:

$$\begin{array}{l}{f}^{X_{t_k},Y_{1:k}|T_{1:M}}(x_{t_k},y_{1:k}|t_{1:M})\\ :={\int }_{\left\{{\tilde{x}}_{t_{0:k}}\in {\mathcal{X}}_{t_{0:k}}:{\tilde{x}}_{t_k}=x_{t_k}\right\}}{f}^{X_{t_{0:k}},Y_{1:k}|T_{1:M}}\left({\tilde{x}}_{t_{0:k}},y_{1:k}|t_{1:M}\right)\\ \times \mathrm{d}{\mathcal{L}}_{X_{t_{0:k}}}\left({\tilde{x}}_{t_{0:k}}\right)\end{array}$$

(3)

and the filter density at time t _k with respect to ${\mathcal{L}}_{X_{t_k}}$ can then be computed with Bayes’ theorem:

$$f^{X_{t_k}|Y_{1:k},T_{1:M}}\left(x_{t_k}\right|y_{1:k},t_{1:M}):=\frac{f^{X_{t_k},Y_{1:k}|T_{1:M}}(x_{t_k},y_{1:k}|t_{1:M})}{f^{Y_{1:k}|T_{1:M}}\left(y_{1:k}\right|t_{1:M})}$$

(4)

with

$$\begin{array}{c}{f}^{Y_{1:k}|T_{1:M}}\left(y_{1:k}\right|t_{1:M})\hfill \\ :={\int }_{{\mathcal{X}}_{t_{0:k}}}{f}^{X_{t_{0:k}},Y_{1:k}|T_{1:M}}\left(x_{t_{0:k}},y_{1:k}|t_{1:M}\right)\mathrm{d}{\mathcal{L}}_{X_{t_{0:k}}}\left(x_{t_{0:k}}\right).\hfill \end{array}$$

(5)

For general (non-linear) models, the practical computation of the filter density is very difficult. Nevertheless, the particle filter computes a Monte Carlo approximation using the fact that the filter densities $f^{X_{t_k}|Y_{1:k},T_{1:M}}$ can be computed recursively. This is done in two steps. We consider the filter distribution at time t_k-1 given by the probabilities

$$\begin{array}{c}\mathsf{\text{P}}\left(X_{t_{k-1}}\in B|Y_{1:k-1}=y_{1:k-1},T_{1:M}=t_{1:M}\right)\hfill \\ ={\int }_B{f}^{X_{t_{k-1}}|Y_{1:k-1},T_{1:M}}\left(x_{t_{k-1}}|y_{1:k-1},t_{1:M}\right)\mathrm{d}{\mathcal{L}}_{X_{t_{k-1}}}\left(x_{t_{k-1}}\right)\hfill \end{array}$$

(6)

for each set $B\in {\mathcal{B}}_{{\mathcal{X}}_{t_{k-1}}}$. We then get first the prediction distribution, i.e. the distribution of $X_{t_k}$ given the data Y_1:k-1 (and t_1:M), by use of the kernel $K_{t_{k-1},t_k}$:

$$\begin{array}{c}\mathsf{\text{P}}(X_{t_k}\in B|Y_{1:k-1}=y_{1:k-1},T_{1:M}=t_{1:M})\hfill \\ ={\int }_B{\int }_{{\mathcal{X}}_{t_{k-1}}}{f}^{X_{t_{k-1}}|Y_{1:k-1},T_{1:M}}\left(x_{t_{k-1}}|y_{1:k-1},t_{1:M}\right)\hfill \\ \times \mathrm{d}{\mathcal{L}}_{X_{t_{k-1}}}\left(x_{t_{k-1}}\right)K_{t_{k-1},t_k}\left(x_{t_{k-1}},\mathrm{d}{x}_{t_k}\right)\end{array}$$

(7)

for each set $B\in {\mathcal{B}}_{{\mathcal{X}}_{t_k}}$. Then we use Bayes’ theorem to get the filter distribution at time t _k :

$$\begin{array}{c}\mathsf{\text{P}}(X_{t_k}\in B|Y_{1:k}=y_{1:k},T_{1:M}=t_{1:M})\hfill \\ ={\int }_B\frac{g_k\left(y_k\right|x_{t_k},t_k)}{f^{Y_k|Y_{1:k-1},T_{1:M}}\left(y_k\right|y_{1:k-1},t_{1:M})}·\hfill \\ {\int }_{{\mathcal{X}}_{t_{k-1}}}{f}^{X_{t_{k-1}}|Y_{1:k-1},T_{1:M}}\left(x_{t_{k-1}}\right|y_{1:k-1},t_{1:M})\hfill \\ \times \mathrm{d}{\mathcal{L}}_{X_{t_{k-1}}}\left(x_{t_{k-1}}\right)K_{t_{k-1},t_k}\left(x_{t_{k-1}},\mathrm{d}{x}_{t_k}\right)\hfill \end{array}$$

(8)

for each set $B\in {\mathcal{B}}_{{\mathcal{X}}_{t_k}}$, with normalizing constant

$$\begin{array}{c}{f}^{Y_k|Y_{1:k-1},T_{1:M}}\left(y_k\right|y_{1:k-1},t_{1:M})\hfill \\ :={\int }_{{\mathcal{X}}_{t_k}}{g}_k\left(y_k\right|x_{t_k},t_k)·\hfill \\ {\int }_{{\mathcal{X}}_{t_{k-1}}}{f}^{X_{t_{k-1}}|Y_{1:k-1},T_{1:M}}\left(x_{t_{k-1}}|y_{1:k-1},t_{1:M}\right)\hfill \\ \times \mathrm{d}{\mathcal{L}}_{X_{t_{k-1}}}\left(x_{t_{k-1}}\right)K_{t_{k-1},t_k}\left(x_{t_{k-1}},\mathrm{d}{x}_{t_k}\right).\hfill \end{array}$$

(9)

Importance sampling

Another ingredient for the particle filter is sequential importance sampling. We assume that a second Markov chain ${\tilde{X}}_{t_{0:M}}$ on the same state space is given with pushforward measures ${\mathcal{L}}_{{\tilde{X}}_{t_j}}$ and kernels ${\tilde{K}}_{t_{j-1},t_j}\left(x_{t_{j-1}},\mathrm{d}{x}_{t_j}\right)$ for j=1,…,M. We assume that for each $x_{t_{j-1}}\in {\mathcal{X}}_{t_{j-1}}$, the measure $K_{t_{j-1},t_j}\left(x_{t_{j-1}},·\right)$ is absolutely continuous with respect to the measure ${\tilde{K}}_{t_{j-1},t_j}\left(x_{t_{j-1}},·\right)$. It follows that the Radon-Nikodym derivative (written as conditional density)

$${\varrho }_{t_j|t_{j-1}}\left(x_{t_j}\right|x_{t_{j-1}}):=\frac{K_{t_{j-1},t_j}(x_{t_{j-1}},\mathrm{d}{x}_{t_j})}{{\tilde{K}}_{t_{j-1},t_j}(x_{t_{j-1}},\mathrm{d}{x}_{t_j})}$$

exists. We further assume that the pushforward measure ${\mathcal{L}}_{X_{t_0}}$ under $X_{t_0}$ is absolutely continuous with respect to the corresponding pushforward measure ${\mathcal{L}}_{{\tilde{X}}_{t_0}}$ under ${\tilde{X}}_{t_0}$ with Radon-Nikodym derivative

$${\varrho }_{t_0}\left(x_{t_0}\right):=\frac{\mathrm{d}{\mathcal{L}}_{X_{t_0}}\left(x_{t_0}\right)}{\mathrm{d}{\mathcal{L}}_{{\tilde{X}}_{t_0}}\left(x_{t_0}\right)}.$$

For sequential importance sampling, we need to be able to sample from the initial measure ${\mathcal{L}}_{{\tilde{X}}_{t_0}}$ and from the kernels

$${\tilde{K}}_{t_{j-1},t_j}\left(x_{t_{j-1}},·\right)$$

for each $x_{t_{j-1}}\in {\mathcal{X}}_{t_{j-1}}$, and to compute ${\varrho }_{t_0}\left(x_{t_0}\right)$ as well as ${\varrho }_{t_j|t_{j-1}}\left(x_{t_j}\right|x_{t_{j-1}})$ pointwise.

Using

$$K_{t_{k-1},t_k}\left(x_{t_{k-1}},\mathrm{d}{x}_{t_k}\right)={\varrho }_{t_k|t_{k-1}}\left(x_{t_k}|x_{t_{k-1}}\right){\tilde{K}}_{t_{k-1},t_k}\left(x_{t_{k-1}},\mathrm{d}{x}_{t_k}\right),$$

we can then write the recursive formula (8) for the filter distribution at time t _k as

$$\begin{array}{c}\mathsf{\text{P}}(X_{t_k}\in B|Y_{1:k}=y_{1:k},T_{1:M}=t_{1:M})\hfill \\ ={\int }_B\frac{g_k\left(y_k\right|x_{t_k},t_k)}{f^{Y_k|Y_{1:k-1},T_{1:M}}\left(y_k\right|y_{1:k-1},t_{1:M})}·\hfill \\ {\int }_{{\mathcal{X}}_{t_{k-1}}}{f}^{X_{t_{k-1}}|Y_{1:k-1},T_{1:M}}\left(x_{t_{k-1}}|y_{1:k-1},t_{1:M}\right){\varrho }_{t_k|t_{k-1}}\left(x_{t_k}|x_{t_{k-1}}\right)\hfill \\ \times \mathrm{d}{\mathcal{L}}_{X_{t_{k-1}}}\left(x_{t_{k-1}}\right){\tilde{K}}_{t_{k-1},t_k}(x_{t_{k-1}},\mathrm{d}{x}_{t_k})\hfill \end{array}$$

(10)

for each $B\in {\mathcal{B}}_{{\mathcal{X}}_{t_k}}$. The direct computation of the normalizing constants $f^{Y_k|Y_{1:k-1},T_{1:M}}\left(y_k\right|y_{1:k-1},t_{1:M})$ (while Y_1:M is assumed to be fixed) is not necessary. Sequential importance sampling is performed as follows. Draw a number N of realizations $x_{t_0}^i$ from ${\mathcal{L}}_{{\tilde{X}}_{t_0}}$ and compute the corresponding unnormalized weights

$$w_{t_0}^i:={\varrho }_{t_0}\left(x_{t_0}^i\right)\text{for all}i=1,\dots ,N.$$

Then, for all k=1,…,M, sample realizations $x_{t_k}^i$ from the kernel ${\tilde{K}}_{t_{k-1},t_k}(x_{t_{k-1}}^i,\mathrm{d}{x}_{t_k})$ for each i=1,…,N and compute the unnormalized weights

$$\begin{array}{l}{w}_{t_k}^i:={\varrho }_{t_k|t_{k-1}}\left(x_{t_k}^i|x_{t_{k-1}}^i\right)g_k\left(y_k|x_{t_k}^i,t_k\right)w_{t_{k-1}}^i\text{for all}\\ i=1,\dots ,N.\hfill \end{array}$$

For suitable integrable functions h (e.g. fulfilling some mild restrictions on how fast h may increase with x, see [23] for details), the expectation of h with respect to the filter density conditioned on the observations Y_1:k=y_1:k, given by

$$\begin{array}{c}\mathsf{\text{E}}\left[h\left(X_{t_k}\right)|Y_{1:k}=y_{1:k},T_{1:M}=t_{1:M}\right]\hfill \\ :={\mathsf{\text{E}}}_{f^{X_{t_k}|Y_{1:k}=y_{1:k},T_{1:M}=t_{1:M}}(·|y_{1:k},t_{1:M})}\left[h\left(X_{t_k}\right)\right]\hfill \\ =\int f^{X_{t_k}|Y_{1:k},T_{1:M}}\left(x_{t_k}\right|y_{1:k},t_{1:M}\left)h\right(x_{t_k})\mathrm{d}{\mathcal{L}}_{X_{t_k}}\left(x_{t_k}\right),\hfill \end{array}$$

(11)

can then be approximated by

$$\begin{array}{c}{\mathsf{\text{E}}}_{t_k,N}\left[h\left(X_{t_k}\right)|Y_{1:k}=y_{1:k},T_{1:M}=t_{1:M}\right]:=\frac{\sum _{i=1}^N{w}_{t_k}^{i}h\left(x_{t_k}^i\right)}{\sum _{i=1}^N{w}_{t_k}^i}\end{array}$$

(12)

where N is the number of particles. In fact, it can be shown that as N approaches infinity, these empirical expectations converge to the filter expectations:

$$\begin{array}{c}\underset{N\to \infty }{lim}{\mathsf{\text{E}}}_{t_k,N}\left[h\left(X_{t_k}\right)|Y_{1:k}=y_{1:k},T_{1:M}=t_{1:M}\right]\\ =\mathsf{\text{E}}\left[h\left(X_{t_k}\right)|Y_{1:k}=y_{1:k},T_{1:M}=t_{1:M}\right].\end{array}$$

(13)

Note that if we can sample from the Markov kernels of $x_{t_j}$, we can choose ${\tilde{X}}_{t_j}=X_{t_j}$ (at least in law), whence ${\varrho }_{t_0}\left(x_{t_0}\right)\equiv 1$ and ${\varrho }_{t_j|t_{j-1}}\left(x_{t_j}\right|x_{t_{j-1}})\equiv 1$. This is a standard choice, but in terms of efficiency of the particle filter algorithm not always the best one. On the other hand, finding a suitable Markov chain ${\tilde{X}}_{t_{0:M}}$ different from $X_{t_{0:M}}$ is not an easy task.

Resampling

If the number N of samples through time is fixed, the samples obtained by sequential importance sampling quickly degenerate since most of the normalized weights decrease rapidly towards 0. The degree of degeneracy is often measured by an estimate of the so-called effective sample size (ESS). This estimate at time t is given by

$$n_{\text{ESS}}:=\frac{1}{\sum _{i=1}^N{\left({\tilde{w}}_t^i\right)}^2}$$

(14)

where

$${\tilde{w}}_t^i:=\frac{w_t^i}{\sum _{\nu =1}^N{w}_t^{\nu }}$$

(15)

are the normalized weights. It obtains its maximal value N if all weights are equal, and it approaches 1 if the variance of the weights and thus the degree of degeneracy increases. To avoid this degeneration of the samples, a resampling step needs to be done when the ESS drops below a threshold N_Threshold (which is usually chosen to be N/2).

Resampling at some time s _ℓ is based on given non-negative (unnormalized) selection weights $v_{s_{\ell }}^i$ for each particle index i: One repeatedly selects particles with probabilities $p_{\ell }^i$ given by the normalized selection weights

$$p_{\ell }^i:=\frac{v_{s_{\ell }}^i}{\sum _{\nu =1}^N{v}_{s_{\ell }}^{\nu }}.$$

(16)

This is multinomial resampling. There exist procedures where each single particle is still selected with probability $p_{\ell }^i$, but with reduced overall variance, for instance stratified resampling or systematic resampling which should be preferred [24, 25]. In any case, resampling defines a selection function ι _ℓ :I→I on the index set I: ={1,…,N}. Resampling is then done by replacing the state samples ${\left(x_{s_{\ell }}^i\right)}_{i=1,\dots ,N}$ by the selected state samples ${\left(x_{s_{\ell }}^{{\iota }_{\ell }\left(i\right)}\right)}_{i=1,\dots ,N}$. Since before selection the probability that the particle i will be chosen is $p_{\ell }^i$ for each draw, the expected number of times that particle i has been chosen after N draws is $N{p}_{\ell }^{{\iota }_{\ell }\left(i\right)}$. To correct for the introduced bias, the normalized weight ${\tilde{w}}_{s_{\ell }}^i$ for each selected particle i needs then to be corrected by replacing it by the weight

$$\frac{{\tilde{w}}_{s_{\ell }}^{{\iota }_{\ell }\left(i\right)}}{N{p}_{\ell }^{{\iota }_{\ell }\left(i\right)}}/\sum _{\nu =1}^N\frac{{\tilde{w}}_{s_{\ell }}^{{\iota }_{\ell }\left(\nu \right)}}{N{p}_{\ell }^{{\iota }_{\ell }\left(\nu \right)}}=\frac{w_{s_{\ell }}^{{\iota }_{\ell }\left(i\right)}}{v_{s_{\ell }}^{{\iota }_{\ell }\left(i\right)}}/\sum _{\nu =1}^N\frac{w_{s_{\ell }}^{{\iota }_{\ell }\left(\nu \right)}}{v_{s_{\ell }}^{{\iota }_{\ell }\left(\nu \right)}}$$

(17)

(using (16)). The necessary correction is therefore achieved if the unnormalized weights ${\left(w_{s_{\ell }}^i\right)}_{i=1,\dots ,N}$ are replaced by the corrected unnormalized weights ${(w_{s_{\ell }}^{{\iota }_{\ell }\left(i\right)}/v_{s_{\ell }}^{{\iota }_{\ell }\left(i\right)})}_{i=1,\dots ,N}$.

Note that in the original particle filter, the selection weights $v_{s_{\ell }}^i$ at time s _ℓ are chosen to be the particle weights (before the replacement), i.e.

$$v_{s_{\ell }}^i=w_{s_{\ell }}^i\text{for}i=1,\dots ,N,$$

such that after the resampling step the unnormalized weights are all equal to 1. Nevertheless, in general their choice is free and may be based on the observations (which is used in the so-called auxiliary particle filter [26]).

Particle filter algorithm

The particle filter computes the state realizations and weights recursively through time. In its standard form, the particle filter can be stated as in algorithm 1.

Note that if one chooses ${\tilde{X}}_{[t_0,\infty )}=X_{[t_0,\infty )}$ (in law), then ${\varrho }_{t_k|t_{k-1}}\left(x_{t_k}^i\right|x_{t_{k-1}}^i)\equiv 1$ and the update of the weights simplifies to

$$w_{t_k}^i=g_k\left(y_k|x_{t_k}^i,t_k\right)w_{t_{k-1}}^i.$$

Algorithm 1 Standard particle filter

Data Likelihood

Model validation or discrimination is generally based on the data likelihood

$$\begin{array}{c}{Z}_{t_k}\left(t_{1:M}\right):=f^{Y_{1:k}|T_{1:M}}\left(y_{1:k}\right|t_{1:M})\hfill \\ =\underset{{\mathcal{X}}_{t_{0:k}}}{\int }{f}^{X_{t_{0:k}},Y_{1:k}|T_{1:M}}\left(x_{t_{0:k}},y_{1:k}|t_{1:M}\right)\mathrm{d}{\mathcal{L}}_{X_{t_{0:k}}}\left(x_{t_{0:k}}\right)\hfill \\ =\mathsf{\text{E}}\left[f^{X_{t_{0:k}},Y_{1:k}|T_{1:M}}(·,y_{1:k}|t_{1:M})\right]\hfill \end{array}$$

(18)

for given observations Y_1:k. Without resampling, the data likelihood could be approximated by the empirical mean of the unnormalized weights, i.e. by

$${\widehat{Z}}_{t_k}\left(t_{1:M}\right):=\frac{1}{N}\sum _{i=1}^N{w}_{t_k}^i$$

(19)

because this is the empirical estimate for the above expectation. After a resampling step, this is not valid any longer. Nevertheless, in any case, the data likelihood can be computed recursively by the following estimate of the ratio $Z_{t_k}\left(t_{1:M}\right)/Z_{t_{k-1}}\left(t_{1:M}\right)$:

$$\begin{array}{c}\frac{\widehat{Z_{t_k}\left(t_{1:M}\right)}}{Z_{t_{k-1}}\left(t_{1:M}\right)}:=\frac{\sum _{i=1}^N{\varrho }_{t_k|t_{k-1}}\left(x_{t_k}^i|x_{t_{k-1}}^i\right)g_k\left(y_k|x_{t_k}^i,t_k\right)w_{t_{k-1}}^i}{\sum _{i=1}^N{w}_{t_{k-1}}^i},\end{array}$$

(20)

with initial estimate ${\widehat{Z}}_{t_0}\left(t_{1:M}\right)=1$ (see e.g. [27]).

MTU particle filter: Uncertain measurement times

We now assume that each observation time T _j is a realization of a random variable T _j . Its distribution is expressed via densities γ _j with respect to the Lebesgue measure ${\lambda }_{[t_0,\infty )}$. The observation times T _j themselves are not observed.

Full model

The full model in this case will include complete continuous state paths, since the observation times are now distributed over the complete time axis [t₀,∞), and the observations may potentially depend on every state $x_{t_j}$ for t _j ∈[t₀,∞). Consider therefore the joint density of the variables $X_{[t_0,\infty )}$, Y_1:M and t_1:M, with respect to the product measure ${\mathcal{L}}_{X_{[t_0,\infty )}}\prod _{j=1}^M{\mu }_{{\mathcal{Y}}_j}\prod _{j=1}^M{\lambda }_{[t_0,\infty )}$:

$$\begin{array}{c}{f}^{X_{[t_0,\infty )},Y_{1:M},T_{1:M}}\left(x_{[t_0,\infty )},y_{1:M},t_{1:M}\right)=\prod _{j=1}^M{g}_j\left(y_j\right|x_{t_j},t_j\left)\prod _{j=1}^M{\gamma }_j\right(t_j).\end{array}$$

(21)

Filter model

The filter at a given time t≥t₀ is again based on a reduced model. This model is given by the joint density of the following variables: $x_{[t_0,t]}$, denoting the state paths until time t; further only those variables y _j for which T _j ≤t; and finally t_1:M. This density is given with respect to the product measure ${\mathcal{L}}_{X_{[t_0,t]}}\prod _{j=1}^M{\mu }_{{\mathcal{Y}}_j}\prod _{j=1}^M{\lambda }_{[t_0,\infty )}$ by:

$$\begin{array}{c}{\widehat{f}}_t^{X_{[t_0,t]},Y_{1:M},T_{1:M}}(x_{[t_0,t]},y_{1:M},t_{1:M})=\prod _{\begin{array}{c}j=1\\ t_j\le t\end{array}}^M{g}_j\left(y_j\right|x_{t_j},t_j\left)\prod _{j=1}^M{\gamma }_j\right(t_j).\end{array}$$

(22)

Note that we cannot use the simple notation of the standard case where for filtering only the first k observations are taken into consideration at time t _k , since neither the observations are ordered in time nor the times T _j are fixed in advance. For this reason we have to include all measurements Y_1:M also into the filter model. Note that even though we use the complete data Y_1:M=y_1:M in the notation, only those y _j have to be known at time t for which T _j ≤t holds. To avoid confusion, we mark all densities connected to the filter model at time t by a hat superscript (and by the index t).

We will now derive formulas for the filter density. Since we assume that the observation times t_1:M are not observed, we use marginalization to get the joint density for $x_{[t_0,t]}$ and Y_1:M only, which is

$$\begin{array}{c}{\widehat{f}}_t^{X_{[t_0,t]},Y_{1:M}}(x_{[t_0,t]},y_{1:M})\hfill \\ ={\int }_{t_0}^{\infty }\cdots {\int }_{t_0}^{\infty }{\widehat{f}}_t^{X_{[t_0,t]},Y_{1:M},T_{1:M}}(x_{[t_0,t]},y_{1:M},t_{1:M})\mathrm{d}{t}_1\cdots \mathrm{d}{t}_M\hfill \\ ={\int }_{t_0}^{\infty }\cdots {\int }_{t_0}^{\infty }\prod _{\begin{array}{c}j=1\\ t_j\le t\end{array}}^M{g}_j\left(y_j\right|x_{t_j},t_j\left)\prod _{j=1}^M{\gamma }_j\right(t_j)\mathrm{d}{t}_1\cdots \mathrm{d}{t}_M\hfill \end{array}$$

(23)

with respect to the product measure ${\mathcal{L}}_{X_{[t_0,t]}}\prod _{j=1}^M{\mu }_{{\mathcal{Y}}_j}$. We will further simplify this density. If we define

$${\bar{g}}_{j,t}\left(y_j|x_{t_j},t_j\right):=\left\{\begin{array}{c}{g}_j\left(y_j\right|x_{t_j},t_j)\text{if}{t}_j\le t\\ 1\text{otherwise},\end{array}\right.$$

then

$$\prod _{\begin{array}{c}j=1\\ t_j\le t\end{array}}^M{g}_j\left(y_j|x_{t_j},t_j\right)=\prod _{j=1}^M{\bar{g}}_{j,t}\left(y_j|x_{t_j},t_j\right)$$

and further

$$\begin{array}{c}{\widehat{f}}_t^{X_{[t_0,t]},Y_{1:M}}(x_{[t_0,t]},y_{1:M})\hfill \\ ={\int }_{t_0}^{\infty }\cdots {\int }_{t_0}^{\infty }\prod _{\begin{array}{c}j=1\\ t_j\le t\end{array}}^M{g}_j\left(y_j|x_{t_j},t_j\right)\prod _{j=1}^M{\gamma }_j\left(t_j\right)\mathrm{d}{t}_1\cdots \mathrm{d}{t}_M\hfill \\ ={\int }_{t_0}^{\infty }\cdots {\int }_{t_0}^{\infty }\prod _{j=1}^M\left({\bar{g}}_{j,t}\left(y_j|x_{t_j},t_j\right){\gamma }_j\left(t_j\right)\right)\mathrm{d}{t}_1\cdots \mathrm{d}{t}_M\hfill \\ =\prod _{j=1}^M{\int }_{t_0}^{\infty }{\bar{g}}_{j,t}\left(y_j\right|x_{t_j},t_j\left){\gamma }_j\right(t_j)\mathrm{d}{t}_j,\hfill \end{array}$$

(24)

where the last step is possible because the factor indexed by j does not depend on $t_{j^{\prime }}$ for j^′≠j. For each j, we can split the integration by T _j at the time point t into two parts and get

$$\begin{array}{c}{\int }_{t_0}^{\infty }{\bar{g}}_{j,t}\left(y_j|x_{t_j},t_j\right){\gamma }_j\left(t_j\right)\hfill \\ ={\int }_{t_0}^t{g}_j\left(y_j|x_{t_j},t_j\right){\gamma }_j\left(t_j\right)\mathrm{d}{t}_j+{\int }_t^{\infty }{\gamma }_j\left(t_j\right)\mathrm{d}{t}_j\hfill \\ =1+{\int }_{t_0}^t\left(g_j\left(y_j\right|x_{t_j},t_j)-1\right){\gamma }_j\left(t_j\right)\mathrm{d}{t}_j\hfill \end{array}$$

where the last step follows from the fact that γ _j is a probability density and therefore

$${\int }_{t_0}^t{\gamma }_j\left(t_j\right)\mathrm{d}{t}_j+{\int }_t^{\infty }{\gamma }_j\left(t_j\right)\mathrm{d}{t}_j={\int }_{t_0}^{\infty }{\gamma }_j\left(t_j\right)\mathrm{d}{t}_j=1$$

holds. Inserting this into (24), we get

$$\begin{array}{c}{\widehat{f}}_t^{X_{[t_0,t]},Y_{1:M}}\left(x_{[t_0,t]},y_{1:M}\right)\hfill \\ =\prod _{j=1}^M\left(1+{\int }_{t_0}^t\left(g_j\left(y_j\right|x_{t_j},t_j)-1\right){\gamma }_j\left(t_j\right)\mathrm{d}{t}_j\right).\hfill \end{array}$$

With a further marginalization, we get the joint density of X _t and Y_1:M for the filter model,

$$\begin{array}{c}{\widehat{f}}_t^{X_t,Y_{1:M}}(x_t,y_{1:M})\hfill \\ :={\int }_{\left\{{\tilde{x}}_{[t_0,t]}\in {\mathcal{X}}_{[t_0,t]}:{\tilde{x}}_t=x_t\right\}}{\widehat{f}}_t^{X_{[t_0,t]},Y_{1:M}}({\tilde{x}}_{[t_0,t]},y_{1:M})\mathrm{d}{\mathcal{L}}_{X_{[t_0,t]}}\left({\tilde{x}}_{[t_0,t]}\right)\hfill \\ ={\int }_{\left\{{\tilde{x}}_{[t_0,t]}\in {\mathcal{X}}_{[t_0,t]}:{\tilde{x}}_t=x_t\right\}}\hfill \\ \prod _{j=1}^M\left(1+{\int }_{t_0}^t\left(g_j\left(y_j\right|{\tilde{x}}_{t_j},t_j)-1\right){\gamma }_j\left(t_j\right)\mathrm{d}{t}_j\right)\mathrm{d}{\mathcal{L}}_{X_{[t_0,t]}}\left({\tilde{x}}_{[t_0,t]}\right)\hfill \end{array}$$

(25)

which is with respect to the product measure ${\mathcal{L}}_{X_t}\prod _{j=1}^M{\mu }_{{\mathcal{Y}}_j}$. From this density, we finally can compute the filter density with respect to ${\mathcal{L}}_{X_t}$ by applying Bayes’ theorem:

$${\widehat{f}}_t^{X_t|Y_{1:M}}\left(x_t\right|y_{1:M})=\frac{{\widehat{f}}_t^{X_t,Y_{1:M}}(x_t,y_{1:M})}{{\widehat{f}}_t^{Y_{1:M}}\left(y_{1:M}\right)}$$

(26)

where

$${\widehat{f}}_t^{Y_{1:M}}\left(y_{1:M}\right)={\int }_{{\mathcal{X}}_t}{\widehat{f}}_t^{X_t,Y_{1:M}}({\tilde{x}}_t,y_{1:M})\mathrm{d}{\mathcal{L}}_{X_t}\left({\tilde{x}}_t\right)$$

(27)

is the data likelihood with respect to the measure $\prod _{j=1}^M{\mu }_{{\mathcal{Y}}_j}$.

Effective computation of the filter distributions

In the following paragraph, we will show how the densities of the filter distributions given by (26) can be effectively computed. This is the basis for the formulation of our MTU particle filter method.

Let the observations $y_{1:M}\in {\mathcal{Y}}_{1:M}$ with $f^{Y_{1:M}}\left(y_{1:M}\right)>0$ be given. For each time t∈ [t₀,∞) and for each j∈{1,…,M}, we define random variables

$$W_{j,t}:\Omega \to {\mathbb{R}}_{\ge 0}\text{and}{W}_t:\Omega \to {\mathbb{R}}_{\ge 0}$$

by the following system of ODEs

$$\begin{array}{ll}\mathrm{d}{W}_{1,t}\left(\omega \right)& =\left(g_1\left(y_1\right|X_t\left(\omega \right),t)-1\right){\gamma }_1\left(t\right)\mathrm{d}t\\ ⋮\\ \mathrm{d}{W}_{M,t}\left(\omega \right)& =\left(g_M\left(y_M\right|X_t\left(\omega \right),t)-1\right){\gamma }_M\left(t\right)\mathrm{d}t\end{array}$$

(28)

for each ω∈Ω with initial values

$$W_{1,t_0}\left(\omega \right)=\cdots =W_{M,t_0}\left(\omega \right)=1,$$

(29)

and by

$$W_t=\prod _{j=1}^M{W}_{j,t}.$$

(30)

We will show that for each set A in the σ-algebra generated by the variable X _t , it holds that

$$\begin{array}{c}{\int }_{X_t\left(A\right)}{\widehat{f}}_t^{X_t|Y_{1:M}}\left(x_t\right|y_{1:M}\left)\mathrm{d}{\mathcal{L}}_{X_t}\right(x_t)=\frac{{\int }_A{W}_t\left(\omega \right)\mathrm{d}\mathsf{\text{P}}\left(\omega \right)}{{\int }_{\Omega }{W}_t\left(\omega \right)\mathrm{d}\mathsf{\text{P}}\left(\omega \right)}\end{array}$$

(31)

where ${\widehat{f}}_t^{X_t|Y_{1:M}}$ is the filter density. That means we can use the processes W_j,t and W _t to compute the filter distributions through time. From this, it follows immediately that we can also compute filter expectations. Indeed, for any real-valued measurable function h on $\mathcal{X}$ such that E[|h(X _t )|]<∞, it holds that the expectation of h(X _t ) given Y_1:M=y_1:M with respect to the filtered state X _t defined by

$$\begin{array}{c}{\widehat{\mathsf{\text{E}}}}_t\left[h\right(X_t\left)\right|Y_{1:M}=y_{1:M}]:={\mathsf{\text{E}}}_{{\widehat{f}}_t^{X_t|Y_{1:M}}(·|y_{1:M})}[h\left(X_t\right)]\hfill \\ ={\int }_{{\mathcal{X}}_t}{\widehat{f}}_t^{X_t|Y_{1:M}}\left(x_t\right|y_{1:M}\left)h\right(x_t\left)\mathrm{d}{\mathcal{L}}_{X_t}\right(x_t),\hfill \end{array}$$

(32)

is given by the following equation:

$$\begin{array}{c}{\widehat{\mathsf{\text{E}}}}_t\left[h\left(X_t\right)|Y_{1:M}=y_{1:M}\right]=\frac{{\int }_{\Omega }{W}_t\left(\omega \right)h\left(X_t\right(\omega \left)\right)\mathrm{d}\mathsf{\text{P}}\left(\omega \right)}{{\int }_{\Omega }{W}_t\left(\omega \right)\mathrm{d}\mathsf{\text{P}}\left(\omega \right)}.\end{array}$$

(33)

To show our assertion, we consider the processes W_j,t for j=1,…,M. According to (28) and (29), each W_j,t is defined as

$$W_{j,t}=1+{\int }_{t_0}^t\left(g_j\left(y_j\right|X_{t_j},t_j)-1\right){\gamma }_j\left(t_j\right)\mathrm{d}{t}_j,$$

(34)

so, with (30),

$$\begin{array}{c}{W}_t=\prod _{j=1}^M{W}_{j,t}=\prod _{j=1}^M\left(1+{\int }_{t_0}^t\left(g_j\left(y_j\right|X_{t_j},t_j)-1\right){\gamma }_j\left(t_j\right)\mathrm{d}{t}_j\right)\end{array}$$

(35)

holds. Thus, to prove (31), we have to show that for each set A from the σ-algebra generated by X _t ,

$$\begin{array}{ll}\frac{{\int }_A{W}_t\left(\omega \right)\mathrm{d}\mathsf{\text{P}}\left(\omega \right)}{{\int }_{\Omega }{W}_t\left(\omega \right)\mathrm{d}\mathsf{\text{P}}\left(\omega \right)}& =\left({\int }_{X_t\left(A\right)}{\int }_{\left\{{\tilde{x}}_{[t_0,t]}\in {\mathcal{X}}_{[t_0,t]}:{\tilde{x}}_t=x_t\right\}}\right.\\ \prod _{j=1}^M\left(1+{\int }_{t_0}^t\left(g_j\left(y_j\right|{\tilde{x}}_{t_j},t_j)-1\right){\gamma }_j\left(t_j\right)\mathrm{d}{t}_j\right)\\ \times \mathrm{d}{\mathcal{L}}_{X_{[t_0,t]}}\left({\tilde{x}}_{[t_0,t]}\right)\mathrm{d}{\mathcal{L}}_{X_t}\left(x_t\right))/{\widehat{f}}_t^{Y_{1:M}}(y_{1:M})\end{array}$$

holds (see (25) and (26)). It is enough to show the equality for the numerator, i.e.

$$\begin{array}{c}{\int }_A{W}_t\left(\omega \right)\mathrm{d}\mathsf{\text{P}}\left(\omega \right)={\int }_{X_t\left(A\right)}{\int }_{\left\{{\tilde{x}}_{[t_0,t]}\in {\mathcal{X}}_{[t_0,t]}:{\tilde{x}}_t=x_t\right\}}\hfill \\ \prod _{j=1}^M\left(1+{\int }_{t_0}^t\left(g_j\left(y_j\right|{\tilde{x}}_{t_j},t_j)-1\right){\gamma }_j\left(t_j\right)\mathrm{d}{t}_j\right)\hfill \\ \times \mathrm{d}{\mathcal{L}}_{X_{[t_0,t]}}\left({\tilde{x}}_{[t_0,t]}\right)\mathrm{d}{\mathcal{L}}_{X_t}\left(x_t\right)\hfill \end{array}$$

since the equality of the denominator follows then immediately from the special case A=Ω and from the fact that

$$\begin{array}{ll}{\widehat{f}}_t^{Y_{1:M}}\left(y_{1:M}\right)=& {\int }_{{\mathcal{X}}_t}{\int }_{\left\{{\tilde{x}}_{[t_0,t]}\in {\mathcal{X}}_{[t_0,t]}:{\tilde{x}}_t=x_t\right\}}\\ \prod _{j=1}^M\left(1+{\int }_{t_0}^t\left(g_j\left(y_j\right|{\tilde{x}}_{t_j},t_j)-1\right){\gamma }_j\left(t_j\right)\mathrm{d}{t}_j\right)\\ \times \mathrm{d}{\mathcal{L}}_{X_{[t_0,t]}}\left({\tilde{x}}_{[t_0,t]}\right)\mathrm{d}{\mathcal{L}}_{X_t}\left(x_t\right).\end{array}$$

Using the variable transformation $\left(X_{[t_0,t]}\right(\omega ),X_t(\omega \left)\right)=({\tilde{x}}_{[t_0,t]},x_t)$, we get

$$\begin{array}{c}{\int }_A{W}_t\left(\omega \right)\mathrm{d}\mathsf{\text{P}}\left(\omega \right)\hfill \\ ={\int }_A\prod _{j=1}^M\left(1+{\int }_{t_0}^t\left(g_j\left(y_j\right|X_{t_j}\left(\omega \right),t_j)-1\right){\gamma }_j\left(t_j\right)\mathrm{d}{t}_j\right)\mathrm{d}\mathsf{\text{P}}\left(\omega \right)\hfill \\ ={\int }_{X_t\left(A\right)}{\int }_{\left\{{\tilde{x}}_{[t_0,t]}\in {\mathcal{X}}_{[t_0,t]}:{\tilde{x}}_t=x_t\right\}}\hfill \\ \prod _{j=1}^M\left(1+{\int }_{t_0}^t\left(g_j\left(y_j\right|{\tilde{x}}_{t_j},t_j)-1\right){\gamma }_j\left(t_j\right)\mathrm{d}{t}_j\right)\hfill \\ \times \mathrm{d}{\mathcal{L}}_{X_{[t_0,t]}}\left({\tilde{x}}_{[t_0,t]}\right)\mathrm{d}{\mathcal{L}}_{X_t}\left(x_t\right).\hfill \end{array}$$

This is what we wanted to show.

Weights

Since for each t and for each j=1,…,M the random variables W_j,t and W _t depend only on the process $x_{[t_0,t]}$ until time t, we can define functions $w_{j,t}:X_{[t_0,t]}\left(\Omega \right)\to {\mathbb{R}}_{\ge 0}$ and $w_t:X_{[t_0,t]}\left(\Omega \right)\to {\mathbb{R}}_{\ge 0}$ by setting for each $x_{[t_0,t]}\in X_{[t_0,t]}\left(\Omega \right)\subset {\mathcal{X}}_{[t_0,t]}$:

$$\begin{array}{c}{w}_{j,t}\left(x_{[t_0,t]}\right):=W_{j,t}\left(\omega \right)\text{and}{w}_t\left(x_{[t_0,t]}\right):=W_t\left(\omega \right)\\ \text{for some}\omega \in \Omega \text{with}{X}_{[t_0,t]}\left(\omega \right)=x_{[t_0,t]}.\end{array}$$

It follows from (34) that

$$w_{j,t}\left(x_{[t_0,t]}\right)=1+{\int }_{t_0}^t\left(g_j\left(y_j\right|x_{t_j},t_j)-1\right){\gamma }_j\left(t_j\right)\mathrm{d}{t}_j$$

(36)

for each j and from (35) that

$$w_t\left(x_{[t_0,t]}\right)=\prod _{j=1}^M{w}_{j,t}\left(x_{[t_0,t]}\right).$$

(37)

The values of W _t will serve as weights in the MTU particle filter. We will call W_j,t the partial weights. Since in each discretization scheme which is applied to solve the integral in the formula (36) for W_j,t the integrand has to be evaluated, we may run into practical problems if we use it as it is written in the formula. If the density $g_j\left(y_j\right|x_{t_j},t_j)$ evaluates to a very small value and if we work with fixed-precision numbers, the subtraction of 1 will result in a value which may be practically equal to -1. If this error accumulates over time, we may end up with wrong values for $w_{j,t}\left(x_{[t_0,t]}\right)$. To reduce the computational error, the integral could be split up in the following way:

$$\begin{array}{l}{w}_{j,t}\left(x_{[t_0,t]}\right)=1-{\int }_{t_0}^t{\gamma }_j\left(t_j\right)\mathrm{d}{t}_j+{\int }_{t_0}^t{g}_j\left(y_j\right|x_{t_j},t_j\left){\gamma }_j\right(t_j)\mathrm{d}{t}_j\\ =1-{\bar{\gamma }}_{j,t}+{\bar{w}}_{j,t}\left(x_{[t_0,t]}\right)\hfill \end{array}$$

(38)

where the cumulative distribution function

$${\bar{\gamma }}_{j,t}:={\int }_{t_0}^t{\gamma }_j\left(t_j\right)\mathrm{d}{t}_j$$

(39)

is independent of $x_{[t_0,t]}$ and where the part

$${\bar{w}}_{j,t}\left(x_{[t_0,t]}\right):={\int }_{t_0}^t{g}_j\left(y_j\right|x_{t_j},t_j\left){\gamma }_j\right(t_j)\mathrm{d}{t}_j$$

(40)

depends on the path $x_{[t_0,t]}$. It is even more convenient to compute ${\bar{\gamma }}_{j,t}$ by evaluating the antiderivative of γ _j , if it is computationally available.

Note that the definition of the filter distribution is dependent on the reference measure ${\mu }_{{\mathcal{Y}}_j}$. A suitable change of this measure may help to further increase the efficiency of the algorithm. This issue still has to be explored.

Resampling

Special attention is needed for the computation of the weights after resampling steps have been applied. As mentioned earlier, resampling at time s _ℓ is done by randomly generating a selection function ι _ℓ :I→I (with index set I={1,…,N}) based on given non-negative (unnormalized) selection weights $v_{s_{\ell }}^i$ for each particle index i. The state samples ${\left(x_{s_{\ell }}^i\right)}_{i=1,\dots ,N}$ have then to be replaced by the selected state samples ${\left(x_{s_{\ell }}^{{\iota }_{\ell }\left(i\right)}\right)}_{i=1,\dots ,N}$, and the unnormalized weights ${\left(w_{s_{\ell }}^i\right)}_{i=1,\dots ,N}$ by the corrected weights ${\left(w_{s_{\ell }}^{{\iota }_{\ell }\left(i\right)}/v_{s_{\ell }}^{{\iota }_{\ell }\left(i\right)}\right)}_{i=1,\dots ,N}$. We assume that resampling steps at times s₁,…,s _ℓ with t₀≤s₁⋯<s _ℓ ≤t have occurred, states and weights have been replaced at these times, and within this paragraph, we denote them by $x_t^{s_1,\dots ,s_{\ell };i}$ and $w_t^{s_1,\dots ,s_{\ell }}\left(x_{[t_0:t]}^{s_1,\dots ,s_{\ell };i}\right)$, respectively, for each particle i. By definition, we have at t=s _ℓ

$$x_{s_{\ell }}^{s_1,\dots ,s_{\ell };i}=x_{s_{\ell }}^{s_1,\dots ,s_{\ell -1};{\iota }_{\ell }\left(i\right)},$$

and

$$w_{s_{\ell }}^{s_1,\dots ,s_{\ell }}\left(x_{[t_0:s_{\ell }]}^{s_1,\dots ,s_{\ell };i}\right)=\frac{w_{s_{\ell }}^{s_1,\dots ,s_{\ell -1}}\left(x_{[t_0:s_{\ell }]}^{s_1,\dots ,s_{\ell };i}\right)}{v_{s_{\ell }}^{{\iota }_{\ell }\left(i\right)}}.$$

For each time t≥s _ℓ and for each particle i, the corrected weights are then recursively given by

$$\begin{array}{c}{w}_t^{s_1,\dots ,s_{\ell }}\left(x_{[t_0:t]}^{s_1,\dots ,s_{\ell };i}\right)=\frac{w_t^{s_1,\dots ,s_{\ell -1}}\left(x_{[t_0:t]}^{s_1,\dots ,s_{\ell };i}\right)}{v_{s_{\ell }}^{{\iota }_{\ell }\left(i\right)}}=\frac{w_t\left(x_{[t_0:t]}^{s_1,\dots ,s_{\ell };i}\right)}{{\prod }_{\lambda =1}^{\ell }{v}_{s_{\lambda }}^{{\iota }_{\lambda :\ell }}}\end{array}$$