Path sampling for particle filters with application to

arXiv:1009.2108v1 [[Link]] 10 Sep 2010

multi-target tracking
Vasileios Maroulas
Department of Mathematics
University of Tennessee
Knoxville, TN 37996
and
Panos Stinis
Department of Mathematics
University of Minnesota
Minneapolis, MN 55455

Abstract
In recent work [15], we have presented a novel approach for improv-
ing particle filters for multi-target tracking. The suggested approach
was based on Girsanov’s change of measure theorem for stochastic dif-
ferential equations. Girsanov’s theorem was used to design a Markov
Chain Monte Carlo step which is appended to the particle filter and
aims to bring the particle filter samples closer to the observations. In
the current work, we present an alternative way to append a Markov
Chain Monte Carlo step to a particle filter to bring the particle fil-
ter samples closer to the observations. Both current and previous
approaches stem from the general formulation of the filtering prob-
lem. We have used the currently proposed approach on the problem of
multi-target tracking for both linear and nonlinear observation models.
The numerical results show that the suggested approach can improve
significantly the performance of a particle filter.

Introduction
Multi-target tracking is a central and difficult problem arising in many sci-
entific and engineering applications including radar and signal processing,
air traffic control and GPS navigation [12]. The tracking problem consists
of computing the best estimate of the targets’ trajectories based on noisy

1
measurements (observations). Several strategies have been developed for
addressing the multi-target tracking problem, see e.g. [1, 6, 5, 16, 8, 11, 12,
13, 14, 22].
As in our recent work [15], in this paper we also focus on particle filter
techniques [5, 16]. The popularity of the particle filter method has increased
due to its flexibility to handle cases where the dynamic and observation
models are non-linear and/or non-Gaussian. The particle filter approach is
an importance sampling method which approximates the target distribution
by a discrete set of weighted samples (particles). The weights of the samples
are updated when observations become available in order to incorporate
information from the observations.
Despite the particle filter’s flexibility, it is often found in practice that
most samples will have a negligible weight with respect to the observation,
in other words their corresponding contribution to the target distribution
will be negligible. Therefore, one may resample the weights to create more
copies of the samples with significant weights [8]. However, even with the
resampling step, the particle filter might still need a lot of samples in order
to approximate accurately the target distribution. Typically, a few samples
dominate the weight distribution, while the rest of the samples are in sta-
tistically insignificant regions. Thus, some authors (see e.g. [7, 23]) have
suggested the use of an extra step, after the resampling step, which can help
move more samples in statistically significant regions.
The extra step for the particle filter is a problem of conditional path sam-
pling for stochastic differential equations (SDEs). In [19], a new approach
to conditional path sampling based on Girsanov’s theorem was presented.
In that paper, it was also shown how the algorithm can be used to per-
form the extra step of a particle filter. In [15], we applied the conditional
path sampling algorithm from [19] to perform the extra step of a particle
filter for the problem of multi-target tracking. The numerical results in [15]
suggested that the approach can improve significantly the performance of a
particle filter for multi-target tracking. In the current work, we show yet
another way of how to perform the extra step of a particle filter. Both the
current approach and the one in [15] stem from the general formulation of
the filtering problem. The details of the currently proposed implementation
of the extra step for a single target are given in Section 1.3 and for multiple
targets in Section 1.4. The relative merits of the proposed approach in this
paper and the one proposed in [15] are briefly discussed in Section 3. A
more detailed comparison will be presented elsewhere.
To address the target-observation association problem we have used a
simple Metropolis Monte Carlo algorithm which first appeared in [15]. This

2
algorithm effects a probabilistic search of the space of possible associations
to find the best target-observation association. Of course, one can use more
sophisticated association algorithms (see [16] and references therein) but the
Monte Carlo algorithm performed very well in the numerical experiments.
The paper is organized as follows. Sections 1.1 and 1.2 provide a brief
presentation of particle filters for single and multiple targets (more details
can be found in [8, 5, 10, 16]), which will serve to highlight the versatility
and drawbacks of this popular filtering method. Sections 1.3 and 1.4 demon-
strate how one can use an extra step to improve the performance of particle
filters for single and multiple targets. Section 2 presents numerical results
for multi-target tracking for the cases of linear and nonlinear observation
models. Finally, Section 3 contains a discussion of the results as well as
directions for future work.

1 Particle filtering
Particle filters are a special case of sequential importance sampling methods.
In Sections 1.1 and 1.2 we discuss the generic particle filter for a single and
multiple targets respectively. In Sections 1.3 and 1.4 we discuss the addition
of an extra step to the generic particle filter for the cases of a single and
multiple targets respectively.

1.1 Generic particle filter for a single target

Suppose that we are given an SDE system and that we also have access
to noisy observations ZT1 , . . . , ZTK of the state of the system at specified
instants T1 , . . . , TK . The observations are functions of the state of the sys-
tem, say given by ZTk = G(XTk , ξk ), where ξk , k = 1, . . . , K are mutually
independent random variables. For simplicity, let us assume that the distri-
bution of the observations admits a density g(XTk , ZTk ), i.e., p(ZTk |XTk ) ∝
g(XTk , ZTk ).
The filtering problem consists of computing estimates of the conditional
expectation E[f (XTk )|{ZTj }kj=1 ], i.e., the conditional expectation of the
state of the system given the (noisy) observations. Equivalently, we are
looking to compute the conditional density of the state of the system given
the observations p(XTk |{ZTj }kj=1 ). There are several ways to compute this
conditional density and the associated conditional expectation but for prac-
tical applications they are rather expensive.
Particle filters fall in the category of importance sampling methods. Be-
cause computing averages with respect to the conditional density involves

3
the sampling of the conditional density which can be difficult, importance
sampling methods proceed by sampling a reference density q(XTk |{ZTj }kj=1 )
which can be easily sampled and then compute the weighted sample mean
N
1 X p(XTnk |{ZTj }kj=1 )
E[f (XTk )|{ZTj }kj=1 ] ≈ f (XTnk )
N
n=1
q(XTn |{ZTj }kj=1 )
k

or the related estimate

PN p(XTn |{ZTj }kj=1 )
n
n=1 f (XTk ) q(X n |{Z
k
k
Tj }j=1 )
E[f (XTk )|{ZTj }kj=1 ] ≈
Tk
n k
, (1)
PN p(XT |{ZTj }j=1 )
k
n=1 q(XTn |{ZTj }kj=1 )
k

where N has been replaced by the approximation

N
X p(XTnk |{ZTj }kj=1 )
N≈ .
n=1
q(XTnk |{ZTj }kj=1 )

Particle filtering is a recursive implementation of the importance sampling

approach. It is based on the recursion

p(XTk |{ZTj }kj=1 ) ∝ g(XTk , ZTk )p(XTk |{ZTj }j=1

k−1
), (2)
Z
k−1 k−1
where p(XTk |{ZTj }j=1 ) = p(XTk |XTk−1 )p(XTk−1 |{ZTj }j=1 )dXTk−1 . (3)

If we set
q(XTk |{ZTj }kj=1 ) = p(XTk |{ZTj }j=1
k−1
),
then from (2) we get

p(XTk |{ZTj }kj=1 )

∝ g(XTk , ZTk ).
q(XTk |{ZTj }kj=1 )

The approximation in expression (1) becomes

PN n n
n=1 f (XTk )g(XTk , ZTk )
E[f (XTi )|{ZTj }kj=1 ] ≈ PN (4)
n
n=1 g(XTk , ZTk )

From (4) we see that if we can construct samples from the predictive distri-
k−1
bution p(XTk |{ZTj }j=1 ) then we can define the (normalized) weights WTnk =

4
 g(XTn ,ZTk )
PN k
n , use them to weigh the samples and the weighted samples
n=1 g(XT ,ZTk )
k
will be distributed according to the posterior distribution p(XTk |{ZTj }kj=1 ).
In many applications, most samples will have a negligible weight with
respect to the observation, so carrying them along does not contribute sig-
nificantly to the conditional expectation estimate (this is the problem of
degeneracy [10]). To create larger diversity one can resample the weights to
create more copies of the samples with significant weights. The particle filter
with resampling is summarized in the following algorithm due to Gordon et
al. [8].

Particle filter for a single target

k−1
1. Begin with N unweighted samples XTnk−1 from p(XTk−1 |{ZTj }j=1 ).

2. Prediction: Generate N samples XT′nk from p(XTk |XTk−1 ).

3. Update: Evaluate the weights

g(XT′nk , ZTk )
WTnk = PN .
′n
n=1 g(XTk , ZTk )

4. Resampling: Generate N independent uniform random variables

′j
{θ n }N n
n=1 in (0, 1). For n = 1, . . . , N let XTk = XTk where

j−1
X j
X
WTl k ≤θ <j
WTl k
l=1 l=1

where j can range from 1 to N.

5. Set k = k + 1 and proceed to Step 1.

The particle filter algorithm is easy to implement and adapt for dif-
ferent problems since the only part of the algorithm that depends on the
specific dynamics of the problem is the prediction step. This has led to the
particle filter algorithm’s increased popularity [5]. However, even with the
resampling step, the particle filter can still need a lot of samples in order to
describe accurately the conditional density p(XTk |{ZTj }kj=1 ). Snyder et al.
[18] have shown how the particle filter can fail in simple high dimensional
problems because one sample dominates the weight distribution. The rest
of the samples are not in statistically significant regions. Even worse, as we

5
will show in the numerical results section, there are simple examples where
not even one sample is in a statistically significant region. In the next sub-
section we present how an extra step can be used to push samples closer to
statistically significant regions.

1.2 Generic particle filter for multiple targets

Suppose that we have λ = 1, . . . , Λ targets. Also, for notational simplicity,
assume that the λth target comes from the λth observation. Even when this
is not the case, we can relabel the observations to satisfy this assumption.
The targets are assumed to evolve independently so that the observation
weight of a sample of the vector of targets is the product of the individual
observation weights of the targets [16]. The same is true for the transition
density of the vector of targets between observations. We denote the vector
of targets at observation Tk by

XTk = (X1,Tk , . . . , XΛ,Tk )

and the observation vector at Tk by

ZTk = (Z1,Tk , . . . , ZΛ,Tk ).

Also, we can have different observation weight densities gλ , λ = 1, . . . , Λ for

different targets. However, in the numerical examples we have chosen the
same observation weight density for all targets.
Following [16] we can write the particle filter for the case of multiple
targets as

Particle filter for multiple targets

k−1
1. Begin with N unweighted samples XTnk−1 from p(XTk−1 |{ZTj }j=1 )=
QΛ k−1
λ=1 p(Xλ,Tk−1 |{Zλ,Tj }j=1 ).

2. Prediction: Generate N samples XT′nk from

Λ
Y
p(XTk |XTk−1 ) = p(Xλ,Tk |Xλ,Tk−1 ).
λ=1

3. Update: Evaluate the weights

′n
QΛ
λ=1 gλ (X λ,Tk , Zλ,Tk )
WTnk = PN QΛ .
g (X ′n , Z )
n=1 λ=1 λ λ,Tk λ,Tk

6
 4. Resampling: Generate N independent uniform random variables
′j
{θ n }N n
n=1 in (0, 1). For n = 1, . . . , N let XTk = XTk where

j−1
X j
X
WTl k ≤ θ j < WTl k
l=1 l=1

where j can range from 1 to N.

5. Set k = k + 1 and proceed to Step 1.

1.3 Particle filter with MCMC step for a single target

Several authors (see e.g. [7, 23]) have suggested the use of a MCMC step
after the resampling step (Step 4) in order to move samples away from
statistically insignificant regions. There are many possible ways to append
an MCMC step after the resampling step in order to achieve that objective.
The important point is that the MCMC step must preserve the conditional
density p(XTk |{ZTj }kj=1 ). In the current section we show that the MCMC
step constitutes a case of conditional path sampling.
We begin by noting that one can use the resampling step (Step 4) in the
particle filter algorithm to create more copies not only of the good samples
according to the observation, but also of the values (initial conditions) of the
samples at the previous observation. These values are the ones who have
evolved into good samples for the current observation (see more details in
[23]). The motivation behind producing more copies of the pairs of initial
and final conditions is to use the good initial conditions as starting points
to produce statistically more significant samples according to the current
observation. This process can be accomplished in two steps. First, Step 4
of the particle filter algorithm is replaced by

Resampling: Generate N independent uniform random variables {θ n }N n=1

in (0, 1). For n = 1, . . . , N let (XTnk−1 , XTnk ) = (XT′jk−1 , XT′jk )where

j−1
X j
X
WTl k ≤θ <j
WTl k
l=1 l=1

To motivate the second step we use the recursive particle filter formulas (2),

7
(3) to write an update equation for conditional expectations. We have

E[f (XTk )|{ZTj }kj=1 ] =

k−1
g(XTk , ZTk )p(XTk−1 |{ZTj }j=1 )dXTk−1 dP
Z
f (XTk ) R k−1
, (5)
g(XTk , ZTk )p(XTk−1 |{ZTj }j=1 )dXTk−1 dP

where P is the Wiener measure. Formula (5) can be approximated (using

the initial conditions of the samples produced by the modified resampling
step above) as follows:
N Z
1 X g(XTnk , ZTk )dP
E[f (XTk )|{ZTj }kj=1 ] ≈ n R
f (XTk ) , (6)
N g(XTk , ZTk )dP
n=1

i.e., we have approximated the integration over the initial conditions XTk−1
appearing in (5) by an average over the initial conditions produced by the
modified resampling step. Note that we have not approximated the integra-
tion over XTk−1 appearing in the denominator because the denominator is
a normalization constant which will not be needed in the MCMC sampling.
Formula (6) allows us to perform the MCMC step which will bring the
samples closer to the current observation. In particular, starting from the
initial conditions at the previous observation which were picked by the
resampling step, we can sample the density g(XTnk , ZTk ) with respect to
the Wiener measure dP and produce samples XTnk which are more signifi-
cant with respect to the current observation. The sampling of the density
g(XTnk , ZTk ) with respect to the Wiener measure dP corresponds to the sam-
pling of the noise process (Wiener measure) so that the point XTnk has high
probability with respect to g(XTnk , ZTk ). In other words, sampling the density
g(XTnk , ZTk ) given the initial conditions at the previous observation is noth-
ing else but a case of conditional sampling of a path that starts from some
given initial condition XTnk−1 and ends at a point XTnk which has high prob-
ability with respect to g(XTnk , ZTk ). Of course, the sampled path must also
obey the dynamic model (see also [23] for a related approach for improving
particle filters).
It is obvious that this sampling procedure preserves the conditional den-
sity p(XTk |{ZTj }kj=1 ) since there was no approximation involved (other than
the necessary approximation of the integration over the initial conditions by
a discrete summation). The details of how to perform the density sampling
will be explained in Section 2.1.
Finally, note that the approximation in equation (6) involves an inte-
gration over the Wiener measure. This is because for each of the initial

8
conditions at the previous observation one can create different paths de-
pending on the choice of the Brownian path. The integration over dP can
be approximated by an average over appropriately weighted Brownian paths
and we find
N M
1 XX
E[f (XTk )|{ZTj }kj=1 ] ≈ f (YTn,m ), (7)
N M n=1 m=1 k

The approximation in (7) has an error of order O( √N1 M ). This means that
even if M is kept finite, the approximation still converges in the infinite
limit for N. In the numerical implementation we pick M = 1 and focus our
computational resources on the averaging over initial conditions.
We are now in a position to present the particle filter with MCMC step
algorithm

Particle filter with MCMC step for a single target

k−1
1. Begin with N unweighted samples XTnk−1 from p(XTk−1 |{ZTj }j=1 ).

2. Prediction: Generate N samples XT′nk from p(XTk |XTk−1 ).

3. Update: Evaluate the weights

g(XT′nk , ZTk )
WTnk = PN .
′n
n=1 g(XTk , ZTk )

4. Resampling: Generate N independent uniform random variables

′j ′j
{θ n }N n n
n=1 in (0, 1). For n = 1, . . . , N let (XTk−1 , XTk ) = (XTk−1 , XTk )
where
j−1
X X j
WTl k ≤ θ j < WTl k
l=1 l=1

where j can range from 1 to N.

5. MCMC step: Construct a Markov chain {YTn,l

k
}Λ
l=0 with initial value
n,0 n
YTk = XTk and stationary distribution

g(Y, ZTk )dP

R .
g(XTk , ZTk )dP

6. Set XTnk = YTn,Λ

k
.

9
 7. Set k = k + 1 and proceed to Step 1.

Since the samples XTnk = YTn,Λ

k
are constructed by starting from different
sample paths, they are independent. Also, note that the samples XTnk are
unweighted. However, we can still measure how well these samples approx-
imate the posterior density by comparing the effective sample sizes of the
particle filter with and without the MCMC step. For a collection of N
samples the effective sample size ess(Tk ) is defined by

N
ess(Tk ) =
1 + Ck2

where
v
u N N
1 u
t1
X 1 X
Ck = (g(XTnk , ZTk ) − Wk )2 and Wk = g(XTnk , ZTk ).
Wk N n=1 N n=1

The effective sample size can be interpreted as that the N weighted samples
N
are worth of ess(Tk ) = 1+C 2 i.i.d. samples drawn from the target density,
k
which in our case is the posterior density. By definition, ess(Tk ) ≤ N. If the
samples have uniform weights, then ess(Tk ) = N. On the other hand, if all
samples but one have zero weights, then ess(Tk ) = 1.

1.4 Particle filter with MCMC step for multiple targets

We discuss now the case of multiple, say Λ, targets. Instead of the obser-
vations for a single target now we have a collection of observations for all
the targets {ZTj }kj=1 = {(ZT1j , . . . , ZTλj )}kj=1 . There are two cases which we
now examine. First, the function whose conditional expectation we want to
compute depends only on one of the targets. Second, the function whose
conditional expectation we want to compute depends at least on two targets.
We begin with the first case. Recall that the targets evolve independently
and after the appropriate association between a target and an observation,
we have for the conditional expectation estimate for the λth target (λ =

10
1, . . . , Λ)

E[f (Xλ,Tk )|{ZTj }kj=1 ]

QΛ k−1
λ′ =1 gλ′ (Xλ′ ,Tk , Zλ′ ,Tk )p(Xλ′ ,Tk−1 |{Zλ′ ,Tj }j=1 )dXλ′ ,Tk−1 dPλ′
Z
= f (Xλ,Tk ) R QΛ k−1
λ′ =1 gλ′ (Xλ′ ,Tk , Zλ′ ,Tk )p(Xλ′ ,Tk−1 |{Zλ′ ,Tj }j=1 )dXλ′ ,Tk−1 dPλ′
k−1
gλ (Xλ,Tk , Zλ,Tk )p(Xλ,Tk−1 |{Zλ,Tj }j=1 )dXλ,Tk−1 dPλ
Z
= f (Xλ,Tk ) R k−1
gλ (Xλ,Tk , Zλ,Tk )p(Xλ,Tk−1 |{Zλ,Tj }j=1 )dXλ,Tk−1 dPλ
(8)

where the second identity comes from integrating out all the targets except
for the λth target. Due to this simplification we are back at the case of
a single target as in equation (5). So, we can perform the MCMC step
individually for each target. The particle filter with MCMC step for the
case of multiple targets is

Particle filter with MCMC step for multiple targets

k−1
1. Begin with N unweighted samples XTnk−1 from p(XTk−1 |{ZTj }j=1 )=
QΛ k−1
λ=1 p(Xλ,Tk−1 |{Zλ,Tj }j=1 ).

2. Prediction: Generate N samples XT′nk from

Λ
Y
p(XTk |XTk−1 ) = p(Xλ,Tk |Xλ,Tk−1 ).
λ=1

3. Update: Evaluate the weights

′n
QΛ
λ=1 gλ (X λ,Tk , Zλ,Tk )
WTnk = PN Q Λ .
′n
n=1 λ=1 gλ (X λ,Tk , Zλ,Tk )

4. Resampling: Generate N independent uniform random variables

′j ′j
{θ n }N n n
n=1 in (0, 1). For n = 1, . . . , N let (XTk−1 , XTk ) = (XTk−1 , XTk )
where
j−1
X X j
WTl k ≤ θ j < WTl k
l=1 l=1

where j can range from 1 to N.

11
 n,l L
5. MCMC step: For λ = 1, . . . , Λ, construct a Markov chain {Yλ,T }
k l=0
n,0 n
with initial value Yλ,Tk
= Xλ,Tk
and stationary distribution
gλ (Yλ , ZTk )dPλ
R . (9)
gλ (Xλ,Tk , Zλ,Tk )dPλ
n n,L
6. Set Xλ,Tk
= Yλ,Tk
.
7. Set k = k + 1 and proceed to Step 1.
For a collection of N samples the effective sample size essΛ (Tk ) for Λ
targets is
N
essΛ (Tk ) = 2
1 + CΛ,k
where
v
u N Y Λ
1 ut1
X
n ,Z 2
CΛ,k = ( gλ (Xλ,T λ,Tk ) − WΛ,k )
WΛ,k N k
n=1 λ=1
N Λ
1 XY n
and WΛ,k = gλ (Xλ,T , Zλ,Tk ).
N k
n=1 λ=1

Finally, we mention the necessary modifications to compute conditional

expectation estimates of a function h(X1,Tk , . . . , XΛ,Tk ) = h(XTk ) that de-
pends, in general, on all targets.

E[h(XTk )|{ZTj }kj=1 ]

QΛ k−1
λ′ =1 gλ′ (Xλ′ ,Tk , Zλ′ ,Tk )p(Xλ′ ,Tk−1 |{Zλ′ ,Tj }j=1 )dXλ′ ,Tk−1 dPλ′
Z
= h(XTk ) R QΛ k−1
λ′ =1 gλ′ (Xλ′ ,Tk , Zλ′ ,Tk )p(Xλ′ ,Tk−1 |{Zλ′ ,Tj }j=1 )dXλ′ ,Tk−1 dPλ′
N Z QΛ n
1 X n λ′ =1 gλ′ (Xλ′ ,Tk , Zλ′ ,Tk )dPλ′
≈ h(XTk ) R QΛ .
N λ′ =1 gλ′ (Xλ′ ,Tk , Zλ′ ,Tk )dPλ′
n=1
In order to carry out the MCMC step in this case one has to construct a
Markov chain with stationary density
QΛ
gλ′ (Y, Zλ′ ,Tk )dPλ′
ρ(Y ) = R QΛ λ =1
′
(10)
λ′ =1 gλ′ (Xλ′ ,Tk , Zλ′ ,Tk )dPλ′
for each sample n = 1, . . . , N. The algorithm described above needs to be
modified only in the MCMC step where one needs to construct a Markov
chain with stationary density given by (10).

12
2 Numerical results
We present numerical results for multi-target tracking using the particle
filter with an MCMC step. We have synthesized tracks of targets moving on
the xy plane using a 2D near constant velocity model [1]. At each time t we
have a total of Kt targets and the evolution of the kth target (k = 1, . . . , Kt )
is given by

xk,t = Axk,t−1 + Bvk,t (11)

= [xk,t , ẋk,t , yk,t , ẏk,t]T ,

where (xk,t , ẋk,t ) and (yk,t , ẏk,t) are the xy position and velocity of the kth
target at time t. The matrices A and B are given by
   2 
1 T 0 0 T /2 0
 0 1 0 0   T 0 
A=  0 0 1 T  and B =  0
 
2
, (12)
T /2 
0 0 0 1 0 T

where T is the time between observations. For the experiments we have set
T = 1, i.e., noisy observations of the model are obtained at every step of the
model (11). The model noise vk,t is a collection of independent Gaussian
random variables with covariance matrix Σv defined as
 2 
σx 0
Σv = . (13)
0 σy2

In the experiments we have σx2 = σy2 = 1. Also, we have considered two

possible cases for the observation model, one linear and one nonlinear. Due
to the different possible combinations of targets to observations we use a
different index m to denote the obsevations. Since we do not assume any
clutter we have m = 1, . . . , Kt . If the mth observation zm,t at time t comes
from the kth target we have
 
xk,t
zm,t = + wm,t (14)
yk,t

for the linear observation model and

y
" #
arctan( xk,t )
zm,t = k,t
+ wm,t (15)
(x2k,t + yk,t
2 )1/2

13
for the nonlinear observation model. As is usual in the literature, the non-
linear observation model consists of the bearing θ and range r of a target.
The observation noise wm,t is white and Gaussian with covariance matrix
 2
σobs,x 0

Σw = 2 (16)
0 σobs,y

for the linear observation model and

 2 
σθ 0
Σw = (17)
0 σr2

for the nonlinear observation model. For the numerical experiments with
2
the linear observation model we chose σobs,x 2
= σobs,y = 1. For the numerical
experiments with the nonlinear observation model we chose σθ2 = 10−4 and
σr2 = 1. These values make our example comparable in difficulty to examples
appearing in the literature (see e.g. [16, 21, 22]).
The synthesized target tracks were created by specifying a certain sce-
nario, to be detailed below, of surviving, newborn and disappearing targets.
According to this scenario we evolved the appropriate number of targets
according to (11) and recorded the state of each target at each step. For
the surviving targets we created an observation by using the state of the
target in the observation model. Thus, for the linear observation model,
the observations were created directly in xy space by perturbing the xy
position of the target by (14). For the nonlinear observation model, the ob-
servations were created in bearing and range space θ, r by using (15). The
perturbed bearing and range were transformed to xy space by the trans-
formation x = r cos θ, y = r sin θ to create a position for the target in xy
space.
The newborn targets for the linear model were created in xy space di-
rectly by sampling uniformly in [−100, 100]. Afterwards, the observations of
the newborn targets were constructed by perturbing the x, y positions using
(14). The newborn targets for the nonlinear model were created in xy space
by sampling uniformly in [−100, 100]. Afterwards, we transformed the x, y
positions to the bearing and range space θ, r and perturbed the bearing and
range according to (15). The perturbed bearing and range were again trans-
formed back to xy space to create the position of the newborn target. Note
that both observation models do not involve the velocities. The newborn
target velocities were sampled uniformly in [−1, 1].
The number of targets at each observation instant is: K0 = 2, K1 = 2,
K2 = 1, K3 = 2, K4 = 3, K5 = K6 = . . . = K200 = 4. So, for the majority

14
of the steps we have 4 targets which makes the problem of tracking rather
difficult.

2.1 Implementation of the MCMC step

We present now the implementation of the extra (MCMC) step for the par-
ticle filter. First, we construct a continuous time model for which (11) is a
discrete approximation. This is not necessary. One can work with the dis-
crete model (11). However, we want to show that the proposed approach is
more general and can be applied to continuous time models as well. Second,
we show how the extra step leads to a sampling problem which we address
via MCMC sampling, in particular via Hybrid Monte Carlo (HMC).

2.1.1 Continuous time reformulation of the dynamic model

The 2D near constant velocity model (11) is a simplified discrete approxi-
mation of the formal linear SDE system

ẍt = σx u̇x,t (18)

ÿt = σy u̇y,t ,

where ux,t , uy,t are independent Brownian motions and σx u̇x,t , σy u̇y,t are
Gaussian white noises with covariances Rx (t) = σx2 δ(t) and Ry (t) = σy2 δ(t)
respectively and δ(t) is the delta function. We can define velocities in the x
and y directions, px,t = ẋt and py,t = ẏt and we can rewrite (18) as

ẋt = px,t,
ṗx,t = σx u̇x,t , (19)
ẏt = py,t ,
ṗy,t = σy u̇y,t .

The formal system (19) can be written rigorously

dxt = px,tdt,
dpx,t = σx dux,t , (20)
dyt = py,t dt,
dpy,t = σy duy,t .

The equations for the positions have zero noise and the equations for the
velocities have zero drift. We can write the SDE system (20) in matrix form

15
as      
dxt px,t 0 0  
 dpx,t   0
 dt +  σx 0  dux,t .
  
 dyt  =  py,t (21)
  
  0 0  duy,t
dpy,t 0 0 σy
Define the 4×1 vectors zt = [z1,t , . . . , z4,t ]T , a(zt ) = [a1,t , . . . , a4,t ]T , the 2×1
vector wt = [w1,t , w2,t ]T and a 4 × 2 constant matrix σ by
     
xt px,t   0 0
 px,t   0  ux,t  σx 0 
zt = 
 yt  , a(zt ) =  py,t  , wt = uy,t
   and σ =   0 0 .


py,t 0 0 σy

Note that      
z1,t xt z2,t
 z2,t   px,t   0 
 z3,t  =  yt and a(zt ) = 
 z4,t  .
    

z4,t py,t 0
With these definitions we can rewrite the SDE system (21) as

dzt = a(zt )dt + σdwt (22)

We suppose that we are at observation time T and we have N samples. After

the resampling step at observation time 0, we have obtained a collection of
pairs (z0n , zTn ), for n = 1, . . . , N. The system (22) can be solved explicitly. In
particular, we have for the position and velocity in the x direction
Z T
′ ′ ′
z1,T = z1,0 + z2,0 T + σx ux,s ds (23)
0
′ ′
z2,T = z2,0 + σx ux,T .

where we have used the fact that the Brownian motion ux,t is zero at t = 0.
′n = z n and z ′n = z n we have
For the nth sample, if we use z1,0 1,0 2,0 2,0
Z T
′n n n
z1,T = z1,0 + z2,0 T+ σx unx,s ds (24)
0
′n n
z2,T = z2,0 + σx unx,T ,

where unx,s , s ∈ [0, T ] is a new Brownian path for each sample.

We should note here that if the continuous time dynamic model is non-
linear and cannot be solved explicitly we can still use the proposed approach

16
by utilizing a numerical discretization of the dynamic model e.g. the Euler-
Maruyama or higher order schemes [9].
Since the targets evolve independently of one another and we are in-
terested in computing conditional expectation estimates of functions that
depend only on one target, we need to formulate the MCMC step only for
the case of a single target. The conditional density that we have to sam-
ple for the linear observation model is given by (see (9) in particle filter
algorithm with MCMC step for multiple targets)

′n ′n
gx (Z1,T , z1,T )gy (Z3,T , z3,T )dP
 ′n )2 
(Z1,T − z1,T
 ′n )2 
(Z3,T − z3,T
∝ exp − 2 exp − 2 dP (25)
2 ∗ σobs,x 2 ∗ σobs,y

where Z1,T the observation value of the x position of the target and Z3,T the
observation value of the y position of the target. Note that the x, y positions
′n , z ′n of the nth sample are determined by
z1,T 3,T

′n n n
z1,T ≈ z1,0 + z2,0 T + T σx unx,T (26)
′n n n
z3,T ≈ z3,0 + z4,0 T + T σy uny,T (27)

where we have used the fact that we go in one stepRof size T from one ob-
T
servation to the next to make the approximations 0 σx unx,s ds ≈ T σx unx,T
RT
and 0 σy uny,s ds ≈ T σy uny,T . Due to (26), (27) and since the initial con-
n , . . . , z n are fixed by the resampling step, the only quantities
ditions z1,0 4,0
that we need to sample with the conditional density (25) are the Brow-
nian values unx,T , uny,T . With the help of (26), (27) we can approximate
′n )g (Z
gx (Z1,T , z1,T ′n
y 3,T , z3,T )dP as

′n ′n
gx (Z1,T , z1,T )gy (Z3,T , z3,T )dP
  n − z n T − T σ un )2
(Z1,T − z1,0 2,0 x x,T
≈ exp − 2
2 ∗ σobs,x
n − z n T − T σ un )2
(Z3,T − z3,0 (unx,T )2 (uny,T )2

4,0 y y,T
+ 2 + + , (28)
2 ∗ σobs,y 2 2

where the last two terms in the exponential come from approximating the
′n )g (Z
Wiener measure dP. The expression (28) for gx (Z1,T , z1,T ′n
y 3,T , z3,T )dP is
n n
a Gaussian density for ux,T , uy,T . We do not need HMC to sample it. We can
rewrite it as the product of two Gaussian densities, one for unx,T and one for

17
uny,T and sample it with standard methods for Gaussian densities. However,
we show how HMC is implemented because for the nonlinear observation
model, the density to be sampled will no longer be Gaussian.

2.1.2 Hybrid Monte Carlo formulation

We present briefly the hybrid Monte Carlo (HMC) formulation that we have
used to sample the conditional density for each target. We start with the
linear observation model. Define the potential Vǫl (unx,T , uny,T ) by

n − z n T − T σ un )2
(Z1,T − z1,0 2,0 x x,T
Vǫl (unx,T , uny,T ) = 2
2 ∗ σobs,x
n − z n T − T σ un )2
(Z3,T − z3,0 4,0 y y,T
+ 2
2 ∗ σobs,y
(unx,T )2 (uny,T )2
+ + (29)
2 2
′n )g (Z
and the density gx (Z1,T , z1,T ′n
y 3,T , z3,T )dP becomes
 
′n ′n
gx (Z1,T , z1,T )gy (Z3,T , z3,T )dP ≈ exp −Vǫl (unx,T , uny,T ) .

Consider unx,T , uny,T as the position variables of a Hamiltonian system. We

define the 2D position vector q = [q1 , q2 ]T with q1 = unx,T and q2 = uny,T .
To each of the position variables we associate a momentum variable and we
write the Hamiltonian
pT p
Hǫl (q, p) = Vǫl (q) + ,
2
where p = [p1 , p2 ]T is the momentum vector. Thus, the momenta variables
are Gaussian distributed random variables with mean zero and variance 1.
The equations of motion for this Hamiltonian system are given by Hamilton’s
equations
dqi ∂Hǫl dpi ∂Hǫl
= and =− for i = 1, . . . , 2.
dτ ∂pi dτ ∂qi
HMC proceeds by assigning initial conditions to the momenta variables
T
(through sampling from exp(− p 2 p )), evolving the Hamiltonian system in
fictitious time τ for a given number of steps of size δτ and then using the

18
solution of the system to perform a Metropolis accept/reject step (more de-
tails in [10]). After the Metropolis step, the momenta values are discarded.
The most popular method for solving the Hamiltonian system, which is the
one we also used, is the Verlet leapfrog scheme. In our numerical imple-
mentation, we did not attempt to optimize the performance of the HMC
algorithm. For the sampling we used 100 Metropolis accept/reject steps
and 1 HMC step of size δτ = 10−1 to construct a trial path.
For the nonlinear observation model

gθ (Zθ,T , zT′n )gr (Zr,T , zT′n )

(Zθ,T − θ(zT′n ))2 (Zr,T − r(zT′n ))2
   
∝ exp − exp − (30)
2 ∗ σθ2 2 ∗ σr2
where Zθ,T , Zr,T are the bearing and range observation values for the target
and  ′n 
′n
z3,T ′n 2 ′n 2 1/2
θ(zT ) = arctan ′n and r(zT′n ) = (z1,T + z3,T )
z1,T
are the bearing and range values for the nth sample. We can use the same
procedure as in the linear observation model to define a Hamiltonian system
and its associated equations. We omit the details.

2.2 Linear observations

We start the presentation of our numerical experiments with results for the
linear observation model (14). Figures 1 and 2 show the evolution in the
xy space of the true targets, the observations as well as the estimates of the
improved particle filter. It is obvious from the figures that the improved
particle filter follows accurately the targets and there is no ambiguity in the
identification of the target tracks.
The performance of the improved particle filter with 100 samples is com-
pared to the performance of the generic particle filter with 120 samples in
Figure 3 by monitoring the evolution in time of the RMS error per target.
The RMS error per target (RMSE) is defined with reference to the true
target tracks by the formula
v
u
u 1 X Kt
RM SE(t) = t kxk,t − E[xk,t |Z1 , . . . , Zt ]k2 (31)
Kt
k=1

where k·k is the norm of the position and velocity vector. Note that the state
vector norm involves both positions and velocities even though the obser-
vations use information only from the positions of a target. xk,t is the true

19
 0

-1000
Y

-2000

-3000

-4000
-2000 -1000 0 1000 2000 3000
X

Figure 1: Linear observation model. The solid lines denote the true target
tracks, the crosses denote the observations and the dots the conditional
expectation estimates from the improved particle filter. We have plotted the
conditional expectation estimates every 5 observations to avoid cluttering
in the figure.

250

200

150
Y

100

50

0 50 100 150 200

X

Figure 2: Linear observation model. Detail of Figure 1.

20
 PF with MCMC step nsamples=100
2000 PF nsamples=120

RMS Error per target

1500

1000

500

0
50 100 150 200
Time steps

Figure 3: Linear observation model. Comparison of RMS error per target

for the improved particle filter and the generic particle filter.

state vector for target k. E[xk,t |Z1 , . . . , Zt ] is the conditional expectation

estimate calculated with the improved or generic particle filter depending
on whose filter’s performance we want to calculate.
The improved particle filter has a computational overhead of the order of
a few percent compared to the generic particle filter. We have thus used the
generic particle filter with more samples than the improved particle filter.
This additional number of samples more than accounts for the computational
overhead of the improved particle filter. As can be seen in Figure 3 the
generic particle filter’s accuracy deteriorates quickly. On the other hand,
the improved particle filter maintains an O(1) RMS error per target for the
entire tracking interval. The average value of the RMS error over the entire
time interval of tracking is about 2.5 with standard deviation of about 0.5.
For the generic particle filter, the average of the RMS error over the time
interval of tracking is about 800 with standard deviation of about 760.
Figure 4 compares the effective sample size for the generic particle filter
and the improved particle filter. Because the number of samples is different
for the two filters we have plotted the effective sample size as a percentage
of the number of samples. We have to note that, after about 50 steps,
the generic particle filter started producing observation weights (before the
normalization) which were numerically zero. This makes the normalization
impossible. In order to allow the generic particle filter to continue we chose

21
 1

Effective sample size as percentage of number of samples

PF with MCMC step nsamples=100
PF nsamples=120
0.8

0.6

0.4

0.2

0
0 50 100 150 200
Time steps

Figure 4: Linear observation model. Comparison of effective sample size for

the improved particle filter and the generic particle filter.

at random one of the samples, since all of them are equally bad, and assigned
all the weight to this sample. We did that for all the steps for which the
observation weights were zero before the normalization. As a result, the
effective sample size for the generic particle filter drops down to 1 sample
after about 50 steps. Once the generic particle filter deviates from the true
target tracks there is no mechanism to correct it. Also, we tried assigning
equal weights to all the samples when the observation weights dropped to
zero. This did not improve the generic particle’s performance either. On
the other hand, the improved particle filter maintains an effective sample
size which is about 25% of the number of samples.

2.3 Nonlinear observations

We continue with results for the nonlinear observation model (15). Figures 5
and 6 show the evolution in the xy space of the true targets, the observations
as well as the estimates of the improved particle filter. Again, as in the
case of the linear observation model, the improved particle filter follows
accurately the targets and there is no ambiguity in the identification of the
target tracks.
The case of the nonlinear observation model is much more difficult than
the case of the linear observation model. The reason is that for the nonlinear

22
 0

-1000
Y

-2000

-3000

-4000
-2000 -1000 0 1000 2000 3000
X

Figure 5: Noninear observation model. The solid lines denote the true target
tracks, the crosses denote the observations and the dots the conditional
expectation estimates from the improved particle filter. We have plotted the
conditional expectation estimates every 5 observations to avoid cluttering
in the figure.

200

100
Y

-100
-100 0 100 200
X

Figure 6: Nonlinear observation model. Detail of Figure 5.

23
observation model, the observation errors, though constant in bearing and
range space, they become position dependent in xy space. In particular,
when x and/or y are large, the observation errors can become rather large.
This is easy to see by Taylor expanding the nonlinear transformation from
bearing and range space to xy space around the true target values. Suppose
that the true target bearing and range are θ0 , r0 and its xy space position
is x0 = r0 cos θ0 , y0 = r0 sin θ0 . Also, assume that the observation error in
bearing and range space is, respectively, δθ and δr. The xy position of a
target that is perturbed by δθ and δr in bearing and range space is (to first
order)

x = x0 − y0 δθ − δr cos θ0
y = y0 + x0 δθ − δr sin θ0 .

Thus, the perturbation in xy space can be significant even if δθ and δr are

small. In our example we have σθ = 10−2 . So, when the true target x and
y values become of the order of 103 as happens for some of the targets, the
observation value in bearing and range space can be quite misleading as far
as the xy space position of the target is concerned. As a result, even if one
does a good job in following the observation in bearing and range space, the
conditional expectation estimate of the xy space position can be inaccurate.
With this in mind, we have used 200 samples for the improved particle
filter and 220 samples for the generic particle filter. Again, the extra samples
used for the generic particle filter more than account for the computational
overhead of the improved particle filter. The performance of the improved
particle filter is compared to the performance of the generic particle filter in
Figure 7 by monitoring the evolution in time of the RMS error per target.
The generic particle filter’s accuracy again deteriorates rather quickly. The
error for the improved particle filter is larger than in the linear observation
model but never exceeds about 80 even after 200 steps when the targets have
reached large values of x and/or y. The average value of the RMS error over
the entire time interval of tracking is about 22 with standard deviation of
about 21. For the generic particle filter, the average of the RMS error over
the time interval of tracking is about 760 with standard deviation of about
770.
Figure 8 compares the effective sample size for the generic particle filter
and the improved particle filter. After about 60 steps, the generic parti-
cle filter, started producing observation weights (before the normalization)
which were numerically zero. This makes the normalization impossible. In
order to allow the generic particle filter to continue we chose at random

24
 PF with MCMC step nsamples=200
PF nsamples=220
2000
RMS Error per target

1500

1000

500

0
50 100 150 200
Time steps

Figure 7: Nonlinear observation model. Comparison of RMS error per target

for the improved particle filter and the generic particle filter.

1
Effective sample size as percentage of number of samples

PF with MCMC step nsamples=200

PF nsamples=220
0.8

0.6

0.4

0.2

0
0 50 100 150 200
Time steps

Figure 8: Nonlinear observation model. Comparison of effective sample size

for the improved particle filter and the generic particle filter.

25
one of the samples, since all of them are equally bad, and assigned all the
weight to this sample. We did that for all the steps for which the observation
weights were zero before the normalization. As a result, the effective sample
size for the generic particle filter drops down to 1 sample after about 60
steps. Once the generic particle filter deviates from the true target tracks
there is no mechanism to correct it. Also, we tried assigning equal weights
to all the samples when the observation weights dropped to zero. This did
not improve the generic particle’s performance either. On the other hand,
the improved particle filter maintains an effective sample size which is about
25% of the number of samples.

3 Discussion
We have presented an algorithm for multi-target tracking which builds on the
existing particle filter methodology for multi-target tracking by appending
an MCMC step after the particle filter resampling step. The purpose of the
addition of the MCMC step is to bring the samples closer to the observation.
Even though the addition of an MCMC step for a particle filter has been
proposed and used before [7], to the best of our knowledge, the currently
proposed implementation of the MCMC step is novel (see also [23] for a
related approach).
We have tested the performance of the algorithm on the problem of
tracking multiple targets evolving under the near constant velocity model [1].
We have examined two cases of observation models: i) a linear observations
model involving the positions of the targets and ii) a nonlinear observation
model involving the bearing and range of the targets. For both cases the
proposed improved particle filter exhibited a significantly better performance
than the generic particle filter. Since the improved particle filter requires
more computations than the generic particle filter it is bound to be more
expensive. However, the computational overhead of the improved particle
filter is rather small, of the order of a few extra samples worth for the generic
particle filter.
In [15] we proposed another way of performing the extra MCMC step
of a particle filter. That approach was based on modifying the drift of the
dynamic model and then accounting for the modification via Girsanov’s the-
orem. In the current work, we use the original drift of the dynamic model
without any modification. For the case of multi-target tracking with obser-
vations at every time step both algorithms perform equally well. Thus, at
first sight it would appear that there is no need for the extra complication

26
of modifying the drift of the dynamic model and then accounting for the
modification by Girsanov’s theorem as was done in [15]. However, in cases
where there are only sparse observations, the sampling of the conditional
density needed for the extra step can be much more difficult (and conse-
quently expensive) for the original dynamic model than for the modified
dynamic model. With this in mind, the approach in [15] seems to have
wider applicability. A detailed comparison of the algorithm proposed in
the current work and the one proposed in [15] will be presented in a future
publication.

Acknowledgements
We are grateful to Prof. J. Weare for many discussions. Also, we would like
to thank the Institute for Mathematics and its Applications in the University
of Minnesota for its support.

References
[1] Bar-Shalom Y. and Fortmann T.E., Tracking and Data Association,
Academic Press, 1988.

[2] Bar-Shalom Y. and Blair W.D., Eds. Multitarget-Multisensor Track-

ing: Applications and Advances, vol. III, Norwood, MA, Artech House,
2000.

[3] Blackman S. and Popoli R., Design and Analysis of Modern Tracking
Systems, Norwood, MA, Artech House, 1999.

[4] Chorin, A.J. and Tu X., Implicit sampling for particle filters, Proc. Nat.
Acad. Sc. USA 106 (2009) pp. 17249-17254.

[5] Doucet A., de Freitas N. and Gordon N. (eds.) , Sequential Monte Carlo
Methods in Practice, Springer NY, 2001.

[6] Fortmann T. E., Bar-Shalom Y. and Scheffe M., Sonar tracking of mul-
tiple targets using joint probabilistic data association, IEEE J. Ocea.
Eng., vol.8 (1983) pp.173-184.

[7] Gilks W. and Berzuini C., Following a moving target. Monte Carlo
inference for dynamic Bayesian models, J. Royal Stat. Soc. B 63 (1)
(1999) pp. 2124-2137.

27
 [8] Gordon N.J., Salmond D.J. and Smith A.F.M., Novel approach to
nonlinear/non-Gaussian Bayesian state estimation, Proc. Inst. Elect.
Eng. F 140(2) (1993) pp. 107-113.

[9] Kloeden P. and Platen E., Numerical Solution of Stochastic Differential

Equations, Springer Berlin 1992.

[10] Liu J.S., Monte Carlo Strategies in Scientific Computing, Springer NY,
2001.

[11] Liu J.S. and Chen R. Sequential Monte Carlo Methods for Dynamic
Systems. Journal of the American Statistical Association, vol.93 no.
443 (1993) pp. 1032-1044.

[12] Mahler R. P., Statistical Multisource-Multitarget Information Fusion,

Artech House Publishers MA, 2007.

[13] Mahler R. P. S., Multitarget Bayes filtering via first-order multitarget

moments, IEEE Trans. Aero. Elect. Sys., Vol. 39 no. 4 (2003) pp. 1152-
1178.

[14] Mahler R.P.S. and Maroulas V. Tracking Spawning Objects, (2010)

submitted.

[15] Maroulas V. and Stinis P., A Girsanov Monte Carlo approach to particle
filtering for multi-target tracking, (2010), arXiv:1006.3100v1 .

[16] Ng W., Li J.F., Godsill S.J. and Vermaak J., A hybrid approach for
online joint detection and tracking for multiple targets, Proc. IEEE
Aerospace Conference (2005).

[17] Oksendal B., Stochastic Differential Equations, An Introduction with

Applications, Sixth Edition Springer Heidelberg 2005.

[18] Snyder C., Bengtsson T., Bickel P. and Anderson J., Obstacles to High-
dimensional Particle Filtering, Mon. Wea. Rev., Vol. 136 (2008) pp.
4629-4640.

[19] Stinis P., Conditional path sampling for stochastic differential equations
by drift relaxation, (2010), arXiv:1006.2492v1.

[20] Stoer J, and Bulirsch R., Introduction to Numerical Analysis, Third

Edition, Springer 2002.

28
[21] Vermaak J., Godsill S. and Perez P., Monte Carlo filtering for multi-
target tracking and data association, IEEE Trans. Aero. Elect. Sys.,
41(1) (2005) pp. 309-332.

[22] Vo B-N., Singh S. and Doucet A., Sequential Monte Carlo Methods for
Multi-Target Filtering with Random Finite Sets, IEEE Trans. Aero.
Elect. Sys., 41(4) (2005) pp. 1224-1245.

[23] Weare J., Particle filtering with path sampling and an application to a
bimodal ocean current model, J. Comp. Phys. 228 (2009) pp. 4312-4331.

29