The Hitchhiker's Guide to Nonlinear Filtering

TL;DR

This paper provides a comprehensive tutorial on nonlinear filtering, covering static Bayesian inference, discrete-time models, and continuous-time theory, with insights into probability measures and particle filtering algorithms.

Contribution

It offers an accessible overview connecting various aspects of nonlinear filtering, serving as an introduction for scientists and engineers to advanced literature.

Findings

Unified view of discrete and continuous filtering

Insights into probability measure changes in filtering

Guidance on particle filtering algorithm construction

Abstract

Nonlinear filtering is the problem of online estimation of a dynamic hidden variable from incoming data and has vast applications in different fields, ranging from engineering, machine learning, economic science and natural sciences. We start our review of the theory on nonlinear filtering from the simplest `filtering' task we can think of, namely static Bayesian inference. From there we continue our journey through discrete-time models, which is usually encountered in machine learning, and generalize to and further emphasize continuous-time filtering theory. The idea of changing the probability measure connects and elucidates several aspects of the theory, such as the parallels between the discrete- and continuous-time problems and between different observation models. Furthermore, it gives insight into the construction of particle filtering algorithms. This tutorial is targeted at scientists and engineers and should serve as an introduction to the main ideas of nonlinear filtering, and as a segway to more advanced and specialized literature.

Equations

$$p(X|Y)=\frac{p(Y|X)p(X)}{p(Y)}.$$

$${\mathbb{E}}_{\mathbb{P}}\left[\varphi (X)\right]$$

$${\mathbb{E}}_{\mathbb{P}}[\varphi (X)|Y]$$

$$L(x,y)=\frac{p(x,y)}{q(x,y)}=\frac{p(y|x)p(x)}{p(x)q(y)}=\frac{p(y|x)}{q(y)},$$

$${\mathbb{E}}_{\mathbb{P}}[\varphi (X)]$$

$${\mathbb{E}}_{\mathbb{P}}[\varphi (X)|Y]$$

$$Z=\sum _{i}p(Y|X^i)=\sum _i{w}_i.$$

$$p(X_{t_n}|X_{t_{0:n-1}})$$

$$p(X_{t_n}|Y_{0:t_n})$$

$$=$$

$$p(X_{t_n}=0|X_{t_{n-1}}=0)=\alpha ,p(X_{t_n}=1|X_{t_{n-1}}=1)=\beta .$$

$$P^{\top }=\left(\begin{array}{cc}\alpha & 1-\beta \\ 1-\alpha & \beta \end{array}\right),$$

$$p_{t_n}$$

$$p(Y_{t_n}=1|X_{t_n}=0)=\delta ,p(Y_{t_n}=0|X_{t_n}=1)=\delta .$$

$$p(X_{t_n}|Y_{t_0:t_n})$$

$$e(Y_{t_n})$$

$$p(X_{t_n}|X_{t_{n-1}})$$

$$p(Y_{t_n}|X_{t_n})$$

$$p(X_{t_{n-1}}|Y_{0:t_{n-1}})$$

$${\mu }_{t_n}$$

$${\Sigma }_{t_n}$$

$$K_t$$

$${\mathbb{E}}_{\mathbb{P}}\left[\varphi (X_{t_n})|Y_{0:t_n}\right]$$

$$L(x_{0:t_n},y_{0:t_n})$$

$$q(x_{0:t_n}|y_{0:t_n})$$

$$p(x_{0:t_n}|y_{0:t_n})$$

$$L(x_{0:t_n},y_{0:t_n})$$

$$X_{t_n}^{(i)}$$

Full text

††journal: Journal of Mathematical Psychology

1 Introduction: A Guide to the Guide

*“The introduction begins like this:

Space, it says, is big. Really big.

You just won’t believe how vastly hugely mind-bogglingly big it is.”

*— Douglas Adams

Filtering is the problem of estimating a dynamically changing state, which cannot be directly observed, from a stream of noisy incoming data. To give a concrete example, assume that you are a Vogon in charge of a spaceship. Since you had a particularly bad day, you decide to destroy a small asteroid to make you feel better. Before you push the red button, you need to know the current position of the asteroid, which corresponds to the hidden state $X_{t}$. You have some idea about the physics of movement in space, but there is also a stochastic component in the movement of your target. Overall, the asteroid’s movement is described by a stochastic dynamical model. In addition, you cannot directly observe its position (because you like to keep your safe distance), so you have to rely on your own ship’s noisy measurements $Y_{t}$ of the position of the asteroid. Because of these uncertainties, it would not only be useful to have an estimate of the asteroid’s current position $X_{t}$ based on the history of measurements $Y_{0:t}=\{Y_{0},Y_{1},\dots,Y_{t}\}$, but also an estimate of the uncertainty of the estimate. Thus generally, the conditional probability density $p(X_{t}|Y_{0:t})$ is the complete solution to your problem (and the beginning of the problem of how to find this solution).

These sorts of problems are not only relevant for bad-tempered Vogons, but in fact are encountered in a wide variety of applications from different fields. Initial applications of filtering were centered mostly around engineering. After the seminal contributions to linear filtering problems by Kalman [1960], Kalman & Bucy [1961], the theory was largely applied to satellite orbit determination, submarine and aircraft navigation as well as space flight [Jazwinski, 1970]. Nowadays, applications of (nonlinear) filtering range from engineering, machine learning [Bishop, 2006], economic science (in particular mathematical finance, some examples are found in Brigo & Hanzon, 1998) and natural sciences such as geoscience [Van Leeuwen, 2010], in particular data assimilation problems for weather forecasting, neuroscience and psychology. As a particular example for its usefulness in neuroscience, the modeling of neuronal spike trains as point processes [Brillinger, 1988, Truccolo, 2004] has led to interesting filtering tasks, such as the problem of decoding a stimulus from the spiking activity of neurons (e.g. Koyama et al., 2010, Macke et al., 2011). In psychology, nonlinear filtering techniques are not only used for data analysis, but can also provide qualitative insight into psychological processes such as perception [Wolpert et al., 1995, Körding et al., 2007] or decision making [Drugowitsch et al., 2014, Glaze et al., 2015, Radillo et al., 2017, Veliz-Cuba et al., 2016, Piet et al., 2018]. To tackle these kinds of questions, knowledge about nonlinear filtering is indispensable. Theoretical understanding can further help in connecting and unifying specific applications of filters and is paramount for understanding more advanced topics in filtering [Jazwinski, 1970, Section 1.2].

The aim of this tutorial is to present – in an easily accessible and intuitive way – the theoretical basis for continuous-time nonlinear filtering with diffusion and point-process observations. The tutorial highlights the change of measure as a powerful tool to derive the fundamental equations of nonlinear filtering as well as numerical approximations. In addition, the unification provided by the concept of change of measure provides a solid basis for diving into the huge body of literature on nonlinear filtering. Our tutorial complements the more advanced theoretical literature (e.g. Jazwinski, 1970, Bain & Crisan, 2009, Bremaud, 1981) or more specialized tutorials, e.g. on particle filtering [Arulampalam et al., 2002, Doucet & Johansen, 2009, Speekenbrink, 2016], Hidden Markov Models [Visser, 2011, Rabiner, 1989] or variational Bayes for latent linear systems [Ostwald et al., 2014].

2 A view from space: from Bayes’ rule to filtering

*“Even the most seasoned star tramp can’t help but shiver

at the spectacular drama of a sunrise seen from space,

but a binary sunrise is one of the marvels of the Galaxy.”

*— Douglas Adams

Suppose that we observe a random variable $Y$ and want to infer the value of an (unobserved) random variable $X$. Bayes’ rule tells us that the conditional distribution of $X$ given $Y$, the so-called posterior, can be computed in terms of three ingredients: the prior distribution $p(X)$, the likelihood $p(Y|X)$, and the marginal likelihood $P(Y)$ (which acts as a normalizing constant):

[TABLE]

This tutorial is concerned with the application of the above idea to a situation where $X,Y$ are continuous-time stochastic processes and we want to perform the inference online as new data from $Y$ comes in. In this section, we want to gradually build up the stage: as some readers might be more familiar with discrete-time filtering due to its high practical relevance and prevalence, we will start our journey from there, picking up important recurring concepts as we make our way to continuous-time filtering.

2.1 Changes of measure

Before we talk about dynamic models, let us briefly highlight a concept in Bayesian inference that will be very important in the sequel: that of changing the probability measure. A probability measure is a function that assigns numbers (‘probabilities’) to events. If we have two such measures $\mathbb{P}$ and $\mathbb{Q}$, then $\mathbb{P}$ is called absolutely continuous wrt. $\mathbb{Q}$ if every nullset of $\mathbb{Q}$ is a nullset of $\mathbb{P}$. Moreover, $\mathbb{P}$ and $\mathbb{Q}$ are called equivalent if they have the same nullsets. In other words, if $A$ denotes an event, and $P(A)$ denotes its probability, then equivalence means that $\mathbb{Q}(A)=0$ if and only if $\mathbb{P}(A)=0$.

But why would we want to change the measure in the first place? Changing the measure allows us to compute expectations of a measurable function $\phi(x)$ with respect to a measure $\mathbb{Q}$, which were originally expressed with respect to another measure $\mathbb{P}$. To see this, consider the two measures $\mathbb{P}$ and $\mathbb{Q}$ for some real-valued random variable $X$, and write them in terms of their densities $p,q$ (with respect to the Lebesgue measure).111In this section, all ‘densities’ are with respect to the Lebesgue measure. We then have

[TABLE]

where we introduced the likelihood ratio $L(x):=\frac{p(x)}{q(x)}$ and $\mathbb{E}_{\mathbb{Q}}$ denotes expectation under the distribution $q$. Thus, changing the measure proves to be very useful whenever expectations under $\mathbb{Q}$ are easier to compute than under $\mathbb{P}$.

A fundamental problem in filtering is that of computing a conditional expectation (i.e. an expected value under the posterior distribution) of this sort:

[TABLE]

for some function $\phi$, and we want to use Eq. (1) to compute $p(X|Y)$. We therefore have to compute the two integrals here

[TABLE]

but the structure of the model (interactions between $X$ and $Y$) might make it very hard to compute the integrals, either analytically or numerically. Thus, we again change the measure to a reference measure $\mathbb{Q}$ with joint density $q(x,y)$, and rewrite Eq. (4):

[TABLE]

where now the likelihood ratio $L(x,y)=\frac{p(x,y)}{q(x,y)}$ is a function of both $x$ and $y$.

The hope is that we can pick a reference measure $\mathbb{Q}$ such that both $L(x,y)$ and $q(x,y)$ are simple enough to make Eq. (5) more tractable than Eq. (3). For instance, some simplification might be achieved by switching from a model $p(x,y)$ of $\mathbb{P}$ in which $X$ and $Y$ are coupled, i.e. statistically dependent, to a model $p(x)q(y)$ of $\mathbb{Q}$ where they are independent (while preserving the distribution of $X$), i.e. under model $\mathbb{Q}$ we find $q(x,y)=p(x)q(y)$. A potential added advantage of changing measure is when the distribution $q(y)$ is computationally simple. Then, the likelihood ratio $L(x,y)$ reads

[TABLE]

and conditional expectations under $\mathbb{Q}$ can simply be taken with respect to the prior probability $p(x)$.

Please take a moment to appreciate the value of this idea: the change of measure has allowed us to replace the expectation with respect to the posterior $p(x|y)$ of $\mathbb{P}$ (which might be hard to get) with an expectation with respect to the prior $p(x)$ of $\mathbb{Q}$ (which might be easy to compute). This ‘trick’ will become the central theme of this manuscript.

2.1.1 Importance sampling

A numerical example where a ‘measure change’ is directly used is importance sampling. Here, the goal is to approximate expectations with respect to a distribution $p(x)$ (under $\mathbb{P}$) by using $M$ empirical samples $X^{i}\sim p(x)$, such that

[TABLE]

However, there might be situations where we cannot draw samples from $p(x)$, but only from another distribution $q(x)$ (under $\mathbb{Q}$). Thus, we first perform a change of measure to $\mathbb{Q}$, and then use the samples $X^{i}\sim q(x)$ to approximate the expectation:

[TABLE]

In this context, the likelihood ratio $L(X^{i})=\frac{p(X^{i})}{q(X^{i})}=:w_{i}$ is referred to as (unnormalized) importance weight. Hence, the target distribution $p(x)$ is not only represented by the density of empirical samples (or ‘particles’), i.e. how many samples can be found in a specific interval in the state space, but also by their respective importance weights (see simple example in Figure 1).

Similarly, we can use importance sampling to approximate a posterior expectation $\mathbb{E}_{\mathbb{P}}[\phi(X)|Y]$ with samples from the prior distribution $p(x)$. For this, consider changing to a measure $\mathbb{Q}$ with density $q(x,y)=p(x)q(y)$, such that with the likelihood ratio in Eq. (6) we find

[TABLE]

where the unnormalized importance weights are given by $w_{i}=p(Y|X^{i})$, and the normalization $Z$ is given by

[TABLE]

Thus, in order to approximate a posterior with empirical samples from the prior, each sample $X^{i}\sim p(x)$ has to be weighted according to how likely it is that this particular sample has generated the observed value of the random variable $Y$ by evaluating the likelihood $p(Y|X^{i})$ for this sample. Generalizing this to a dynamical inference setting will give rise to the bootstrap particle filter, and we will show in Section 6 how changing the measure in empirical sampling for a dynamical system results in dynamical equations for the particles and weights.

2.2 Filtering in discrete time - an introductory example

The inference problems in the previous section were of purely static nature.222If you ask yourself why we needed 2.2 pages to get to this point, please bear with us: the concept of changing the measure is very straightforward in a static setting, and might help to grasp the (seemingly) more complicated applications in a dynamic setting later on. However, this being the Hitchhiker’s guide to nonlinear filtering, let us now start to consider dynamical models for filtering.

Filtering means computing the conditional distribution of the hidden state $X_{t}$ at time $t$ using the observations up to that time $Y_{0:t}=\{Y_{0},\dots,Y_{t}\}$. There are two important ingredients to this problem: first, the signal model describes the dynamics of the hidden state $X_{t}$. In order to perform the inference recursively, the usual minimum assumption for the hidden, or latent, state process $X_{t}$ with state space $S$ is that it is a first-order Markov process, which, roughly speaking, means that the (probability of the) current state just depends on the last state, rather than on the whole history. In discrete time,333In discrete time $t_{n}=n\Delta t$. we can write more formally

[TABLE]

Thus, the dynamics of the whole process is captured by the transition probability $p(X_{t_{n}}|X_{t_{n-1}})$, which is assumed to be known.

Second, the observation model describes the (stochastic) generation of the observation process $Y_{t_{n}}$, and is captured by the emission probability $p(Y_{t_{n}}|X_{t_{n}})$, which is also assumed to be known. Together, the transition and emission probability form a so-called state space model (SSM).444Somewhat oddly, the name ‘state space model’ usually refers to a model with continuous state space, i.e. $X_{t}\in\mathbb{R}^{n}$, which is distinct from models with finite state space such as the Hidden Markov Model below. In this tutorial, the state space can both be discrete or continuous, and if necessary, will be further clarified in the text. With these ingredients, the filtering problem in discrete time reduces to a simple application of Bayes’ rule (Eq. 1) at each time step, which may be written recursively:

[TABLE]

The simplest dynamic model for filtering is a Hidden Markov Model (HMM). To see that you don’t need rocket science for hitchhiking and applying Eq. (13), let us consider an HMM with two hidden states and two observed states, i.e. $X_{t_{n}}$ and $Y_{t_{n}}$ can take values of 0 or 1 for each time $t_{n}$. The transition probabilities for $X_{t_{n}}$ are given by

[TABLE]

Thus, $\alpha$ is the probability of staying in state 0, whereas $\beta$ is the probability of staying in state $1$, and leaving those states has to have probability $1-\alpha$ and $1-\beta$, respectively, where we assume that $0<\alpha,\beta<1$ such that each state is visited. This can be represented by a matrix.555Here $P^{\top}$ is used to denote the transpose of the matrix $P$. It will become clear later why we define the transition matrix as $P^{\top}$.

[TABLE]

which recursively determines the distribution of the hidden Markov chain at each time: if $p_{t_{n-1}}=(p_{t_{n-1}}^{(1)},p_{t_{n-1}}^{(2)})^{\top}$ is a two-dimensional vector, denoting probability of state occupancy at time $t-1$, i.e. $p_{t_{n-1}}^{(1)}=P(X_{t_{n-1}}=0)$ and $p_{t_{n-1}}^{(2)}=P(X_{t_{n-1}}=1)$, the corresponding vector at time $t_{n}$ is given by

[TABLE]

In our example, the emission probabilities of $Y$ are given by a binary symmetric channel (random bit flip) with error probability $0<\delta<1$

[TABLE]

The structure of this model is depicted in Figure 2a.

For filtering, we can directly apply Eq. (13), and since the state space is discrete, the integral reduces to a sum over the possible states [math] and $1$. Thus, Eq. (13) may be expressed as

[TABLE]

where $\text{diag}(v)$ is a diagonal matrix with the vector $v$ along the diagonal and $e$ is a vector encoding the emission likelihood:

[TABLE]

Figure 2b shows a sample trajectory of the hidden state $X_{t_{n}}$, the corresponding observations $Y_{t_{n}}$ as well as the filtered probabilities $p_{t_{n}}$ and the estimated state $\hat{X}_{t_{n}}$. Even though what is presented here is a very simple setting (discrete time finite number of states), it illustrates nicely that the filter takes into account both the dynamics of the hidden states as well as the reliability of the observations.

2.3 Continuous (state) space

Remarkably, in the previous example the filtering problem could be solved in closed form because it was formulated in discrete time for a discrete state space. We will now continue our journey towards more complex filtering problems involving a continuous state space. For this, we have to go back to Eq. (13), which is actually the filtering recursion for any state space. While being straightforward to write down - is it possible to solve it in closed form? Depending on the specific transition and emission densities, the integrals in Eq. (13) might not admit a closed-form solution. In fact, this is almost always the case! Except…

2.3.1 The Kalman filter

… if the transition and emission probabilities are Gaussians and linear, i.e.

[TABLE]

where we consider $X_{t_{n}}\in\mathbb{R}^{k}$ and $Y_{t_{n}}\in\mathbb{R}^{l}$ to be vector-valued random processes. Further, $A\in\mathbb{R}^{k\times k}$ and $B\in\mathbb{R}^{l\times k}$ are the transition and emission matrices, respectively, and $\Sigma_{x}\in\mathbb{R}^{n\times n}$ and $\Sigma_{y}\in\mathbb{R}^{l\times l}$ are state and observation noise covariances, respectively.

Let us assume that at time $t_{n-1}$ the posterior is given by a Gaussian

[TABLE]

We can immediately plug Eq. (20) and (21) together with this assumption into Eq. (13). After a bit of tedious but straightforward algebra [see Bishop, 2006, Section 13.3.1], we find that the posterior is also a Gaussian $\mathcal{N}(X_{t};\mu_{t_{n}},\Sigma_{t_{n}})$. The famous Kalman filter equations give us update rules for its mean and variance:

[TABLE]

where $\tilde{\Sigma}_{t_{n-1}}=A\Sigma_{t_{n-1}}A^{\top}+\Sigma_{x}$ is the variance of $p(X_{t}|Y_{0:{t_{n-1}}})$, obtained after performing the marginalization over the state transition. The so-called Kalman gain $K_{t}$ is given by

[TABLE]

The immediate implication of this result is that for this particular model, given that the initial distribution is a Gaussian, the posterior stays Gaussian at all times.

2.3.2 Particle filtering in discrete time

In those cases where transition and emission probabilities are not Gaussian, we cannot expect Eq. (13) to take an analytically accessible form. In other words: as time goes by (in terms of time steps $n$), we will have to keep track of an ever-growing amount of integrals, which is clearly not desirable. Alternatively, we can try to approach this task numerically, by considering empirical samples and propagating these samples through time to keep track of this posterior. This idea is the very basis of particle filters (PF).

The only remaining problem is that direct sampling from the true posterior is usually not possible. In Section 2.1.1 we have motivated importance sampling for a static setting from a change of measure perspective, and we will now use the same reasoning to motivate sequential importance sampling. In other words: we will replace samples from the true posterior (under $\mathbb{P}$) by weighted samples from a proposal density under $\mathbb{Q}$. Importantly, a ‘sample’ $i$ here refers a single realization of the whole path $X_{0:t_{n}}=\{X_{0},\dots,X_{t_{n}}\}$, and the measure change needs to be done with respect to the whole sequence of state and observations.

Let us first note that the posterior expectation can be understood as an expectation with respect to the whole sequence

[TABLE]

where in the last step we simply used that $\int_{S}dx_{0:t_{n-1}}\,p(x_{0:t_{n-1}}|Y_{0:t_{n}})=1$. Now, we perform the measure change according to Eq. (5) :

[TABLE]

with

[TABLE]

where $p$ and $q$ denote densities of $\mathbb{P}$ and $\mathbb{Q}$, respectively. Let us now choose the measure $\mathbb{Q}$ such that the conditional density $q(x_{0:t_{n}}|y_{0:t_{n}})$ factorizes, i.e.

[TABLE]

Further, we can rewrite the conditional density $p(x_{0:t_{n}}|y_{0:t_{n}})$ using the structure of the SSM

[TABLE]

Thus, using that all factors independent of the state variable can be taken out of the expectations in Eq. (28) and cancel subsequently, we find

[TABLE]

In analogy to Section 2.1.1, we now take $M$ i.i.d. samples from the proposal density, i.e. we sample $X^{i}_{0:t_{n}}\sim q(X_{0:t_{n}}|Y_{0:t_{n}})$, and weigh them according to the value of the likelihood ratio evaluated at the particle positions (cf. Eq. 9). Since the proposal in Eq. (29) was chosen to factorize, both the sampling process as well as the evaluation of the unnormalized importance weights $w_{t_{n}}^{(i)}$ (according to Eq. 31) can be done recursively. More specifically, the problem of sampling (and weighing) the whole sequences $X^{(i)}_{0:t_{n}}$ is replaced by sampling just a single transition $X_{t_{n}}^{(i)}$ for each of the $M$ particles at each time step $n$ and updating the associated particle weights.

[TABLE]

such that the posterior expectation is approximated by

[TABLE]

with $Z_{t_{n}}=\sum_{i=1}^{P}w_{t_{n}}^{(i)}$.

A simple (but not necessarily efficient) choice is to use the transition probability $p(X_{t_{n}}|X_{t_{n-1}})$ as the proposal function in Eq. (32). Then, computation of the unnormalized weights simplifies to

[TABLE]

This scheme is the basis of the famous Bootstrap PF (BPF, Gordon et al., 1993)666Although technically, the BPF requires a resampling step at every iteration step.. Doucet et al. [2000] state that the BPF is “inefficient in simulations as the state space is explored without any knowledge of the observations”. To account for this, alternative proposal densities can be crafted in discrete time, which may take into account the observations in the particle transitions (e.g. the ‘optimal proposal’ in Doucet et al., 2000).

3 Knowing where your towel is: setting the stage for continuous-time models

*“A towel is about the most massively useful thing an interstellar hitchhiker can have.

Partly it has great practical value. More importantly,

a towel has immense psychological value.”

*— Douglas Adams

So far, we have made our journey from Bayes’ theorem to discrete-time filtering first for discrete state spaces and then made the transition towards continuous state space models. The next logical step would be the transition to continuous time models. In the following three sections, we will see that the mindset is very similar to the approaches taken before, just in their respective continuous-time limit, i.e. $dt=t_{n}-t_{n-1}\to 0$. In particular, we will use the change of measure approach to derive the filtering equations, i.e. dynamical equations for the posterior expectations $\mathbb{E}[\phi(X_{t})|Y_{0:t}]$ or, equivalently, the posterior density $p(X_{t}|Y_{0:t})$.

But let us take a step back here and first explain the model assumptions under which we will present continuous-time filtering theory. For the purpose of this tutorial, we have seen that a generative model consists of two parts:

A signal model or hidden process model that describes the dynamics of some system whose states we want to estimate. In continuous-time, we will consider two very general classes of signal model, namely continuous-time Markov chains (countable or finite state space) and jump-diffusion processes (continuous state space). 2. 2.

An observation model that describes how the system generates the information that we can observe and utilize in order to estimate the state. We will elaborate the filtering theory for two types of observation noise, namely continuous-time Gaussian noise and Poisson noise.

3.1 Signal models

As in Section 2, we will restrict ourself to the treatment of Markovian processes for the signal, i.e. $p(X_{t}|X_{0:t-dt})=p(X_{t}|X_{t-dt})$. Our goal in this subsection will be to obtain dynamical equations that fully describe the signal process.

3.1.1 Markov chain

An important example is when $X_{t}$ is a continuous-time time-homogeneous Markov chain with a finite number of states, i.e. $S=\{1,...,m\}$. In this case we may represent the function $\phi:\{1,...,m\}\to\mathbb{R}$ as a vector $\phi=(\phi(1),...,\phi(m))^{\top}$ and we have a transition probability matrix $P(t)^{\top}$. The entry $P_{ji}(t)$ gives the probability of going from state $j$ to state $i$ within a time interval of length $t$, so it is a time-dependent generalization of Eq. (15). This allows us to compute the distribution at time $t$, $p(t)$ from the initial distribution $p(0)$ as $p(t)=P(t)^{\top}p(0)$. We therefore have two equivalent ways of computing the expectation of $\phi$:

[TABLE]

In the first one, the observable is fixed while the distribution changes as a function of time, while in the second, the distribution is fixed to the initial distribution, and the observable evolves in time, i.e. $\phi(t)=P(t)\phi$.

By differentiating with respect to time, we obtain differential equations for the distribution $p(t)$ and the observable $\phi(t)$,

[TABLE]

The Markov property ensures that $P(t+s)=P(t)P(s)=P(s)P(t)$. Further, since $P(0)=\mathbb{I}$ is the unit matrix777Because $p(t=0)=P(0)^{\top}p(0)$ is only fullfilled if $P(0)=\mathbb{I}$., the time derivative of the matrix $P(t)$ can be simplified to

[TABLE]

We denote $A=\dot{P}(0)$ and then get

[TABLE]

Equivalently, we find for the time derivative of the expectation

[TABLE]

So conceptually, the whole temporal evolution of the stochastic process $X_{t}$ is encapsulated in the matrix $A$, the so-called generator matrix. In other words, the generator matrix, together with the initial distribution, is all we need to completely characterize the Markov chain process.

3.1.2 Jump-diffusion process

Intuitively, in order to make the transition to a continuous state space, we have to exchange “sums by integrals and matrices by linear operators”. We will now see that this holds for the hidden state dynamics by characterizing a continuous-time stochastic process with continuous state space $S$ similarly to the equations (41) and (42) above.

An important signal model, which is a generalization of the classical diffusion model in continuous time, is a hidden state $X_{t}$ that is a jump-diffusion process, i.e. it evolves according to a stochastic differential equation (SDE) in $S=\mathbb{R}^{n}$,

[TABLE]

Here, $f:\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{n}$, $G:\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{n\times n}$, and $J:\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{n\times k}$ are called the drift, diffusion, and jump coefficients of $X_{t}$, respectively. The process noise is modelled by two types of noise sources: $W_{t}\in\mathbb{R}^{n}$ is a vector Brownian motion that models white Gaussian noise in continuous time, and we may consider $dW_{t}\sim\mathcal{N}(0,dt)$. $N_{t}$ is a $k$-dimensional point process with $k$-dimensional rate (or intensity) vector $\lambda(X_{t})$, i.e. $dN_{t}\sim\text{Poisson}(\lambda(X_{t})dt)$. Note that $dN_{t}$ takes only values 0 or 1, because in the limit $dt\to 0$, the Poisson distribution becomes a Bernoulli distribution. In Figure 3, we show example trajectories from Eq. (43), one being a drift-diffusion (where the jump term vanishes), one being a pure jump process, and the last one being a jump-diffusion process.

Dealing with this type of SDE model is considerably more technical than the Markov chains above. Therefore, we will outline the theory of diffusion processes for readers who are new to them. Unless stated otherwise, derivations presented here roughly follow Gardiner [2009] and Bremaud [1981].

We can choose to describe the process in terms of transition probability densities $p(x,t|x^{\prime},s)$, which give the probability density at a point $x\in\mathbb{R}^{n}$ at time $t$ conditioned on starting at a point $x^{\prime}\in\mathbb{R}^{n}$ at time $s<t$. This transition density can be combined with the initial density $p_{0}(x^{\prime})$ by integrating in order to compute an expectation:

[TABLE]

in complete analogy with the vector-matrix-vector product for the Markov chains in Eq. (36). Taking this analogy further, differentiating with respect to time gives rise to two different (equivalent) ways of writing the time evolution of the expected value:

[TABLE]

where in analogy to Eq. (42) we have introduced the adjoint operator $\mathcal{A}^{{\dagger}}$ that describes the time evolution of the probability density. Thus, in analogy to Eq. (37) we can set out to find the appropriate form of the generator $\mathcal{A}$, which generalizes the generator matrix $A$ of the Markov chain, and then, by integration by parts, we may derive the corresponding adjoint operator $\mathcal{A}^{{\dagger}}$.

Itô lemma for jump diffusions

The form of the generator $\mathcal{A}$ can be obtained by changing the variables in Eq. (43) from the random variable $X_{t}$ to the random variable $\phi_{t}:=\phi(X_{t})$. The following calculation will be performed for a scalar process $X_{t}$.888Generalization to a multivariate state process $X_{t}$ is straightforward. Consider an infinite Taylor expansion of its increment $d\phi_{t}$ around $dX_{t}=0$ up to $\mathcal{O}(dt)$:

[TABLE]

with $\phi_{t}^{(n)}:=\left(\partial_{x}^{n}\phi(x)\right)|_{x=X_{t}}$.

In a deterministic differential, Taylor-expanding up to first order would suffice since $dt^{n}=0$ $\forall n>1$. In Eq. (43), the additional stochastic terms add additional orders of $dt$. Particularly, since the variance of the Brownian motion process grows linearly in time, we have $dW_{t}^{2}=dt$, and thus the diffusion term has to be expanded up to second order. For the jump term, all order up to infinity have to be considered: $N_{t}$ is not a continuous process, with jumps of always size 1 irrespectively of the infinitesimally small time interval $dt$. Therefore, any power of this jump will have the same magnitude, i.e. $dN_{t}^{n}=dN_{t}$, $\forall n$. Thus, we find for a scalar process

[TABLE]

This formula is called Itô lemma. In the last step, we have defined

[TABLE]

where $\mathcal{A}$ is the infinitesimal generator of the stochastic process $X_{t}$. The stochastic process $M_{t}^{\phi}$ is a so-called martingale999Loosely speaking, a martingale is a sequence of random variables, whose conditional expectation in the next time step is equal to the value of the random variable at the current time step. and the contribution from its increment vanishes upon taking expectations, i.e. $\mathbb{E}[M_{t}^{\phi}]=0$. Thus, taking expectations on both sides of Eq. (47) we find indeed

[TABLE]

which is the continuous state space analogue to Eq. (42).

The multivariate version is completely analogous:

[TABLE]

where now the infinitesimal generator of the stochastic process is given by

[TABLE]

and the martingale part reads

[TABLE]

The operator $\mathcal{A}$, just like the generator matrix $A$ of the Markov chain, together with the initial distribution, completely characterizes the Markov process and allows us to describe its time evolution on an abstract level. Or in other words: even though the particular form of $\mathcal{A}$ might be different for each of these models presented here, the structure of the mathematics remains the same, and can therefore be generalized to arbitrary $\mathcal{A}$ when the need arises.

The evolution of the probability density

With the explicit form in Eq. (51) of the generator $\mathcal{A}$, we can go back to Eq. (45) and perform the integration by parts to find the adjoint operator $\mathcal{A}^{{\dagger}}$, which will take the role of $A^{\top}$ in the Markov chain case.

Plugging Eq. (51) into Eq. (45), we obtain

[TABLE]

The first two integrals can be dealt with by ordinary integration by parts101010Here, we make the assumption that the density $p(x,t)$ and all its derivatives vanish at infinity., i.e.

[TABLE]

and

[TABLE]

For the third integral in Eq. (53), we perform a change of variables $x\mapsto x+J_{i}(x,t)$ (where we assume the integral boundaries to not be affected by the substitution), thereby obtaining

[TABLE]

where $\det{\frac{\partial(J_{1i},\dots,J_{ji})}{\partial(x_{1},\dots,x_{j})}}$ denotes the Jacobi determinant of the entries of the column vector $J_{i}$. Combining all of these integrals (including the sums) and comparing with Eq. (45), we can therefore read off the form of the adjoint operator:

[TABLE]

Using Eq. (45) once more, we find the evolution equation for the density $p(x,t)$:

[TABLE]

If we leave out the jump terms, this is called the Fokker-Planck equation or Kolmogorov forward equation. With the jump terms, it is often referred to as the Master equation.

3.2 Observation model

In the previous section, we have encountered various signal models, which are the processes we want to infer. The knowledge about how these processes evolve in time, formally given by the generator $\mathcal{A}$, serves as prior knowledge to the inference task. Equally important, we need measurements, or observations, to update this prior knowledge. In particular, an observation model describes how the signal gets corrupted during measurement. This may comprise both a lossy transformation (e.g. only certain components of a vector-valued process are observed), and some stochastic additive noise that randomly corrupts the measurements. Roughly speaking, the measurements $Y_{t}$ are given by

[TABLE]

but we will need to be careful about the precise way in which the noise is added in order to make sense in continuous time.

In the following, we will consider two types of noise: Gaussian and Poisson. The simplicity of noise of these two noise models greatly simplifies the formal treatment of the filtering problem, and while the two types of noise seem very different, there is a common structure that will emerge.

When considering more general noise models than the ones below, the technique of Section 4 (change of measure) can be applied whenever the observation noise (whatever is added to the deterministic transformation) is additive and independent of the hidden state.

3.2.1 Continuous-time Gaussian noise

The simplest noise model is often white Gaussian noise. For continuous-time observations, however, one cannot simply take an observation model $Y_{t}=h(X_{t})+\eta_{t}$ with independent Gaussian $\eta_{t}$ because for a reasonably well-behaved process $X_{t}$, an integration of $Y_{t}$ over a finite time interval would completely average out the noise and therefore allow one to perfectly recover the transformed signal $h(X_{t})$.111111If the observations are made at discrete times $t_{1},t_{2},...$, this is not problematic. Filtering of a continuous-time hidden process with discrete-time observations is reviewed in Jazwinski [1970, p. 163ff]. If the observation model in Eq. (59) is discretized, one gets back to a discrete-time observation model with Gaussian noise. The filtering problem would therefore be reduced to simply inverting $h$.

One way of resolving the problem of finding a (nontrivial) model of white Gaussian noise is to switch to a differential form and use increments of the Wiener process as a noise term. One therefore obtains an SDE for the observation process $Y_{t}$:

[TABLE]

Here, $h:\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{l}$ is a vector-valued function that links the hidden state (and time, if time-dependence is explicit) with the deterministic drift of the observations. Further, $V_{t}\in\mathbb{R}^{l}$ is a vector Brownian motion process and $\Sigma_{y}(t):\mathbb{R}\to\mathbb{R}^{l\times l}$ is the time-dependent observation noise covariance.

In the standard literature, one usually finds the special case $\Sigma_{y}=\mathbb{I}^{m\times m}$, which is equivalent to Eq. (59) if the increment of the observation process $Y_{t}$ is rescaled accordingly:

[TABLE]

where $\tilde{h}(x,t)=\Sigma_{y}(t)^{-1/2}h(x,t)$ is the rescaled observation function.

3.2.2 Poisson noise

In many fields, observations come in the form of a series of events. Examples include neuroscience (neural spike trains), geoscience (earthquakes, storms), financial transactions, etc. This suggests a point process (whose output is a series of event times) or counting process (which counts the number of events) observation model. A simple, but versatile, model for events is a Poisson process $N_{t}$ with time-varying and possibly history-dependent intensity. As an observation model, this doubly-stochastic Poisson process (also known as Cox process, Cox 1955) has an intensity that depends on its own history as well as the hidden state. We can describe this process by

[TABLE]

where the intensity processes $R^{i}_{t}$ are nonnegative processes that can be computed from the current value of $X_{t}$ and the history of observations $N_{0:s}$ for $s<t$.

To keep the notation simple, we will assume that the vector $R_{t}$ of intensities is given by a function of the hidden state,

[TABLE]

but history-dependence in the form $R_{t}=h(X_{t},N_{0:t^{-}})$ does not significantly increase the difficulty of the filtering (given the observation, any history-dependence of the intensity is deterministic and can be factored out of the conditional expectation of the intensity).

4 The answer to life, the universe and (not quite) everything: the filtering equations

*“I’ve just been created. I’m completely new to the Universe in all respects.

Is there anything you can tell me?”

*— Douglas Adams

The filtering problem is to compute the posterior (or conditional) density of the hidden state conditioned on the whole sequence of observations up to time $t$, $Y_{0:t}$, or equivalently, to compute the posterior expectation (if it exists) of any real-valued measurable function $\phi:\mathbb{R}^{n}\to\mathbb{R}$,

[TABLE]

where we use subscript $\mathbb{P}$ to indicate expectations with respect to the original probability measure $\mathbb{P}$.

That this is not an easy problem should be clear by now, because already the discrete-time filtering task (e.g. in Eq. 13) involved a computation of as many integrals as there are time steps. In continuous time, this would amount to an infinite number of integrals. This continuous-time problem has already been recognized and tackled by mathematicians in the 60s and 70s of the last century, providing formal solutions for the posterior density in terms of stochastic partial differential equations [Kushner, 1962, Zakai, 1969]. In the following, we will derive these equations, using what we have been using in the previous sections as our ultimate “Point of View Gun for nonlinear filtering”:121212The Point of View Gun is a weapon that causes its target to see things from the side of the shooter. Actually, it never appeared in any of Douglas Adams’ novels, but it was featured in the 2005 movie. the change of probability measure method.131313There are other methods to arrive at the same equations, for instance the innovations approach [Bain & Crisan, 2009, Chpt. 3.7] or the more heuristic continuum limit approach originally taken by Kushner [1962].

4.1 Changes of probability measure - once again

Let us once more revisit the change of probability measure in the context of filtering. The goal is to pass from the original probability measure $\mathbb{P}$ (under which the processes behave as our signal and observation model dictates), to an equivalent measure $\mathbb{Q}$, called reference measure, under which the observation process becomes simpler and decouples from the signal process. Here, we will finally generalize our introductory treatment from Section 2.1 to stochastic processes. Unsurprisingly, the calculations are quite similar.

If $\mathbb{P}$ is a probability measure and we have a collection of processes ($X_{t}$ and $Y_{t}$), the measure $\mathbb{P}_{t}$ is the restriction of $\mathbb{P}$ to all events that can be described in terms of the behavior of $X_{s}$ and $Y_{s}$ for $0\leq s\leq t$. If $\mathbb{P}$ and $\mathbb{Q}$ are equivalent, also their restrictions $\mathbb{P}_{t}$ and $\mathbb{Q}_{t}$ are equivalent.141414 For stochastic processes, equivalence implies having the same noise covariance. The Radon-Nikodym theorem [Klebaner, 2005, Theorem 10.6, p. 272ff] then states that a random variable $L_{t}$ exists, such that for all functions $\phi$

[TABLE]

where $L_{t}=\frac{d\mathbb{P}_{t}}{d\mathbb{Q}_{t}}$ is called the Radon-Nikodym derivative or density of $\mathbb{P}_{t}$ with respect to $\mathbb{Q}_{t}$. This is the generalization of Eq. (2) in Section 2.1.

In analogy to Eq. (64), also the conditional expectations can then be rewritten in terms of a reference probability measure $\mathbb{Q}$:

[TABLE]

Equation (65) is known as a Bayes’ formula for stochastic processes (compare Eq. 5) or Kallianpur-Striebel formula. Here, we require a time-dependent normalization $Z_{t}:=\mathbb{E}_{\mathbb{Q}}[L_{t}|Y_{0:t}]$, and $\Braket{\phi_{t}L_{t}}_{\mathbb{Q}}:=\mathbb{E}_{\mathbb{Q}}[\phi_{t}\,L_{t}|Y_{0:t}]$ was introduced to keep the notation concise. This generalizes Eq. (5) above.

But wait: what exactly does the Radon Nikodym derivative $L_{t}$ look like? This really depends on the measure change we are about to perform, but it helps to recall that in a discrete-time (or actually already static) setting the equivalent of the Radon-Nikodym derivative is nothing else than the ratio between two probability distributions. For the filtering problem below, we will choose a reference measure $\mathbb{Q}$ such that the path of the observations $Y_{0:t}$ (or equivalently the set of the increments $dY_{0:t}$) becomes independent of the path of the state process $X_{0:t}$, i.e. $q(X_{0:t},dY_{0:t})=p(X_{0:t})q(dY_{0:t})$. This is very convenient, as this allows us to compute expectations with respect to the statistics of the state process (and we know how to do that). Equation (6) then directly tells us what the likelihood ratio has to look like for this measure change:

[TABLE]

where now $x_{0:t}$ and $dy_{0:t}$ are variables reflecting the whole path of the random variable $X_{t}$ and the set of infinitesimal increments $dY_{0:t}$. Importantly, this particular measure change is agnostic to how the hidden state variable $X_{t}$ evolves in time, but just takes into account how the observations are generated via $p(dy_{t}|x_{t})$.

Let us first consider Gaussian observation noise, as encountered in Section 3.2.1. From Eq. (59), we know that $dY_{t}\sim\mathcal{N}(dY_{t};h(X_{t})dt,\Sigma_{y}dt)$. Further, we choose $q(dY_{t})=\mathcal{N}(dY_{t};0,\Sigma_{y}dt)$ under the reference measure $\mathbb{Q}$. Thus, the Radon-Nikodym derivative $L_{t}=L(X_{0:t},dY_{0:t})$ can be written as

[TABLE]

where in the last step we took the continuum limit $\lim_{dt\to 0}$. Consistently, we would have obtained this result if we had simply (and mindlessly) applied a theorem called Girsanov’s theorem, choosing the reference measure $\mathbb{Q}$ under which the rescaled observations process $\tilde{Y}_{t}$ is a Brownian motion process (Klebaner, 2005, Chapter 10.3, on p. 274, in particular Remark 10.3).

Similarly, we can compute the Radon-Nikodym derivative for observations corrupted by Poisson noise as in Eq. (61). Here, we choose $\mathbb{Q}$ such that $q(dN_{0:t})$ is Poisson with a constant reference rate $\lambda_{0}$. The corresponding density reads

[TABLE]

Again, the same result could have been obtained with a Girsanov theorem [see Bremaud, 1981, Chapter VI, Theorems T2 and T3].

4.2 Filtering equations for observations corrupted by Gaussian noise

We are now equipped with the necessary tools to tackle the derivation of the filtering equations. Here, the derivation will be briefly outlined [for a more detailed and formal derivation, see Van Handel, 2007, Bain & Crisan, 2009].

As we stated in the beginning of this Section, we want to find a convenient reference measure which decouples the signal and observations and at the same time turns the observations into something simple. Recall that the Radon-Nikodym derivative (expressed for the rescaled observations process $\tilde{Y}_{t}$) then takes the form

[TABLE]

which evolves according to the following SDE (see A.1 for calculation steps):

[TABLE]

Under $\mathbb{Q}$, the stochastic differential can be taken inside the expectation [see Van Handel, 2007, Chapter 7, Lemma 7.2.7], and we therefore obtain using Itô’s lemma151515Recall that this corresponds to a Taylor expansion up to second order (i.e. first order in $dt$, since $\mathcal{O}(dW_{t})=dt^{1/2}$) for diffusion processes, which is why we have to consider the product of differentials (product rule for stochastic differentials).

[TABLE]

where introduced the short-hand notation $\mathbb{E}_{\mathbb{Q}}\left[\cdot|Y_{0:t}\right]=\Braket{\cdot}_{\mathbb{Q}}$ for the conditional expectation. We recall from Section 3 that for both of the signal models, we may write the time evolution of $\phi_{t}=\phi(X_{t})$ as

[TABLE]

where we denote $\mathcal{A}\phi_{t}=\mathcal{A}\phi(X_{t})$. $M^{\phi}_{t}$ is a martingale that is independent of the observations under $\mathbb{Q}$, and thus $\Braket{dM_{t}^{\phi}}_{\mathbb{Q}}=0$ as well as $\Braket{L_{t}\,dM_{t}^{\phi}}_{\mathbb{Q}}=0$. Therefore, we only retain the $dt$ term under the conditional expectation. The first term in Eq. (72) can then be computed using Eq. (71):

[TABLE]

which is the SDE of the unnormalized posterior expectation. Here we further used that $\Braket{(d\phi_{t})\,(dL_{t})}=0$, because the noise in state and observations are independent.

Note that the evolution equation of the normalization constant $Z_{t}$ in Eq. (65), $dZ_{t}=d\Braket{L_{t}}_{\mathbb{Q}}$, corresponds to Eq. (73) with the constant function $\phi=1$. The differential $dZ_{t}^{-1}$ is obtained by consecutive application of Itô’s lemma (Eq. 47). By plugging Eq. (73) and $dZ_{t}^{-1}$ into Eq. (72) and rewriting everything in terms of expectations under $\mathbb{P}$ using Eq. (65), one finally obtains the evolution equation for the posterior expectation, the so-called Kushner-Stratonovich equation (KSE, Bain & Crisan, 2009, p. 68, Theorem 3.30, cf. A.2 for calculation steps):

[TABLE]

Equivalently, it can straightforwardly be expressed in terms of the original observation process in Eq. (59):

[TABLE]

In analogy to the calculations in Section 3, one may also pass from an evolution equation for the expectations to an adjoint equation for the conditional probability density,

[TABLE]

Writing Eq. (74) and (76) in terms of the (adjoint of the) infinitesimal generator of the signal process, allows us to use any signal process for which $\mathcal{A}$ is known. For instance, if the signal process is a Markov chain on a finite set $S$, the expression $p(x|Y_{0:t})$ can be interpreted as the vector of posterior probabilities $\hat{p}(t)$, with entries $\hat{p}_{i}(t)$ denoting the probability to be in state $i$ given all observations up to time $t$. The generator $\mathcal{A}^{{\dagger}}$ is then represented by the matrix $A^{\top}$ that has appeared in the evolution equation for the prior density, Eq. (41). Specifically, $\hat{p}_{i}(t)$, evolves as

[TABLE]

where $h_{i}=h(i)\in\mathbb{R}^{l}$, $i=1,...,n$ and $h$ is an $l\times n$-matrix whose columns are the $h_{i}$’s. Eq. (77) is known as the Wonham filter [Wonham, 1964], and it is a finite-dimensional SDE that completely solves the filtering problem.

Equation (76) is a stochastic integro-differential equation, known as Kushner-Stratonovich equation (KSE) [Stratonovich, 1960, Kushner, 1962], and its solution is in general infinite-dimensional. This fundamental problem is easily illustrated by considering the time evolution of the first moment, using $\phi(x)=x$:

[TABLE]

For non-trivial (i.e. non-constant) observation functions $h$, any moment equation will depend on higher-order moments due to the posterior covariance between the observation function $h$ and the function $\phi$, which effectively amounts to a closure problem when $f(x)$ is nonlinear. This is not surprising; even the Fokker-Planck equation (58) (the evolution equation for the prior distribution) presents such a closure problem. In some cases (e.g. when using kernel or Galerkin methods), it is more convenient to use the evolution equation of the unnormalized posterior density $\varrho(X_{t}|Y_{0:t})$:

[TABLE]

which is a linear stochastic partial differential equation (SPDE), the Zakai equation [named after Zakai, 1969].

In very rare cases, under specific signal and observation models, the moment equations close, e.g. in the Kalman-Bucy filter [Kalman & Bucy, 1961, see Section 4.3 below] or the Beneš filter [Benes, 1981]. Other finite-dimensional filters include the Daum filter [Daum, 1986] for continuous-time processes and discrete-time measurements. However, in most cases that occur in practice, the KSE needs to be approximated using a finite-dimensional realization. For instance, one could use the KSE as a starting point for these approximations, e.g. Markov-chain approximation methods [Kushner & Dupuis, 2001] or projection onto a finite-dimensional manifold [Brigo et al., 1998, 1999], which can be shown to be equivalent to assumed density filtering (ADF), or Galerkin-type methods with specific metrics and manifolds [see Armstrong & Brigo, 2013]. Other numerical algorithms associated with overcoming the numerical burden of solving the Kushner or Zakai equation rely on a Fourier approximation of the involved densities, and the fact that convolutions correspond to simple products in Fourier space [Mikulevicius & Rozovskii, 2000, Brunn et al., 2006, Jia & Xin, 2010].

4.3 A closed-form solution for a linear model: Kalman-Bucy filter

In models where the hidden drift function $f(X)$ and the observation function $h(X)$ are linear, i.e. $f(x)=Ax$ and $h(x)=Bx$, and the initial distribution is Gaussian, there exists a closed-form solution to the filtering problem. In this case, the KSE (Eq. 74) closes with the second posterior moment $\Sigma_{t}$, i.e. the evolution equation for $\Sigma_{t}$ becomes independent of the observations process, and the posterior density corresponds to a Gaussian with time-varying mean $\mu_{t}$ and variance $\Sigma_{t}$. The dynamics of these parameters are given by the Kalman-Bucy filter (KBF, Kalman & Bucy, 1961) and form a set of coupled SDEs:

[TABLE]

The posterior variance follows a differential Riccati equation and, due to its independence from observations as well as from the posterior mean, it can be solved offline.

4.4 Filtering equations for observations corrupted by Poisson noise

In analogy to the previous section, the formal solution to the filtering problem with observations corrupted by Poisson noise (Eq. 61) can also be derived with the change of probability measure method. We will very briefly outline the derivation, referring to similarities to continuous-time derivations.161616A very detailed derivation is offered in Bremaud [1981, p. 170ff.] or, more intuitively, in Bobrowski et al. [2009, SI] and in Surace [2015, p. 41ff].

The idea is again to make use of the Kallianpur-Striebel formula (Eq. 65) and to rewrite the posterior expectation under the original measure $\mathbb{E}_{\mathbb{P}}[\phi_{t}|N_{0:t}]$ in terms of an expectation under a reference measure $\mathbb{Q}$, under which hidden process $X_{t}$ and observation process $N_{t}$ are decoupled. Using a Girsanov-type theorem for point processes [see Bremaud, 1981, Chapter VI, Theorems T2 and T3], the measure is changed to the reference measure $\mathbb{Q}$ under which all point processes have a constant rate $\lambda_{0}$. The Radon-Nikodym derivative reads

[TABLE]

which solves the SDE

[TABLE]

We can now repeat the calculations of the previous section. First, we obtain

[TABLE]

Here, we used again that under $\mathbb{Q}$, differentiation and expectation can be interchanged.

Using $\phi_{t}=1$ gives us the evolution of the time-dependent normalization $Z_{t}$ in the Kallianpur-Striebel formula (65). We can use these, together with Itô’s lemma (Eq. 47) to compute $dZ_{t}^{-1}=d\Braket{L_{t}}_{\mathbb{Q}}^{-1}$, to obtain a point-process observations analogue to the KSE for the normalized posterior estimate:171717See Appendix A.3 for detailed derivation steps.

[TABLE]

where $\text{diag}(x)$ denotes a diagonal matrix with diagonal entries given by the vector $x$. The adjoint form of this equation, i.e. the evolution equation for the posterior density $p(x|N_{0:t})$, reads:

[TABLE]

Note the structural similarity to the Kushner equation (Eq. 76): it also relies on a Fokker-Planck term denoting the prediction, and a correction term that is proportional to the posterior density, the ‘innovation’ $dN_{t}-\Braket{h(X_{t})}dt$, as well as a local correction $h(x)-\Braket{h(X_{t})}$. The difference is that the observation noise covariance $\Sigma_{y}$ in the Kushner equation has been replaced by a diagonal matrix whose components are proportional to the rate function in each observable dimension. Considering that the observations are Poisson processes, this is not surprising: for Poisson processes, the variance is proportional to the instantaneous rate, and thus, analogously, the correction term in this equation has a similar proportionality.

Similarly, we find for the unnormalized posterior density $\varrho(x|N_{0:t})$:

[TABLE]

Analogously to Eqs. (76) and (79), these equations are obtained by integrating the equation for the unnormalized posterior estimate (Eq. 85) and the normalized posterior estimate (Eq. 86) twice.

4.5 Down to Earth - an example from decision making

In this section, we have seen that the filtering equations can be derived by changing from the original measure to a measure under which the signal and observation processes are independent. Interestingly, what we have seen here also gives us a recipe of how to treat a filtering problem in general: all we need is a characterization of the state dynamics in terms of the infinitesimal generator $\mathcal{A}$. Further, any information about the observations is carried by the Radon-Nikodym derivative.

To illustrate this, let us consider the following example from dynamical decision making. In order to make a decision that relies on an unobserved, but dynamic, variable, animals have to be able to integrate the (noisy) observations about this variable in a way that discounts older evidence in favor of more recent evidence, depending on the dynamics of the variable. In other words: they have to do filtering. In Piet et al. [2018] an experimental paradigm is considered where the hidden variable can be in one of two states, i.e. $X_{t}\in\{S_{1},S_{2}\}$, which switches states with a hazard rate $a$ and which influences the click rates $r_{+}$ and $r_{-}$ of two speakers. More precisely, whenever the system is in state $S_{1}$, speaker 1 will have click rate $r_{+}$ and speaker 2 will have click rate $r_{-}$, and vice versa if the system is in state $S_{2}$ (Fig. 4a). Presented with clicks from these speakers, rats have to make a decision about the state the environment is in at the end of the trial. The optimal solution to the problem is then the filtering distribution conditioned on the clicks at the end of the trial, and the optimal decision corresponds to the state with the highest posterior probability. This study has shown that rats are able to perform optimal evidence discounting.181818We are glad they didn’t use mice, as these animals, according to the Hitchhiker’s guide, are the most intelligent species on planet Earth and as such would surely have outperformed the optimal solution. Rats are close enough, though.

To come up with the solution for the filtering task, Piet et al. [2018] consider the evolution of the log-likelihood as the decision variable (based on Veliz-Cuba et al., 2016) and derive this by taking the continuum-limit of the discrete-time process. In principle, this approach is perfectly equivalent to a change of measure in the log domain. Here, we will re-derive their result for the optimal solution of the filtering task by directly applying our ‘recipe’.

Let us first consider the state process, which is a discrete-state Markov model. Without observations, the probabilities of the hidden states evolve as $d\tilde{p}_{t}=\mathcal{A}^{\dagger}\tilde{p}_{t}\,dt$, where $\tilde{p}_{t}^{i}=p(X_{t}=S_{i})$. It is easy to check that for this model the adjoint generator matrix is given by

[TABLE]

where $a$ denotes the switching rate, i.e. $a\,dt$ is the probability of switching the state.

The observation model is a two-dimensional point process, and $N_{t}^{j}$ denotes the number of clicks in speaker $j$ up to time $t$. Let $h_{ji}=h_{j}(X=S_{i})$ be the rate of speaker $j$ if the hidden variable is in state $S_{i}$. Thus, the evolution of the posterior probability $p_{t}$ is given by (cf. Eq. 87):

[TABLE]

where $\langle h_{j}\rangle=\sum_{i}h_{ji}p_{t}^{i}$. Since this particular system is 2-dimensional, i.e. $p_{t}^{2}=1-p_{t}^{1}$, we can substitute this in the expression for the first component of the posterior and get a one-dimensional equation for $p_{t}^{1}$.

[TABLE]

where $\langle h_{j}\rangle=h_{j1}p_{t}^{1}+h_{j2}(1-p_{t}^{1})$. At the end of the trial, $S_{1}$ is the optimal choice whenever $p_{t}^{1}>1/2$.

For several reasons, it might be desirable to define the decision variable as the log likelihood ratio of being in state $S_{1}$ as opposed to being in $S_{2}$. Let $\alpha_{t}=\log\frac{p_{t}^{1}}{p_{t}^{2}}=\log\frac{p_{t}^{1}}{1-p_{t}^{1}}$. In order to derive its evolution equation, we can directly apply Itô’s lemma for point processes to Eq. (91) with $\phi(x)=\frac{x}{1-x}$ and after straightforward, but tedious, algebra arrive at the desired SDE.

[TABLE]

Note that if $h_{11}-h_{12}=h_{22}-h_{21}=r_{+}-r_{-}$ (which is the experimental setting in Piet et al., 2018), this equation becomes very simple:

[TABLE]

and resembles Eq. 9 in Piet et al. [2018]. Without this symmetry, there is a drift term modulated by the difference in click rates, indicating that the absence of clicks is informative for estimating the current state. In Figure 4, we plotted the log likelihood ratio, both computed from the posterior probabilities as in Eq. (90), and directly from running Eq. (93). Unsurprisingly, both plots lie exactly on top of each other.

Note that this result was obtained by simply plugging in a model, i.e. the signal and the observation process, into the solution to the filtering equation, and making use of the fact that for a finite state space, the Kushner equation for point processes, Eq. (87), becomes finite-dimensional. Unlike in Piet et al. [2018] or in Veliz-Cuba et al. [2016], we did not have to explicitely carry out the continuum limit - in fact, this is implicitely taken care of by using the appropriate Radon-Nikodym derivative for this observation model. This allows for much more flexibility when the models and/or experimental settings become more complicated, for instance if we want to increase the number of states, modify the state dynamics or modify the properties of the speakers.

A fully annotated code for this example is available in our github repository [Kutschireiter, 2019].

5 Don’t panic: Approximate closed-form solutions

*“It is a mistake to think you can solve any major problems just with potatoes.”

*— Douglas Adams

If the signal model is a jump-diffusion, the KSE (Eq. 74, Eq. 86) is infinite-dimensional, a fact known as ‘closure problem’. In other words, except for some important exceptions for very specific models, such as the KBF or the Beneš filter [Benes, 1981], solutions to the general filtering problems are not analytically accessible. Furthermore, unlike for observations following a diffusion process, no closed-form filter for point-process observations is known. However, there exist important approximate closed-form solutions, which address the closure problem by approximating the posterior density in terms of a set number of sufficient statistics.

Here, we will briefly outline some important examples: first, the Extended Kalman-Bucy Filter and related methods for point-process observations that rely on a series expansion of the functions in the generative model, such that the posterior is approximated by a Gaussian density. We will further describe Assumed Density Filters, that choose a specific form of the posterior and propagate the KSE according to this approximation.

5.1 The Extended Kalman Bucy Filter and related approaches

Based on the Kalman-Bucy filter (Section 4.3), the extended Kalman-Bucy filter (EKBF) is an approximation scheme for nonlinear generative models of the form

[TABLE]

The EKBF approximates the posterior by a Gaussian with mean $\mu_{t}$ and variance $\Sigma_{t}$, whose dynamics are derived by local linearization (around the mean) of the nonlinearities in the model [Jazwinski, 1970, p. 338, Example 9.1]:

[TABLE]

where $F_{ij}=\frac{\partial f_{i}}{\partial x_{j}}$ and $H_{ij}=\frac{\partial h_{i}}{\partial x_{j}}$ denote the Jacobian of the hidden drift function and the observation function, respectively. For models with multimodal posteriors, this approximation often breaks down: e.g. if the noise covariance $\Sigma_{y}$ is large, the mean of the EKBF tends to ‘get stuck’ in one of the modes.

Similar approximations exist for point-process observations. One way to achieve this would be to simply construct an EKBF by assuming Gaussian noise in the observations, together with the appropriate linearization [see paragraph below Eq. 17 in Eden, 2007]. Another way that allows the point-process observations to directly enter the expressions for mean and variance relies on a Taylor expansion in the log domain of the approximated posterior up to second order (see Eden et al., 2004 for discrete-time and Eden & Brown, 2008 for the continuous-time models). The continuous-time approximate filter for point processes can be seen as the point-process analogue of the EKBF [cf. Eden & Brown, 2008, extended by nonlinearity in the hidden state process]:

[TABLE]

with

[TABLE]

5.2 Assumed density filtering

The idea of assumed density filters (ADF) is to specify a set of sufficient statistics, which is supposed to approximate the posterior density, derive evolution equations from the KSE, i.e. from Equations (74) and (86), and approximate expectations within these evolution equations under the initial assumptions. To be less abstract, consider approximating the posterior density by a Gaussian density. Then it suffices to derive evolution equations for mean $\boldsymbol{\mu}_{t}$ and variance $\Sigma_{t}$ of the approximated Gaussian posterior. In these evolution equations, higher-order moments will enter, which in turn can be expressed in terms of mean and variance for a Gaussian.

As a concrete example, let us consider a Gaussian ADF for point-process observations (the treatment for diffusion-observations is completely analogous). Consider the SDEs for the first two moments of the posterior (cf. Eq. 86, detailed derivation in Appendix A.4):

[TABLE]

The effective realization of the ADF will crucially depend on the specifics of the signal model defined by Eq. (43) and the observation model (61), respectively. For example, Pfister et al. [2009] consider an exponential rate function $h(x)\propto\exp(\beta x)$, which leads to a variance update term that is independent of the spiking process. Of particular interest for decoding tasks in neuroscience are ADFs with Gaussian-shaped rate function (e.g. Harel et al., 2018).

However, for some models ADFs cannot be computed in closed form. Consider for simple example a rate function that is a soft rectification of the hidden process, e.g. $h(x)=\log(\exp(x)+1)$, which, when taking expectations with respect to a Gaussian, does not admit a closed-form expression in terms of mean $\mu_{t}$ and variance $\Sigma_{t}$.

6 Approximations without Infinite Improbability Drive: Particle Methods

*“Each particle of the computer, each speck of dust

held within itself, faintly and weakly, the pattern of the whole.”

*— Douglas Adams

Particle filtering (PF) is a numerical technique to approximate solutions to the filtering problem by a finite number of samples, or ‘particles’, from the posterior. Thus, they serve as a finite-dimensional approximation of the KSE, overcoming the closure problem. The true posterior is approximated by the empirical distribution formed by the particle states $X_{t}^{(i)}$, i.e. a sum of Dirac-delta functions at the particle positions $\delta(x-X_{t}^{(i)})$, and, if it is a weighted PF, weighted by their corresponding importance weights $w_{t}^{(i)}$,

[TABLE]

with $\sum_{i}w_{t}^{(i)}=1$ ensuring normalization. Consequently,

[TABLE]

The rationale is based on a similar idea as using the Euler-Maruyama scheme to numerically solve the Fokker-Planck equation and its associated equation for the posterior moments. As a numerical recipe [for instance provided by Doucet et al., 2000, Doucet & Johansen, 2009, for discrete-time models], it is easily accessible, because in principle no knowledge of the Fokker-Planck equation, nonlinear filtering theory or numerical methods for solving partial differential equations is needed.

In this section, weighted particle filters will be introduced from a continuous-time perspective based on the change of probability measure formalism [roughly following Bain & Crisan, 2009, Chapt. 9.1, and extending this to point-process observations]. From this formalism, we derive dynamics for the weights and link these to the ‘curse of dimensionality’. Finally, to give context for readers more familiar with discrete-time particle filtering, the continuous-time perspective will be linked to the ‘standard’ particle filter (PF).

6.1 Particle filtering in continuous time

Based on sequential importance sampling, both samples (or ‘particles’) $X_{t}^{(i)}$ as well as their respective weights are propagated through time. As we have seen before in Section 2.1.1, importance sampling amounts to a change of measure from the original measure $\mathbb{P}$ to a reference measure $\mathbb{Q}$, from which sampling is feasible. Here, the idea is to change to a measure under which the observation processes are independent of the hidden process, effectively enabling us to sample from the hidden process. In other words, the particle positions evolve as specified by the infintesimal generator $\mathcal{A}$ of the hidden process (e.g. Eq. 43 if the hidden process is a jump-diffusion process).

Why this should be the case is rather intuitive when recalling the Kallianpur-Striebel formula (65):

[TABLE]

If we want to approximate the left-hand side of this equation with empirical samples, it would require us to have access to samples from the real posterior distribution, which is usually not the case. However, since under the measure $\mathbb{Q}$ on the right-hand side the hidden state and observations are decoupled, the estimate is approximated by empirical samples that correspond to realizations of the hidden process:

[TABLE]

$\bar{Z}_{t}=\sum_{i=1}^{M}L_{t}(X_{t}^{(i)})$ is an empirical estimate of the normalization constant.

Thus, we just need to evaluate the Radon Nikodym derivative at the particle states $X_{t}^{(i)}$, giving us the importance weight $w_{t}^{(i)}$ of particle $i$ at time $t$. For observation corrupted by Gaussian noise (cf. Eq. 70), this reads:

[TABLE]

Similiarly, we find for point-process observations (cf. Eq. 82)

[TABLE]

6.1.1 Weight dynamics in continuous time

If one is interested how the weight of particle $i$ changes over time, it is possible to derive an evolution equation for the particle weights. Using Itô’s lemma, we find:

[TABLE]

For continuous-time observations, (cf. Eq. 71) yields

[TABLE]

where $\bar{h}_{t}:=\sum_{i}w_{t}^{(i)}h(X_{t}^{(i)})=\bar{Z}_{t}^{-1}\sum_{i}L_{t}^{(i)}h(X_{t}^{(i)})$ is the weighted estimate of the observation function $h_{t}$ (i.e. under the original measure $\mathbb{P}$). Applying Itô’s lemma on Eq. (109) to obtain $d\bar{Z}_{t}^{-1}$, we find for the dynamics of the weights

[TABLE]

Similarly, with Eq. (83) we find for point-process observations:

[TABLE]

and thus, using Itô’s lemma for point processes to obtain $d\bar{Z}_{t}^{-1}$, with Eq. (107):

[TABLE]

Interestingly, there is a striking similarity to the dynamics of the importance weights and the Kushner equation (76) and the point-process observation equivalent of the Kushner equation (Eq. 87), respectively. The weight dynamics seem to directly correspond to the dynamics of the correction step (with true posterior estimates replaced by their empirical counterparts). This is rather intuitive: since we chose a change of measure under which the particles follow the hidden dynamics, serving as the prediction step, the observation dynamics have to be fully accounted for by the weight dynamics, in a way to be consistent with the Kushner equation.

It is important to note that the weight dynamics inevitably lead to a system, in which all but one weight equals zero, the so-called weight degeneracy. In this degenerate state, the particle system cannot represent the posterior sufficiently, and hence has to be avoided numerically, e.g. by resampling the particles from the weight distribution and resetting the weights to $1/M$. The time scale on which this weight degeneracy happens depends on the number of observable dimensions, in other words it is accelerated as the dimensionality of the system is increased [Surace et al., 2019b]. This is a form of the so-called ’curse of dimensionality’, a common nuisance in weighted particle filters.

6.1.2 Equivalence between continuous-time particle filtering and bootstrap particle filter

The practitioner who is using PF algorithms in their numerical implementations might usually be more familiar with the discrete-time formulation. Further, since the continuous-time formulation of the particle filter based on the measure change formalism seems to be so different from the discrete-time formulation, one might rightfully ask how these two concepts are related and whether they are equivalent. Indeed, we will now quickly show that the continuous-time PF in Section 6.1 corresponds to the Bootstrap PF in a continuous-time limit. More precisely, if we apply the Bootstrap PF to a time-discretized version of our hidden state process model and observation model, and then take the continuum limit, we will regain the equations for the weights as in Section 6.1.

Irrespectively of the generator $\mathcal{A}$ of the hidden process, it is straightforward to write the hidden process in terms of a transition density, with $t-dt$ corresponding to the previous time step. This acts as the proposal density $\pi(X_{t}|X_{0:t-dt}^{(i)},Y_{0:t})=p(X_{t}|X_{t-dt}^{(i)})$. Consider for example a drift-diffusion process (Eq. 43 with $J=0$). Then the particles are sampled from the time-discretized transition density, $X_{t}^{(i)}\sim p(X_{t}|X_{t-dt}^{(i)})$, which is given by:191919Remark: The so-called Euler-Maruyama scheme for numerical implementation of diffusion processes is based on the very same discretization.

[TABLE]

For observations corrupted by Gaussian noise, the emission likelihood is given by the emission probability for the instantaneous increments $dY_{t}$ in Eq. (43), i.e.

[TABLE]

such that

[TABLE]

It is evident that the proposal of continuous-time particle filtering and that of the BPF match, and it remains to show that the same holds for the importance weights. In other words, when taking the continuum limit of Eq. (117), we should be able to recover Eq. (105). Keeping only terms up to $\mathcal{O}(dt)$, we find

[TABLE]

where the term $\propto(dY_{t})^{2}$ was absorbed in the normalization because it is independent of the particle positions. Thus, the continuous-time limit $dt\to 0$ of Eq. (35) reads

[TABLE]

which, up to the normalization constant $\bar{Z}_{t}$, is equivalent to Eq. (105).

For point-process observations, the emission likelihood is given by $p(dN_{t}|X_{t})$, which is defined by the Poisson density in Eq. (61). Neglecting the term that is independent of the particle positions (which is absorbed in the normalization), it can be rewritten as:

[TABLE]

Again, since $\tilde{w}_{t}^{(i)}\propto\prod_{s}p(dN_{s}|X_{s}^{(i)})$, the continuous-time limit of the unnormalized importance weight is

[TABLE]

The explicit dependence of Eq. (106) on the reference rate $\lambda_{0}$ can be absorbed in the normalization constant, yielding equivalent expressions for the normalized weights $w_{t}^{(i)}$.

6.2 The Feedback Particle filter

In contrast to weighted particle filtering approaches, unweighted approaches for particle filtering exist, for example the Ensemble Kalman (Bucy) Filter [Evensen, 1994, Bergemann & Reich, 2012], the Feedback Particle Filter (FBPF, Yang et al., 2013, 2014), the (stochastic) particle flow filter [Daum et al., 2010, de Melo et al., 2015] or the point-process analogue to the FBPF [Surace et al., 2019a]. Since these methods do not rely on importance weights in the first place, there is no weight degeneracy. Unweighted particle filters therefore hold the promise of avoiding the curse of dimensionality [see Surace et al., 2019b].

All of these methods have in common that the posterior is approximated by equally weighted particles, i.e.:

[TABLE]

Consequently, the observations have to enter the particle dynamics in such a way that the particles are moved towards regions of high posterior density. As an example, we will outline how this is achieved in the Feedback particle filter.

In the FBPF, the observations directly enter the particle dynamics, which evolve according to the Itô SDE:

[TABLE]

where $B_{t}^{(i)}$ are uncorrelated vector Brownian motion processes, $K(X_{t}^{(i)},t)$ is the gain matrix, and $\bar{h}_{t}=\frac{1}{N}\sum_{i=1}^{N}h(X_{t}^{(i)})$ denotes the particle estimate of the observation function. The components of the additional vector-valued drift function $\Omega(X_{t}^{(i)},t)$ are given by Yang et al. [2016]

[TABLE]

The gain $K$ is the solution of a boundary value problem that emerges from an optimal control problem [Yang et al., 2016, Eqs. 4, 5]. It is chosen such that it minimizes the Kullback-Leibler divergence between the particle distribution and the posterior filtering distribution (conditioned on proper initialization), which leads to the following conditions:

[TABLE]

with

[TABLE]

In general, the gain matrix $K(x,t)$ cannot be solved for in closed form, and in practical implementations relies on a numerical solution of the Euler-Lagrange boundary value problem [Taghvaei & Mehta, 2016]. For instance, one way to approximate the gain $K(x)$ is using a Galerkin approximation. In particular, choosing the coordinate functions as basis functions, the so-called constant gain approximation reads [Yang et al., 2016, Eq. 20]:

[TABLE]

In this approximation, the gain is constant with respect to the particle positions, i.e. each particle has the same gain, but still changes as a function of time. In this approximation, the additional drift function $\Omega$ in Eq. (123) is zero.

For a linear state space model the FBPF with constant-gain approximation becomes exact and is identical202020up to small numerical difference when computing the gain to the ensemble Kalman-Bucy filter (EnKBF, Bergemann & Reich, 2012, Taghvaei et al., 2018), which can be shown to be asymptotically exact [Künsch, 2013]. More precisely, for a linear model with $f(x)=Ax$, $G(x)=\Sigma_{x}^{1/2}$ and $h(x)=Bx$, the gain $K(x)$ can be solved for in closed form, using the knowledge that the posterior is Gaussian at all times, and is given by $K=\hat{\Sigma}_{t}B$. The particles are thus propagated according to (cf. Eq 123)

[TABLE]

where $\hat{\mu}_{t}$ and $\hat{\Sigma}_{t}$ denote the posterior mean and variance as estimated from the particle positions.

6.3 Particle filters in action

Here, we have seen how weighted particle filters can be constructed for nonlinear, continuous time-filtering. Further, we introduced unweighted particle filters as an alternative, which, for some systems (for instance with a high-dimensional observation model), may be advantageous over standard particle filtering. However, unweighted particle filters come at the cost of having to compute the gain function, which can be numerically expensive. Here, we want to demonstrate with an example how these algorithms can be applied as numerical algorithms to a simple nonlinear filtering problem.

Let us consider a hidden process with drift function $f(x)=-4x(x^{2}-1)$ and diffusion constant $g(x)=\Sigma_{x}^{-1/2}$. The corresponding stationary probability density of this nonlinear model is a bimodal distribution, with peaks at $x=\pm 1$. First, we use observations corrupted by Gaussian noise, e.g. with $h(x)=x$ (note that the model is still nonlinear due to the nonlinearity in the state transition). This filtering problem cannot be solved for in closed form and thus we have to use a finite-dimensional approximation, such as the particle filter,

With Eqs. (115) and (117), the particle transition and weight dynamics for a standard particle filter (Bootstrap particle filter, BPF) for this model is given by

[TABLE]

After each iterative step, the weights need to be normalized according to $w_{t}^{(i)}=\tilde{w}_{t}^{(i)}/\sum_{j}\tilde{w}_{t}^{(j)}$. An example tracking simulation is shown in Fig. 5a and b. For comparison, we also consider an extended Kalman-Bucy filter (see Section 5.1) and a feedback particle filter with constant gain approximation (Fig. 5a,c).

Similarly, we might also consider point-process observations with intensity $g(X_{t})$. For example, let us use a Gaussian-shaped rate function $g(x)=g_{0}\exp(\frac{x-m_{o}}{2s_{0}^{2}})$ for two sensors with peaks at $m_{0}=\pm 1$ and width $s_{0}$. Figure 5d shows that a particle filter is able to track the hidden state reasonably well based only on only the events elicited from these two sensors. Note also the similarity of this model to the two-state HMM model example we used earlier in Section 4.5.

A fully annotated code for this example, which also explains the technical details of this simulation, is available in our github repository [Kutschireiter, 2019].

7 The restaurant at the end of the universe: Take-away messages

*“For where he had expected to find nothing,

there was instead a continuous stream of data.”

*— Douglas Adams

In this tutorial, we reviewed some of the theory of nonlinear filtering theory. In fact, we wanted to emphasize that the change of probability measure approach can be used as a universal tool in this context. Not only does it help to derive the corresponding equations such as the filtering equations, but also leads to a deeper understanding of particle filters.

Let us once more summarize the main take-away messages of this tutorial:

The change of measure method comes in handy whenever expectations are easier to evaluate under a different measure. For example, when computing conditional expectations, it is often easier to compute under a measure in which the random variables are independent. The Radon Nikodym derivative acts as the ‘conversion factor’ between the expectations under the different measures. 2. 2.

The filtering equations can be derived by changing to a reference measure under which signal process and observation process are independent. Since this measure change acts on the observation process and leaves the signal process untouched, the filtering equations have the same structure independently of the signal process (which enters in terms of its infinitesimal generator $\mathcal{A}$). Further, all the information about the observations is carried by the Radon-Nikodym derivative. 3. 3.

For a general continuous-time and continous state space nonlinear filtering problem the filtering equations suffer from a closure problem. Suitable approximations are based on a finite-dimensional representation of the filtering density, e.g. in terms of a finite number of statistics (such as the EKF or ADFs) or a finite number of samples (such as PFs). 4. 4.

Bootstrap particle filtering can be derived by again changing to a reference measure under which the signal and observation processes are decoupled (which is the very same that was used for deriving the filtering equations), and evaluating the expectations empirically. The importance weights correspond to the Radon-Nikodym derivative of this measure change, evaluated at the particle positions. Thus, the contribution of the signal process (to the solution of the filtering problem) enters via the positions of the particles (‘prediction’), while the contribution of the observations enters via the importance weights (‘update’).

So long, and thanks for all the fish: Acknowledgments

The authors would like to express their gratitude towards Christian Horvat, Jannes Jegminat, Luke Rast, Hazem Toutounji, and the anonymous reviewers for taking a test ride and providing constructive feedback. Further, we thank Johannes Bill and Jan Drugowitsch for helpful discussions.

8 Funding

AK, SCS and JPP were supported by the Swiss National Science Foundation grants PP00P3_179060 and PP00P3_150637 (AK, SCS, JPP) and P2ZHP2_184213 (AK).

Appendix A Mostly harmless: detailed derivation steps

*“You ask this of me who have contemplated the very vectors of the atoms of the Big Bang itself?

Molest me not with this pocket calculator stuff.”

*— Douglas Adams

Here, we provide some additional steps that we used in our derivations and left out to keep the main text concise.

A.1 Evolution equation of Radon-Nikodym derivative $L_{t}$ (Eq. 71)

In Section 4.2, we took the form of the Radon-Nikodym derivative in Eq. (70), and from this form deduced the evolution equation in Eq. (71), i.e., Eq. (70) is the solution of Eq. (71). Showing this is a nice application of Itô’s lemma and is therefore outlined here in more detail.

Define

[TABLE]

Then, with Itô’s lemma (recall that for diffusion processes this amounts to a Taylor expansion up to second order in the differential):

[TABLE]

A.2 Kushner-Stratonovich Equation (Eq. 74)

The Kushner-Stratonovich equation describes the SDE for the posterior estimate under the original measure $\mathbb{P}$. With the Kallianpur-Striebel formula 65, we first write the posterior estimate under the reference measure $\mathbb{Q}$. This is given by Eq. (72), revisited here for convenience:

[TABLE]

Let us evaluate these terms separately. The first term can be expressed as

[TABLE]

where we first used the Itô lemma for products, and then made use of $\Braket{dW_{t}=0}_{\mathbb{Q}}$. The term $d\phi_{t}dL_{t}$ equals zero because their noise components are independent. Further, we used the Kallianpur-Striebel formula to rewrite the expressions in terms of expectations under $\mathbb{P}$. This equation is basically the unnormalized measure in Eq. (73) multiplied by $Z_{t}^{-1}$, which we could have used directly.

To obtain $dZ_{t}^{-1}$, we need to apply Itô’s lemma to the SDE of the normalization constant $Z_{t}$ in Eq. (65). $dZ_{t}=d\Braket{L_{t}}_{\mathbb{Q}}$ is given by to Eq. (73) with $\phi=1$.

[TABLE]

Using Itô’s lemma, we find:

[TABLE]

Thus, the second term in Eq. (72) reads:

[TABLE]

Finally, the third term uses the SDE of the unnormalized posterior expectation in Eq. (73) and the SDE in Eq. (138), keeping only terms up to $\mathcal{O}(dt)$.

[TABLE]

Adding up and rearranging the terms, we end up with

[TABLE]

A.3 Kushner-Stratonovich equation for point-process observations (Eq. 86)

The steps are analogous to those taken in the previous section, with a little caveat: here, SDEs are governed by a point process due to the observation process $N_{t}$, so whenever we apply Itô’s lemma, we need to consider an infinite-dimensional Taylor expansion in the differential, since $dN_{t}^{n}=dN$. Also, for simplicity, the following derivation is done for a one-dimensional observation process. However, it is straightforward to be generalized to $l$ dimensions by considering that the observations in each dimension are independent conditioned on all observations up to $t$, which leading to the product in Eq. (82) and the sums in Eq. (83ff).

First, using the Kallianpur Striebel formula, we compute the SDE for the normalized posterior expectation by expressing it in terms of the unnormalized posterior SDE:

[TABLE]

where we used Itô’s lemma for products and the fact that under $\mathbb{Q}$, we can interchange expectation and differentiation.

Again, we compute the terms separately. Using the evolution equation of the unnormalized measure (Eq. 85, we find:

[TABLE]

For the second term, we again write down the SDE for the normalization constant $Z_{t}$ and its inverse. From Eq. (85) with $\phi_{t}=1$ we find:

[TABLE]

The SDE for its inverse is obtained by using Itô’s lemma for point processes (Eq. 47):

[TABLE]

Now we can write

[TABLE]

Finally, for the last term, we only keep terms of $\mathcal{O}(dN)$:

[TABLE]

Adding up and rearranging the terms, we end up with

[TABLE]

A generalization to multivariate $N_{t}$ and $h_{t}$, respectively, is done by treating all observable dimensions separately and summing up.

A.4 ADF for point-process observations (Eq. 99 and 100)

For the ADF in subsection 5.2, we need the SDEs for the first two posterior moments. These can be obtained with the KSE for point processes (Eq. 86). For the mean in Eq. (99), we use $\phi(x)=x$:

[TABLE]

To compute the SDE for the posterior variance, we use Itô’s lemma:

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Adding these up gives us Eq. (100).