Joint Estimation and Classification

Problems involving dynamical continuous and static discrete variables

1. Introduction
In the previous chapters we have presented solutions for dynamic inference problems with either
discrete or continuous unknown variables. In this and coming chapters we will focus on problems
including both discrete and continuous variables, where at least one of the two involves dynamically
changing variables. In this chapter, we concentrate on problems involving static discrete variable(s) and
dynamic continuous.

As an example, consider a sensor producing measurements of a moving object. In the previous chapters
we have seen that the measurements contain (hidden) information about the movement of the object
and we have derived estimators like the Kalman Filter and the Particle Filter to estimate the state of the
object, thereby uncovering the hidden movement information.

Now let the object be of a certain type that is however unknown to the observer, and let the movement
of the object and the measurement of the sensor carry somehow information about the object type as
well. In case the object type is not known, this problem can then be considered as a joint estimation and
classification problem. The classification of the object can be done if the states of the object would be
known and the state of the object could be estimated if the object class would be known, however both
are unknown.

The above problem is just one of many variations of joint estimation and classification problems. In this
chapter we will first introduce a structural solution to this problem and afterwards provide some
example problems that can be solved using the structural solution.

2. Derivation of the basic structural solution

2.1. Problem setting
Assume that the dynamics and measurement equations describing the system can be represented with

𝑠𝑘+1 = 𝑓(𝑠𝑘 , 𝑐, 𝑤𝑘 )
𝑧𝑘 = ℎ(𝑠𝑘 , 𝑐, 𝑛𝑘 )

where the probability density functions of the process and measurement noise and the initial state are
𝑝(𝑤𝑘 ⁄𝑐), 𝑝(𝑛𝑘 ⁄𝑐), and 𝑝(𝑠0 ⁄𝑐) and let them all be dependent on the unknown class. These equations
look similar to the ones used in the formulation of the (regular) estimation problem, involving unknown
dynamic continuous variables. However the equations now include the (unknown) object type variable 𝑐,
which reflects that both object dynamics and measurements, as well as the prior density on the initial
state vector can now all depend in the unknown object class.

1
This material has been prepared by Hans Driessen and is part of the course Sensor Signal and Data Processing (EE5020) provided at EEMCS, TU-
Delft.
Assuming the prior class probabilities to be known and represented with 𝑃(𝑐), the objective is to
compute recursively the posterior joint pdf 𝑝(𝑠𝑘 , 𝑐⁄𝑍𝑘 ).

The fact that we have a joint problem, involving continuous and discrete variables does not imply that
the user will be interested in the estimation of both variables.

Some of the problems will be classification/detection problems where the user may not be interested
directly in the state of the object. Maybe the user is only interested in the class of the object but is not
interested in where the object is.

Vice versa, the user may be interested in the state vector, however not in the class.

2.2. Basic structural solution

After receiving measurements 𝑍𝑘 = {𝑧𝑘 , 𝑧𝑘−1 , … , 𝑧1 }, all information about the unknown variables 𝑠𝑘 , 𝑐
is contained in the joint density 𝑝(𝑠𝑘 , 𝑐⁄𝑍𝑘 ). Many algorithms for computing this joint density can be
derived by decomposing this density as a product of two terms

𝑝(𝑠𝑘 , 𝑐⁄𝑍𝑘 ) = 𝑝(𝑠𝑘 ⁄𝑐, 𝑍𝑘 )𝑃(𝑐⁄𝑍𝑘 )

The first term on the right hand side is the posterior probability density of the continuous state vector
conditioned on the unknown class 𝑐, and the second term is the posterior probability distribution of the
unknown discrete class variable 𝑐. This formulation stresses the jointness of the estimation and
classification problem.

It must be admitted that the algorithms can also be derived from the viewpoint of classification or
estimation alone, where we would start with deriving the marginal distribution 𝑝(𝑐⁄𝑍𝑘 ) on the type
variable or the marginal density of the state variable 𝑝(𝑠𝑘 ⁄𝑍𝑘 ), respectively. Especially the first approach
might seem to make more sense if we have to deal with detection problems, as we will see in one of the
examples. However even if we follow such an approach, we finally get the other term anyway, since we
cannot solve one problem without the other.

2.2.1. The conditional posterior density

The conditional posterior density can be written using the Bayesian evolution and update equations

𝑝(𝑠𝑘+1 ⁄𝑐, 𝑍𝑘 ) = ∫ 𝑝(𝑠𝑘+1⁄𝑠𝑘 , 𝑐)𝑝(𝑠𝑘 ⁄𝑐, 𝑍𝑘 ) 𝑑𝑠𝑘

𝒮
𝑝(𝑧𝑘 ⁄𝑠𝑘 , 𝑐)𝑝(𝑠𝑘 ⁄𝑐, 𝑍𝑘−1 )
𝑝(𝑠𝑘 ⁄𝑐, 𝑍𝑘 ) =
𝑝(𝑧𝑘 ⁄𝑐, 𝑍𝑘−1 )

where 𝑝(𝑧𝑘 ⁄𝑐, 𝑍𝑘−1 ) is the conditional likelihood of the current measurement given the type variable 𝑐
and the previous measurement, which acts as a normalizing constant, and can be computed via

𝑝(𝑧𝑘 ⁄𝑐, 𝑍𝑘−1 ) = ∫ 𝑝(𝑧𝑘 ⁄𝑠𝑘 , 𝑐, 𝑍𝑘−1 )𝑃(𝑠𝑘 ⁄𝑐, 𝑍𝑘−1 ) 𝑑𝑠𝑘
𝒮

2
This material has been prepared by Hans Driessen and is part of the course Sensor Signal and Data Processing (EE5020) provided at EEMCS, TU-
Delft.
This recursion is similar to the recursion we have derived for the problem involving dynamically changing
continuous variables. In that problem we had to deal with one specific model only and we did not
include that condition in the recursion. Now since the model depends on the type variable 𝑐, we have to
include it to make explicit the type. Actually this recursion describes 𝐶 recursions, one for each of the
types we have.

2.2.2. The posterior probability

The posterior probability can be computed via the recursion

𝑝(𝑧𝑘 ⁄𝑐, 𝑍𝑘−1 )𝑃(𝑐⁄𝑍𝑘−1 )

𝑃(𝑐⁄𝑍𝑘 ) =
𝑝(𝑧𝑘 ⁄𝑍𝑘−1 )

where 𝑝(𝑧𝑘 ⁄𝑍𝑘−1 ) is the likelihood of the current measurement given the previous measurement that
can be computed via
𝐶
𝑝(𝑧𝑘 ⁄𝑍𝑘−1 ) = ∑ 𝑝(𝑧𝑘 ⁄𝑐, 𝑍𝑘−1 )𝑃(𝑐⁄𝑍𝑘−1 )
𝑐=1

This recursion is almost equal to the Bayesian Classifier recursion we have derived for the problem with
unknown static discrete variables, except for the term representing the conditional likelihood of the
current measurement 𝑝(𝑧𝑘 ⁄𝑐, 𝑍𝑘−1 ). In the Bayesian Classifier we assumed independence of the
measurements and therefore the condition on the previous measurements could be removed. Here the
dependency is introduced due to the presence of the unknown dynamic continuous variables 𝑠𝑘 .

Notice the important role of the conditional likelihood in the solution to solving the joint problem.
Where in the problem involving only dynamic continuous unknown variable the term just seemed to play
the role of normalizing constant in the recursion for the posterior density, it now is the determining
factor in computing the posterior class probabilities.

2.2.3. Remarks
The basic structural solution is now clear. It comprises running 𝐶 state estimators, each tuned to one of
the specific object classes. These estimators deliver the likelihood on the unknown object class given the
measurements to a sequential Bayesian Classifier that computes the posterior probabilities. A pictorial
representation of the basic overall scheme leading to the solution of the joint problem is given in the
below figure.

3
This material has been prepared by Hans Driessen and is part of the course Sensor Signal and Data Processing (EE5020) provided at EEMCS, TU-
Delft.
 𝑝(𝑧𝑘 ⁄𝑐 = 1, 𝑍𝑘−1 )

Filter 1

Bayes Classifier
𝑃(𝑐 ⁄𝑍𝑘 )
𝑝(𝑠𝑘 ⁄𝑐 = 1, 𝑍𝑘 )
𝑧𝑘

𝑝(𝑧𝑘 ⁄𝑐 = 𝐶, 𝑍𝑘−1 )
Filter C

𝑝(𝑠𝑘 ⁄𝑐 = 𝐶, 𝑍𝑘 )

For the estimation algorithms (the filters in the figure), we might choose the Kalman Filter (or one of its
variants) when dealing with linear Gaussian problems, or we may choose a Particle Filter in case the
problems are highly nonlinear and/or non-Gaussian. For each of the estimation problems we can make
our own choice depending on which algorithm achieves the best trade-off between accuracy and
computational load.

2.3. Special cases

The above fully joint problem description is generic and comprises some more special cases. Let us
assume that the measurements can be split into two sets 𝑧𝑘 = {𝑧𝑘𝑠 , 𝑧𝑘𝑐 }, where 𝑧𝑘𝑠 represents the
measurement of the state vector 𝑠𝑘 via 𝑧𝑘𝑠

𝑠𝑘+1 = 𝑓(𝑠𝑘 , 𝑐, 𝑤𝑘 )
𝑧𝑘𝑠 = ℎ(𝑠𝑘 , 𝑐, 𝑛𝑘 )

and where the measurements 𝑧𝑘𝑐 represent measurements (or features) that are independent of the
state vector though directly dependent on the class, i.e.

𝑧𝑘𝑐 ~𝑝(𝑧𝑘𝑐 ⁄𝑐 )

assuming temporal independency of the features.

2.3.1. State dependent and independent measurements

Now the structural solution simplifies into

𝑝(𝑠𝑘 , 𝑐⁄𝑍𝑘 ) = 𝑝(𝑠𝑘 ⁄𝑐, 𝑍𝑘 )𝑃(𝑐⁄𝑍𝑘 )

= 𝑝(𝑠𝑘 ⁄𝑐, 𝑍𝑘𝑠 )𝑃(𝑐⁄𝑍𝑘 )
4
This material has been prepared by Hans Driessen and is part of the course Sensor Signal and Data Processing (EE5020) provided at EEMCS, TU-
Delft.
where the posterior 𝑝(𝑠𝑘 ⁄𝑐, 𝑍𝑘𝑠 ) depends only on the state vector related measurements and is
independent of the direct class dependent measurements.

The likelihood in the updating formula for 𝑃(𝑐 ⁄𝑍𝑘 ) can now be written as

𝑝(𝑧𝑘 ⁄𝑐, 𝑍𝑘−1 ) = 𝑝(𝑧𝑘𝑠 ⁄𝑐, 𝑍𝑘−1

𝑠 )𝑝(𝑧 𝑐 ⁄ )
𝑘 𝑐

Where the first part of the right hand side is the likelihood originating via the state dependent
measurements, while the second part are the likelihoods originating from the state independent
measurements.

2.3.2. Independent estimation and classification

Now consider additionally that the state evolution and measurement equations do not depent on the
class. Trivially this will further lead to a simplification in the sense that we now have two completely
independent problems, which is reflected in the equation corresponding with the structural solution
which comprises two independent terms

𝑝(𝑠𝑘 , 𝑐⁄𝑍𝑘 ) = 𝑝(𝑠𝑘 ⁄𝑍𝑘𝑠 )𝑝(𝑐⁄𝑍𝑘𝑐 )

In this case, the continuous variables and state dependent measurements do not contribute to the
solution to the classification problem, and the class and the state independent measurements do not
contribute to the solution of the estimation problem.

2.4. Computing the conditional likelihood 𝒑(𝒛𝒌 ⁄𝒁𝒌−𝟏 , 𝒄)

The class-conditional likelihood can be computed in different ways, depending on whether the
estimation problem is linear or not.

2.4.1. Linear Gaussian case

In the linear Gaussian case, the dynamic and measurement equations that depend on the unknown class
𝑐, are all linear and Gaussian

𝑠𝑘+1 = 𝐹 𝑐 𝑠𝑘 + 𝑤𝑘𝑐
𝑧𝑘 = 𝐻 𝑐 𝑠𝑘 + 𝑛𝑘𝑐

where the densities of process noise and measurement noise are zero mean Gaussian densities

𝑤𝑘𝑐 ~ 𝒩(𝑤𝑘𝑐 ; 0, 𝑊𝑘𝑐 )

𝑛𝑘𝑐 ~ 𝒩(𝑛𝑘𝑐 ; 0, 𝑁𝑘𝑐 )

and density of the initial state

𝑝(𝑠0 ⁄𝑐) = 𝒩(𝑠0 ; 𝑠0𝑐⁄0 , 𝑄0𝑐⁄0 )

All densities in the Bayesian recursion are therefore Gaussian densities and remain Gaussian, and can be
executed by running 𝐶 Kalman Filters that compute the mean and covariance of both predicted and
filtered probability density functions.

5
This material has been prepared by Hans Driessen and is part of the course Sensor Signal and Data Processing (EE5020) provided at EEMCS, TU-
Delft.
 𝑝(𝑠𝑘 ⁄𝑐, 𝑍𝑘 ) = 𝒩(𝑠𝑘 ; 𝑠𝑘𝑐⁄𝑘 , 𝑄𝑘𝑐⁄𝑘 )
𝑝(𝑠𝑘 ⁄𝑐, 𝑍𝑘−1 ) = 𝒩(𝑠𝑘 ; 𝑠𝑘𝑐⁄𝑘−1 , 𝑄𝑘𝑐⁄𝑘−1 )

Also the density of the measurements conditioned on the previous measurements and the class is
Gaussian

𝑝(𝑧𝑘 ⁄𝑍𝑘−1 , 𝑐) = ∫ 𝑝(𝑧𝑘 ⁄𝑠𝑘 , 𝑐 )𝑝(𝑠𝑘 ⁄𝑍𝑘−1 , 𝑐) 𝑑𝑠𝑘

𝒮
= 𝒩(𝑧𝑘 ; 𝜇𝑘𝑐 , 𝑉𝑘𝑐 )

with

𝜇𝑘𝑐 = 𝐻𝑠𝑘𝑐⁄𝑘−1
𝑉𝑘𝑐 = 𝐻𝑄𝑘𝑐⁄𝑘−1 𝐻 𝑇 + 𝑁𝑘

Therefore

1 1
𝑝(𝑧𝑘 ⁄𝑍𝑘−1 , 𝑐) = exp {− (𝑧𝑘 − 𝜇𝑘𝑐 )𝑇 (𝑉𝑘𝑐 )−1 (𝑧𝑘 − 𝜇𝑘𝑐 )}
2
√|2𝜋𝑉𝑘𝑐 |

This can also be derived by looking at the measurement equation which is the sum of two terms 𝐻𝑠𝑘 and
𝑛𝑘 . Conditioned on the previous measurements (and the class), the density of 𝑠𝑘 (𝑝(𝑠𝑘 ⁄𝑍𝑘−1 , 𝑐 )) is
Gaussian with mean and covariance 𝑠𝑘𝑐⁄𝑘−1 , 𝑄𝑘𝑐⁄𝑘−1 , leading to the above result (without the need for
solving the integral).

2.4.2. Nonlinear non-Gaussian case

Assuming we have available a particle cloud representation {𝑠𝑘𝑐,𝑖 , 𝛾𝑘−1
𝑐,𝑖
} of the density 𝑝(𝑠𝑘 ⁄𝑍𝑘−1 , 𝑐)

𝑝(𝑧𝑘 ⁄𝑍𝑘−1 , 𝑐) = ∫ 𝑝(𝑧𝑘 ⁄𝑠𝑘 , 𝑐 )𝑝(𝑠𝑘 ⁄𝑍𝑘−1 , 𝑐) 𝑑𝑠𝑘

𝒮
𝑁𝑝
≈ ∑ 𝑝(𝑧𝑘 ⁄𝑠𝑘𝑐,𝑖 , 𝑐)𝛾𝑘−1
𝑐,𝑖
𝑖=1

Such a particle cloud representation is available in a running particle filter, for instance in the SIR-filter
after predicting all particles towards the time validity of the received measurement.

Thus, in case of highly nonlinear non-Gaussian problems, a particle filter can be applied as part of the
structural solution, where the approximate likelihood can be calculated as detailed above.

3. Examples
3.1. Accurate object dependent tracking
Assume that we are tracking objects of different type where the only difference between these objects is
their randomness in the movement. This can be modeled according to
6
This material has been prepared by Hans Driessen and is part of the course Sensor Signal and Data Processing (EE5020) provided at EEMCS, TU-
Delft.
 𝑠𝑘+1 = 𝐹 𝑠𝑘 + 𝑤𝑘𝑐
𝑧𝑘 = 𝐻 𝑠𝑘 + 𝑛𝑘

where the dependency on the object type is only in the process noise 𝑤𝑘𝑐 .

We might not be interested in the object type itself, but we might be interested in accurate tracking.
Accurate tracking is possible only when we know the object type.

This problem falls into the category of problems of joint estimation and classification. Also the generic
solution that we developed applies here: we can apply a set of Kalman filters, each tuned to the
corresponding model of the process noise and apply the Bayesian algorithm to compute the posterior
class probability.

A question that we did not discuss yet is what is the output of this generic algorithm with respect to the
continuous state variable. We probably might want to report
𝐶

𝑝(𝑠𝑘 ⁄𝑍𝑘 ) = ∑ 𝑝(𝑠𝑘 ⁄𝑐, 𝑍𝑘 )𝑝(𝑐⁄𝑍𝑘 )

𝑐=1

For the linear problem this density is a weighted sum of Gaussians. If the user is interested in the mean
and covariance, we can compute those numbers using earlier derived expressions for the mean and
covariance of a weighted sum of Gaussians.

3.2. Moving object detection; plot-based linear Gaussian case

In this example we extend the previously discussed sequential object detection problem with
measurements of the state of the moving object. The derivation in this example does not directly follow
the structural solution presented above. It merely starts deriving the probabilities of detection via the
Sequential Probability Ratio Test. However, as we will see, the unknown object state is part of that
formulation. This can be “solved” by marginalizing the state from the detection part of the problem;
however in order to do so, leads to a (Kalman) filtering element. The final solution is exactly equal to the
solution obtained via the structural solution.

3.2.1. Problem specification

Consider a sensor that reports each scan whether it has detected an object or not (i.e. similarly to the
sequential object detection example with the Bayes Classifier). In addition, if it reports a detection, it
also reports a measurement of the unknown state object, be it its position and or velocity. Let the sensor
report be represented with 𝑧𝑘 ∈ {∅, 𝑦𝑘 }, where 𝑧𝑘 = 𝑦𝑘 represents a detection with 𝑦𝑘 the
measurement, and 𝑧𝑘 = ∅ represents no detection. In case there is an object, let the probability of
reporting a detection be represented with 𝑃𝐷 . Conversely, in case of no object, let the probability of
falsely reporting a detection be represented with 𝑃𝐹 . Assuming that 𝐻1 and 𝐻0 represent the hypothesis
that the object is present and absent, respectively, each valid with known prior probabilities 𝑃(𝐻1 ) and
𝑃(𝐻0 ) = 1 − 𝑃(𝐻1 ).

7
This material has been prepared by Hans Driessen and is part of the course Sensor Signal and Data Processing (EE5020) provided at EEMCS, TU-
Delft.
With the above modeling assumptions the likelihoods of the current measurement conditional to both
object present and absent hypotheses are

1 − 𝑃𝐷 , 𝑧𝑘 = ∅
𝑝(𝑧𝑘 ⁄𝑠𝑘 , 𝐻1 ) = {
𝑃𝐷 𝑝(𝑦𝑘 ⁄𝑠𝑘 , 𝐻1 ), 𝑧𝑘 = 𝑦𝑘

and

1 − 𝑃𝐹 , 𝑧𝑘 = ∅
𝑝(𝑧𝑘 ⁄𝐻0 ) = {
𝑃𝐹 𝑝(𝑦𝑘 ⁄𝐻0 ), 𝑧𝑘 = 𝑦𝑘

Here 𝑝(𝑦𝑘 ⁄𝑠𝑘 , 𝐻1 ) represents the pdf of the measurement given the unknown object state 𝑠𝑘 , and
𝑝(𝑦𝑘 ⁄𝐻0 ) represents the pdf of a false measurement.

It will be clear that forming the likelihood ratio (or taking the ratio of the posterior probabilities, is
complicated due to the presence of the unknown state in the likelihood for the object present
hypothesis.

For solving the above problem, assume further that the object movement and measurement can be
represented with the following stochastic, linear and Gaussian evolution and measurement equations

𝑠𝑘+1 = 𝐹𝑠𝑘 + 𝑤𝑘
𝑦𝑘 = 𝐻𝑠𝑘 + 𝑛𝑘

Where 𝑝(𝑤𝑘 ), 𝑝(𝑛𝑘 ), and 𝑝(𝑠0 ) are known Gaussian probability density functions.

3.2.2. Derivation of the Sequential Probability Ratio Detection algorithm

The above problem can be cast in the joint estimation and classification framework. However it does not
seem to make sense to derive the solution in the exact same way as before, since under the null-
hypothesis 𝐻0 there is no unknown state. Instead we apply a slightly different derivation that however
leads to the same structural solution; illustrating the point we made before.

We want to sequentially compute the posterior probability on the hypotheses conditioned on the
received measurements

𝑝(𝑧𝑘 ⁄𝐻𝑖 , 𝑍𝑘−1 )𝑃(𝐻𝑖 ⁄𝑍𝑘−1 )

𝑃(𝐻𝑖 ⁄𝑍𝑘 ) =
𝑝(𝑧𝑘 ⁄𝑍𝑘−1 )

with normalization factor the likelihood on the measurement conditioned on the previous
measurements 𝑝(𝑧𝑘 ⁄𝑍𝑘−1 ) = ∑1𝑖=0 𝑝(𝑧𝑘 ⁄𝐻𝑖 , 𝑍𝑘−1 )𝑃(𝐻𝑖 ⁄𝑍𝑘−1 ).

In that equation the various likelihoods can be computed via

𝑝(𝑧𝑘 ⁄𝐻0 , 𝑍𝑘−1 ) = 𝑝(𝑧𝑘 ⁄𝐻0 )

which was detailed above, and

8
This material has been prepared by Hans Driessen and is part of the course Sensor Signal and Data Processing (EE5020) provided at EEMCS, TU-
Delft.
 1 − 𝑃𝐷 , 𝑧𝑘 = ∅
𝑝(𝑧𝑘 ⁄𝐻1 , 𝑍𝑘−1 ) = {
𝑃𝐷 𝑝(𝑦𝑘 ⁄𝑌𝑘−1 , 𝐻1 ), 𝑧𝑘 = 𝑦𝑘

In this way we have marginalized the unknown state from the detection problem via

𝑝(𝑦𝑘 ⁄𝑌𝑘−1 , 𝐻1 ) = ∫ 𝑝(𝑦𝑘 ⁄𝑠𝑘 , 𝐻1 ) 𝑝(𝑠𝑘 ⁄𝑌𝑘−1 , 𝐻1 )𝑑𝑠𝑘

𝑆

This conditional likelihood 𝑝(𝑦𝑘 ⁄𝑌𝑘−1 , 𝐻1 ) can be computed from the results of a running Kalman Filter
applied to the evolution and measurement models described above.

𝑝(𝑦𝑘 ⁄𝑌𝑘−1 , 𝐻1 ) = 𝒩(𝑦𝑘 ; 𝜇𝑘 , 𝑉𝑘 )

with

𝜇𝑘 = 𝐻𝑠𝑘⁄𝑘−1
𝑉𝑘 = 𝐻𝑄𝑘⁄𝑘−1 𝐻𝑇 + 𝑁𝑘

that can both be computed from a running Kalman Filter.

With this the probability ratio between the posterior probabilities can be recursively computed via

𝑃(𝐻1 ⁄𝑍𝑘 )
𝑅(𝑍𝑘 ) ≡
𝑃(𝐻0 ⁄𝑍𝑘 )
𝑃(𝑧𝑘 ⁄𝐻1 ) 𝑃(𝐻1 ⁄𝑍𝑘−1 )
=
𝑃(𝑧𝑘 ⁄𝐻0 ) 𝑃(𝐻0 ⁄𝑍𝑘−1 )
= 𝑟(𝑧𝑘 )𝑅(𝑍𝑘−1 )

with

1 − 𝑃𝐷
, 𝑧𝑘 = ∅
1 − 𝑃𝐹
𝑟(𝑧𝑘 ) =
𝒩(𝑦𝑘 ; 𝜇𝑘 , 𝑉𝑘 ) 𝑃𝐷
, 𝑧𝑘 = 𝑦𝑘
{ 𝑝(𝑦𝑘 ⁄𝐻0 ) 𝑃𝐹

This ratio is often used in the Sequential Probability Ratio Test (often implemented in the logarithmic
form).

Notice the similarity with the previously derived expression for the sensor only reporting the detection,
𝒩(𝑦𝑘 ;𝜇𝑘 ,𝑉𝑘 )
however not the measurement 𝑦𝑘 , the only difference is in the additional term 𝑝(𝑦𝑘 ⁄𝐻0 )
.

Often the pdf of the false alarm measurements is taken as uniform, 𝑝(𝑦𝑘 ⁄𝐻0 ) = 𝛾. As long as the
measurements are such that

1
𝒩(𝑦𝑘 ; 𝜇𝑘 , 𝑉𝑘 ) >
𝛾

9
This material has been prepared by Hans Driessen and is part of the course Sensor Signal and Data Processing (EE5020) provided at EEMCS, TU-
Delft.
the term will support the object presence hypothesis and vice versa.

3.2.3. Equality with the structural solution

Although we did not derive the above algorithm using the structural derivation, it appears that the final
solution is mathematically equal to the structural solution.

A (Kalman) filter is running for the object present hypothesis which calculates the posterior and the
conditional likelihood 𝑝(𝑧𝑘 ⁄𝐻1 , 𝑍𝑘−1 ). Under the object absent hypothesis the computation of
𝑝(𝑧𝑘 ⁄𝐻0 , 𝑍𝑘−1 ) = 𝑝(𝑧𝑘 ⁄𝐻0 ) does not require a filter. A Bayesian Classifier is computing the posterior
probabilities 𝑃(𝐻𝑖 ⁄𝑍𝑘 ). The computation of the Probability Ratio is possible and meaningful since this is
a binary decision problem.

3.3. Moving object detection: Track-before-detect

Let us consider problem of object tracking in the context of track-before-detect, i.e. tracking based on
raw video measurements. This is the problem that we discussed applying the particle filter to nonlinear
problems. Here we extend that estimation problem to an estimation and detection problem.

3.3.1. TBD modeling

Let us repeat the modeling here. We have a similar dynamic model as earlier

𝑠𝑘+1 = 𝐹𝑠𝑘 + 𝐺𝑤𝑘

But now introduce a nonlinear measurement equation as

𝑧𝑘 = 𝑎𝑘 ℎ(𝑠𝑘 ) + 𝑛𝑘

The specific nonlinearity could arise in image processing, radar and sonar problems. The vector 𝑧𝑘 can
for instance represent an array of measured intensities (sound, light, radar reflections), ℎ(𝑠𝑘 ) represents
for instance a point spread function (in image processing) or an array response vector in audio, radar or
sonar, where 𝑎𝑘 is an unknown signal strength. The noise 𝑛𝑘 vector is assumed to be distributed as zero
mean Gaussian with covariance 𝜎𝑛2 𝐼.

Assuming that the prior for 𝑎𝑘 is given as a zero Gaussian density 𝑁(𝑎𝑘 ; 0, 𝜎𝑎2 ), the pdf of the
measurements conditioned on the state vector can be written as

1
𝑝(𝑧𝑘 ⁄𝑠𝑘 ) = exp{−𝑧𝑘𝑇 𝑄𝑧−1 𝑧𝑘 }
√|2𝜋𝑄𝑧 |

With
𝑇
𝑄𝑧 = ℎ(𝑠𝑘 )(ℎ(𝑠𝑘 )) 𝜎𝑎2 + 𝜎𝑛2 𝐼

3.3.2. Extension with detection/classification

We can extend this problem to a detection problem when we assume that we do not know a priori
whether the object is present in the measurements 𝑧𝑘 . We can model this using the well-known object

10
This material has been prepared by Hans Driessen and is part of the course Sensor Signal and Data Processing (EE5020) provided at EEMCS, TU-
Delft.
absent and present hypotheses. Under the model present hypothesis we run the Particle Filter, which
approximately computes the conditional likelihood
𝑁𝑝
𝑝(𝑧𝑘 ⁄𝑍𝑘−1 , 𝐻1 ) ≅ ∑ 𝑝(𝑧𝑘 ⁄𝑠𝑘𝑐,𝑖 , 𝑐)𝛾𝑘−1
𝑐,𝑖
𝑖=1

while the conditional likelihood under the object absent hypothesis is simply

𝑝(𝑧𝑘 ⁄𝑍𝑘−1 , 𝐻0 ) = 𝑝(𝑧𝑘 )

Which is the likelihood that the measurement vector just comprises noise.

It will be clear, also given the previous example that this can be implemented as a Sequential Probability
Ratio Test

𝑅(𝑍𝑘 ) = 𝑟(𝑧𝑘 )𝑅(𝑍𝑘−1 )

With
𝑁
𝑝
∑𝑖=1 𝑝(𝑧𝑘 ⁄𝑠𝑘𝑐,𝑖 , 𝑐)𝛾𝑘−1
𝑐,𝑖
𝑟(𝑧𝑘 ) =
𝑝(𝑧𝑘 )

3.3.3. Implementation aspects

[Link]. Likelihood ratio calculation

Whether we directly calculate the posterior probabilities 𝑃(𝐻𝑖 ⁄𝑍𝑘−1 ) or their ratio (in case of two
hypotheses), it is good to observe that all common terms in the likelihoods 𝑝(𝑧𝑘 ⁄𝑍𝑘−1 , 𝐻𝑖 ) vanish, not
just constants but also common data dependent terms. Let us have a look at the total likelihood
𝑁𝑝
𝑝(𝑧𝑘 ⁄𝑍𝑘−1 , 𝐻𝑖 ) = ∑ 𝑝(𝑧𝑘 ⁄𝑠𝑘𝑐,𝑖 , 𝑐)𝛾𝑘−1
𝑐,𝑖
𝑖=1

Notice that I am using a mixed notation, which is not very clear now with 𝑐 and 𝐻𝑖 , but I will resolve this
later. Especially concentrate on the likelihood under the signal presence hypothesis

1 𝑧𝑘𝑇 𝑄𝑧−1 𝑧𝑘
𝑝(𝑧𝑘 ⁄𝑠𝑘 , 𝐻1 ) = exp {− }
√|2𝜋𝑄𝑧 | 2

with the covariance matrix

𝑇
𝑄𝑧 = ℎ(𝑠𝑘 )(ℎ(𝑠𝑘 )) 𝜎𝑎2 + 𝜎𝑛2 𝐼

This covariance matrix actually could be class/hypothesis dependent, but let us ignore this for
convenience here. The inverse and determinant of this matrix can be re-written into simpler forms. Start
with the inverse matrix that can be re-written using the Matrix Inversion lemma
𝑇
𝜎𝑛−2 ℎ(𝑠𝑘 )(ℎ(𝑠𝑘 )) 𝜎𝑛−2
𝑄𝑧−1 = 𝜎𝑛−2 𝐼 − 𝑇
(ℎ(𝑠𝑘 )) 𝜎𝑛−2 ℎ(𝑠𝑘 ) + 𝜎𝑎−2
11
This material has been prepared by Hans Driessen and is part of the course Sensor Signal and Data Processing (EE5020) provided at EEMCS, TU-
Delft.
Using this equality

𝑇 2
𝑧𝑘𝑇 𝑧𝑘 1 (ℎ(𝑠𝑘 )) 𝑧𝑘
𝑧𝑘𝑇 𝑄𝑧−1 𝑧𝑘 = 2 − | |
𝜎𝑛 𝜒 𝜎𝑛2

𝑇
𝜒 = (ℎ(𝑠𝑘 )) 𝜎𝑛−2 ℎ(𝑠𝑘 ) + 𝜎𝑎−2

Then the determinant of the covariance matrix also can be re-written as in

𝑇
|𝑄𝑧 | = |ℎ(𝑠𝑘 )(ℎ(𝑠𝑘 )) 𝜎𝑎2 + 𝜎𝑛2 𝐼|
𝑇 𝜎𝑎2
= |𝜎𝑛2 𝐼| |(ℎ(𝑠𝑘 )) ℎ(𝑠𝑘 ) + 1|
𝜎𝑛2
= (𝜎𝑛2 )𝑁 𝜉

With

𝑇 𝜎𝑎2
𝜉 = ((ℎ(𝑠𝑘 )) ℎ(𝑠𝑘 ) + 1)
𝜎𝑛2

This implies that the likelihood can be written as the product

𝑇 2
1 𝑧𝑘𝑇 𝑧𝑘 1 1 (ℎ(𝑠𝑘 )) 𝑧𝑘
𝑝(𝑧𝑘 ⁄𝑠𝑘 , 𝐻1 ) = exp {− } exp { | | }
√(2𝜋𝜎𝑛2 )𝑁 2𝜎𝑛2 𝜉 2𝜒 𝜎𝑛2

where close inspection of the first two terms on the right-hand side can be recognized as the likelihood
under the assumption that the signal contains noise only, i.e.

1 𝑧𝑘𝑇 𝑧𝑘
𝑝(𝑧𝑘 ⁄𝐻0 ) = exp {− }
√(2𝜋𝜎𝑛2 )𝑁 2𝜎𝑛2

In the computations of now for instance the Sequential Probability Ratio, the likelihood ratio term
simplifies into

𝑇 2
𝑝(𝑧𝑘 ⁄𝑠𝑘 , 𝐻1 ) 1 1 (ℎ(𝑠𝑘 )) 𝑧𝑘
= exp { | | }
⁄
𝑝(𝑧𝑘 𝐻0 ) 𝜉 2𝜒 𝜎𝑛2

Next to the numerical robustness tricks that are discussed in the next subsection, this “trick” has been
used to implement the TBD particle filter in the Matlab example thereby preventing the inverse of a
potentially large covariance matrix, preventing very small likelihood values to appear and enabling the
computation of likelihoods/weights without a loop over all the particles.
12
This material has been prepared by Hans Driessen and is part of the course Sensor Signal and Data Processing (EE5020) provided at EEMCS, TU-
Delft.
 [Link]. Numerically robust calculation of weights and (log) likelihoods
Joint estimation and detection in TBD requires the computation of the particle filter weights and the
total likelihood under H1. In terms of the logarithm of the individual likelihoods

𝜆𝑖 = log{𝑝(𝑧𝑘 ⁄𝑠𝑘𝑖 , 𝐻1 )}

the weights can be written as

𝑖
exp{𝜆𝑖 }
𝑤 =
∑𝑗 exp{𝜆 𝑗 }

and the total likelihood under H1 can be written as

1
𝑝(𝑧𝑘 ⁄𝑍𝑘−1 , 𝐻1 ) = ∑ exp{𝜆𝑗 }
𝑁
𝑗

The problem in implementation is that often exp{𝜆𝑖 } leads to underflow for large negative values of 𝜆𝑖 ,
which prevents the particle filter from running and the classifier to work properly.

Solution:

Compute 𝜆𝑚𝑎𝑥 = max 𝜆𝑖 where then computing exp{𝜆𝑖 − 𝜆𝑚𝑎𝑥 } leads to less underflow, then
𝑖

exp{𝜆 𝑖 − 𝜆𝑚𝑎𝑥 } exp{𝜆𝑖 } exp{−𝜆𝑚𝑎𝑥 }

=
∑𝑗 exp{𝜆 𝑗 −𝜆𝑚𝑎𝑥 } ∑𝑗 exp{𝜆 𝑗 } exp{−𝜆𝑚𝑎𝑥 }
exp{𝜆𝑖 }
=
∑𝑗 exp{𝜆 𝑗 }
= 𝑤𝑖

and

1
log{𝑝(𝑧𝑘 ⁄𝑍𝑘−1 , 𝐻1 )} = log { ∑ exp{𝜆𝑗 }}
𝑁
𝑗

1
= log { ∑ exp{𝜆𝑗 − 𝜆𝑚𝑎𝑥 } exp{𝜆𝑚𝑎𝑥 }}
𝑁
𝑗

1
= log { ∑ exp{𝜆𝑗 − 𝜆𝑚𝑎𝑥 }} + log{exp{𝜆𝑚𝑎𝑥 }}
𝑁
𝑗

1
= log { ∑ exp{𝜆𝑗 − 𝜆𝑚𝑎𝑥 }} + 𝜆𝑚𝑎𝑥
𝑁
𝑗

13
This material has been prepared by Hans Driessen and is part of the course Sensor Signal and Data Processing (EE5020) provided at EEMCS, TU-
Delft.