The Pairwise-Markov Bernoulli Filter

The Bernoulli filter is a general, Bayes-optimal solution for tracking a single disappearing and reappearing target, using a sensor whose observations are corrupted by missed detections and a general, known clutter process. Like virtually all target-tracking algorithms it presumes restrictive independence assumptions, namely a hidden Markov model (HMM) structure on the sensor and target. That is, the current state of the target depends only on its previous state, and the measurement collected from it depends only on its current state. Pieczynski’s pairwise Markov model (PMM) relaxes these restrictions. In it, the current target state can additionally depend on the previous measurement; and the current measurement can additionally depend on the previous measurement and previous target state. In this paper we show how to correctly generalize the PMM to the multitarget (MPMM) case; and use the MPMM to derive a “PMM Bernoulli filter” that obeys PMM rather than restrictive HMM sensor/target statistics.


I. INTRODUCTION
The Bernoulli filter was independently and contemporaneously devised by Vo [1] and Mahler [2,Sec. 14.7]. It is a general and Bayes-optimal solution for tracking a single disappearing and reappearing target, using a sensor whose observations are corrupted by missed detections and a general, known clutter process. It propagates a probability hypothesis density (PHD) D via time-update and measurement-update steps D(x k−1 |Z 1:k−1 ) → D(x k |Z 1:k−1 ) and D(x k |Z 1:k−1 ) → D(x k |Z 1:k ), where Z 1:k : Z 1 ,. . . , Z k is the time-sequence of collected measurement-sets. See Section VI-A for more detail.
Like virtually all target-tracking algorithms, the Bernoulli filter presumes restrictive independence assumptions, namely a hidden Markov model (HMM) structure on the sensor and target. That is, at time t k the target's state x k depends only on its previous state x k−1 with Markov transition density f (x k |x k−1 ); and the measurement y k that the sensor collects from it depends only on x k with measurement density f (y k |x k ). Pieczynski's pairwise Markov model (PMM) [3]- [7] relaxes these restrictions.

A. THE PAIRWISE MARKOV MODEL (PMM)
The PMM generalizes the HMM by treating the target and sensor as a joint dynamical system with joint state (x k , y k ), The associate editor coordinating the review of this manuscript and approving it for publication was Qiangqiang Yuan. which is governed by a Markov transition density f (x k , y k |x k−1 , y k−1 ) = f (x k |x k−1 , y k−1 ) · f (y k |x k , x k−1 , y k−1 ) (1) where the factorization on the right is due to Bayes' rule. In the PMM, the current target state can additionally depend on the previous measurement (as described by f (x k |x k−1 , y k−1 ) (i.e., the target can be non-Markovian); and in that the current measurement can additionally depend on the previous measurement and the previous target state as described by f (y k |x k , x k−1 , y k−1 ) (and thus measurement noise can be colored or correlated with plant noise [7, p. 4487]). See Section III for more detail.
Pieczynski and Desbouvries [6] have described practical Kalman filter-based implementations of PMMs to singletarget tracking. Petetin and Desbouvries [7] proposed a PMM generalization of the probability hypothesis density (PHD) filter of [2,Sec. 16.3]; described concrete practical applications and implementations; and demonstrated that their The PMM Bernoullii filter can be summarized as follows. Let us be given: (i) κ k (Y k ) (the multi-object probability density function of the clutter process); (ii) p S (x k−1 ) (the probability that the target will not disappear at time t k−1 ); (iii) q B k (the probability that the target will reappear at time t k− after having disappeared); (iv) s B k (x k ) (the target's spatial density after reappearance); (v) p D (x k ) (the target's probability of detection); (vi) f (x k , y k |x k−1 , y k−1 ) (the PMM transition density); (vii) f (x k |x k−1 , y k−1 ) (the marginal of f (x k , y k |x k−1 , y k−1 ); (viii) M x k (x k−1 ) = f (x k |x k−1 ) (Markov density associated with transition ({x k−1 }}, Ø) → ({x k }, Ø), see (76); and (ix) L y k (x k ) = f (y k |x k ) (measurement density associated with (Ø, Ø) → ({x k },{y k }), see (89). Define (2) where by convention the summation vanishes if Y k = Ø. Also, if f (x k |x k−1 , y k−1 ) is the marginal of f (x k , y k |x k−1 , if Y k−1 = Ø (and where by convention the second summation vanishes if Y k = Ø); whereas if Y k−1 = Ø, if Y k−1 = Ø whereas if Y k−1 = Ø (and employing the notation defined in (47)), Also, if f (x) is a density function and 0 ≤ h(x) ≤ 1 a unitless function, define Given this, the PMM Bernoulli filter is given by the following single-step recursive update, (10), as shown at the bottom of the next page. This equation is derived in Appendix B. If the PMM is actually an HMM, then (see (36)) from which it follows that in which case, as will be shown in (200), (10) reduces to the single-step HMM Bernoulli filter as given in (111).

D. ORGANIZATION OF THE PAPER
The remainder of the paper is organized as follows: A brief summary of the mathematical theory required to understand the paper (Section II); the PMM (Section III); the MPMM (Section IV); the corrected MPMM transition density (Section V); and the Bernoulli MPMM filter (Section VI). Conclusions can be found in Section VII, and the Bernoulli MPMM and PMM Bernoulli filters are derived in Appendices A and B, respectively. The following notation will be employed hereafter: A : = B means ''A is defined to be B''; and A != B means ''A is an abbreviation of B.''

II. OVERVIEW OF FINITE-SET STATISTICS (FISST)
This section summarizes the theory necessary to understand the remainder of the paper. Greater detail can be found in books [2], [10]- [12], tutorials [13]- [15], and a short survey of advances c. 2015 [16]. Also, systematic investigations of FISST vs. ''point processes'' can be found in [17], [18] and of FISST vs. measurement-to-track approaches in [19]. Significant recent advances can be found in [20], [21]. Specifically, [20] describes an implementation of the generalized labeled multi-Bernoulli (GLMB) filter that is capable of simultaneously tracking over a million 2D targets in significant clutter in real time using off-the-shelf computing equipment, as well as a theoretically rigorous, large-scale track quality measure, ''OSPA (2) ''; and [21] describes a multiscan extension of the GLMB filter.
The section is organized as follows: random finite sets (Section II-A); multitarget calculus (Section II-B); Bernoulli RFSs (Section II-C); and the multitarget recursive Bayes filter (Section II-D).

A. RANDOM FINITE SETS (RFSs)
Let be a single-target state-space with x ∈ , and let be the sensor measurement space with z ∈ . Then the state of a multitarget system is represented as a finite subset X = {x 1 ,. . . , x n } ⊆ with X = Ø for n = 0. The number of elements in X is denoted as |X |. In a Bayesian approach, unknown states are random variables. Thus an unknown multitarget state is a random finite set (RFS) ⊆ .

B. MULTITARGET CALCULUS
A multitarget density function is a function f (X ) ≥ 0 of the finite-set variable X ⊆ such that the unit of measurement of (13) where f n (x 1 ,. . . , x n ) : = f ({x 1 ,. . . , x n })/n! for distinct x 1 ,. . . , x n and f n (x 1 ,. . . , x n ) : = 0 otherwise. Every random finite state-set has a multitarget probability distribution where ι is the unit of measurement of . An MPMM density function is a joint probability density if ∫ f (X , Y )δX δY = 1. If ⊆ and ⊆ are RFSs then , have an MPMM probability density f , (X , Y ) that describes the statistical correlation between them.
The probability generating functional (p.g.fl.) of is, for where h X = 1 if X = Ø and h X = x∈X h(x) otherwise. The simplest nontrivial p.g.fl.'s are where s(x) ≥ 0 is a probability density function on . If 0 ≤ g(z) ≤ 1 for z ∈ then the joint p.g.fl. of , is The intuitive definition of the Volterra functional deriva- where δ x (y) is the Dirac delta function concentrated at x.
(For a rigorous definition see [14].) If X = {x 1 ,. . . , x n } with |X | = n then the iterated functional derivative is δG δX where the '' • '' notation indicates that δ/δX is taken with respect to the variable h and δ/δY with respect to the variable g. When Y = Ø or X = Ø we have: .
(10) VOLUME 8, 2020 The p.g.fl. and distibution of an RFS are related by: Likewise, the bivariate p.g.fl. and bivariate multitarget distibution of RFSs , are related by: C. THE BERNOULLI RFS An RFS of special importance for this paper, the Bernoulli RFS, is most easily described using its p.g.fl. : where 0 ≤ q ≤ 1 and the probability density s(x) are, respectively, the existence probability and spatial distribution of a single target.

D. MULTITARGET RECURSIVE BAYES FILTER
Given a time-sequence Z 1:k :Z 1 , . . . , Z k of collected measurement-sets from a sensor, this is: and where f (X k |X k−1 , Z 1:k−1 ) is the multitarget statetransition density and f (Z k |X k , Z 1:k−1 ) is the sensor multitarget measurement density. It is assumed that In this paper we will be concerned with f (X k |X k−1 ) and f (Z k |X k ) for only the ''standard'' multitarget motion and measurement models, respectively-see Section V-A.

III. THE PAIRWISE MARKOV MODEL (PMM)
The section is organized as follows: single-target recursive Bayes filter (Section III-A); the PMM (Section III-B); and single-target tracking using PMMs (Section III-C).

A. SINGLE-TARGET RECURSIVE BAYES FILTER
The PMM concept is most easily explained via the singletarget recursive Bayes filter. Let x ∈ denote a single-target state and z ∈ a single-target measurement. In the Bayesian approach, the unknown state at time t k is a random variable X k|k ∈ and the measurement process at time t k is a random variable Z k ∈ . Let us be given: 1. the distribution f (x 0 ) of the initial state X 0|0 ; 2. the sequence z 1:k : z 1 ,. . . , z k of measurements collected from the target at times t 1 ,. . . , t k ; 3. the transition density f (x k |x k−1 , z 1:k−1 ), which describes the evolution of X k−1|k−1 at time t k−1 to X k|k−1 at time t k ; and 4. the measurement density f (z k |x k , z 1:k−1 ) at time t k , which characterizes the statistics of Z k if x k is a realization of X k|k−1 . Given this, the recursive Bayes filter is defined by the timeand measurement-update equations It is usually assumed that f (x k |x k−1 , z 1:

B. THE PAIRWISE MARKOV MODEL
Now let the state space be the Cartesian product × rather than . In this case, the unknown quantity at time t k is the joint state of the joint target-measurement system, and is represented as the random pair (X k|k , Y k ) ∈ × . What is unknown is not only X k|k and Y k but also their statistical correlation, as described by the posterior PMM density f (x k , y k |z 1:k−1 ). Let us be given 1. a PMM transition density that describes the evolution of the PMM system; 2. the measurement density f (z k |x k , y k ) = δ y k (z k ) of the PMM system, in which case it follows that the measurement equation is z k = η(x k , y k ) with measurement function η(x, y) : = y-i.e., if the joint system has state (x, y) then y is the only measurement that can be collected from it. Given this, f (x k , y k |z 1:k−1 ) can be recursively derived from f (x k−1 , y k−1 |z 1:k−2 ) as (29)-(31), shown at the bottom of the next page. Here, (29) is the Bayes' filter time-update step for the PMM; (30) incorporates the Bayes' filter measurementupdate step for the PMM; and (31) follows from the fact is an initial target distribution and f (y 1 |x 1 ) is an initial measurement density.
Remark 2: Since Z k = Y k , in the PMM literature the distinction between z j (a collected measurement) and y j (a realization of the unknown random variable Y j ) is notationally suppressed: The estimated measurement density at time t k can be determined from f (x k , y k |z 1:k−1 ) via A PMM reduces to an HMM if, for k > 1, in which case (28) reduces to: Thus PMMs significantly weaken HMM's to encompass non-Markov targets and correlated sensor noise [7, p. 4487].

C. SINGLE-TARGET TRACKING USING PMMs
The single-target posterior distribution f (x k |y 1:k ) can be recursively propagated as follows [7,Eq. 12]: Note that f (x k |y 1:k ) is related to f (x k , y k |y 1:k−1 ) as follows: so that is the MAP estimate of the target state given y 1:k . The predicted target distribution is

IV. MULTITARGET PMM (MPMM)
This is a direct generalization of the single-target PMM. Let Fin( ) denote the set of multitarget states (i.e., the finite subsets X of ) and let Fin( ) denote the set of multitarget measurements (i.e., the finite subsets Z of ). Then the unknown multitarget state at time t k is an RFS k|k ⊆ and the multitarget measurement process is an RFS k ⊆ . Now let the state space be Fin( )×Fin( ). Then the unknown state at time t k is that of the joint multitarget, multimeasurement system, as represented by the random pair ( k|k , k ) ∈ Fin( )×Fin( ). What is unknown is not only k|k and k but also their statistical correlation, as described by the posterior MPMM density f (X k , Y k |Y 1:k−1 ). Let us be given: 1. an MPMM transition density describing the evolution of the MPMM system; Then as with the single-target case, the recursions for the MPMM density f (X k , Y k |Y 1:k−1 ) and the multitarget posterior f (X k |Y 1:k ) are, respectively, as in the next section.

V. MPMM TRANSITION DENSITIES
The section is organized as follows: the ''standard'' multitarget motion and measurement models (Section V-A); the original MPMM transition model (Section V-B); the corrected MPMM transition model (Section V-C); and the evolution of the MPMM pair ({x k−1 },{y k−1 }) according to this model (Section V-D). The remaining subsections address extensions of this basic evolution model: the general evolution of ({x k−1 },{y k−1 }) (Section V-E); the evolution of ({x k−1 }, Ø) (Section V-F); the evolution of (Ø, Ø) (Section V-G); and the evolution of (Ø,{y k−1 }) (Section V-H).
This question was addressed in [9] by endeavoring to infer the form of its p.g.fl.
The ''standard'' multitarget motion model [2, Eq. 14.273], [10, Eq. 5.94], presumes that: (a) individual target motions are statistically independent; (b) the probability that a target with state x k−1 at time t k−1 will survive to time t k is p S (x k−1 ); (c) if so, then f (x k |x k−1 ) is the probability (density) that it will transition to a target with state x k ; and (d) f B (X k ) is the multitarget density of newly-appearing targets, with corresponding p.g.fl. G B [h k ]. This motion model is used to construct the ''standard'' multitarget Markov density f (X k |X k−1 ), with corresponding p.g.fl.
The ''standard'' multitarget measurement model [2, Eq. 14.290], [10, Eq. 5.104] presumes that: (e) all measurements are generated independently of each other; (f) the probability that a target with state x k at time t k will generate a measurement is p D (x k ); (g) if so, then f (y k |x k ) is the probability (density) that the measurement is y k ; and (h) f κ (Z k ) is the multi-object density of the clutter process, with corresponding p.g.fl. G κ [g k ]. This model is used to construct the standard multitarget measurement density f (Y k |X k ), with corresponding p.g.fl.

B. THE ORIGINAL MPMM TRANSITION MODEL
The multitarget analog of (36) is where G κ [g k ] characterizes the clutter process; G B [h k ] characterizes the target-appearance process; and where the evolution model is Equation (53)  In retrospect, (53) cannot be correct for at least two reasons. First, when the underlying PMM is an HMM-i.e., when Thus from (52-54) we see that rather than, as should be the case, 1 Second, suppose that the scenario contains at most a single target obscured by clutter. Then the multi-object system will always be described by an MPMM pair (X k , Y k ) with |X k | ≤ 1. In this case the evolution of the system must consist of transitions . . .
Such an evolution is impossible if it is given as (53) since describes a system that has as many as |Y k−1 | targets.
Accordingly, a corrected model is required, as follows.
Then replace (53) with In this case does properly reduce to (51). Thus the incorrect (58) is replaced with Equation (62) has a physical interpretation. Suppose that p S = 1 (the target does not disappear) and p D = 1 (the target is detected) and |Y k−1 | > 0. Then (62) reduces tõ The corresponding distribution is: 1. If x k−1 evolves to x k and measurement y k is collected from x k then the probability (density) that this event will occur is: 2. If x k−1 evolves to x k but x k is not detected, then the probability (density) that this event will occur is, if f (x k |x k−1 , 3. If x k−1 does not survive to time t k then no (nonempty) measurement can be collected from it and so the probability that this event will occur is: 4. If x k−1 does not survive to time t k and yet measurement y k is collected from it, this is an impossibility and so the probability that this event will occur is: Subsections V-E through V-H will address generalizations and extensions of this basic evolution model.
The dynamics model f (X k , Y k |{x k−1 },{y k−1 }) of the previous section is a special case of a more general model, deduced from an arbitrary MPMM density f (X k , Y k |{x k−1 },{y k−1 }) assuming only that |X k |, |Y k | ≤ 1. Since a nonexistent target cannot generate a nonempty measurement, we may assume that f (Ø, f (x k , y k |x k−1 , y k−1 ) VOLUME 8, 2020 This reduces to the model of Remark 3: To simplify notation, this is what will be assumed later in Section VI-C (though this assumption is not a necessity).
Remark 4: To simplify notation, it will be assumed in Section VI-C that p D (x k |x k−1 , Ø) does not depend on x k−1 and that f (y k |x k , x k−1 , Ø) does not depend on Then it is easily shown that For future reference, the p.g.fl. of f (X k , Y k |{x k−1 }, Ø) is Remark 5: To simplify notation, it will later be assumed in Then it is easily shown that Here,q B k is the ''birth'' probability-i.e., the nonexistent target Ø at time t k−1 transitions to a target with state x k at time t k -and s B k (x k ) is its spatial distribution. Remark 6: The obvious choice for q B k s B k (x k ) is the multiobject version of (40): For future reference, the p.g.fl. of f (X k , Y k |Ø, Ø) is Remark 7: To simplify notation, it will later be assumed in Section VI-C that q B k (Ø, y k−1 ) = q B k and s B k (x k |Ø, y k−1 ) = s B k (x k ) and p D (x k |Ø, y k−1 ) = p D (x k ) and f (y k |x k , Ø, y k−1 ) = f (y k |x k ). As in Section V-G it follows that and that the corresponding p.g.fl. is from which we conclude that

VI. THE BERNOULLI MPMM FILTER
The section is organized as follows: the Bernoulli filter (Section VI-A); the Bernoulli MPMM filter (Section VI-B); transition p.g.fl.'s for the Bernoulli MPMM filter (Section VI-C); summary of the Bernoulli MPMM filter (Section VI-D); derivation of the Bernoulli MPMM filter update when Y k−1 = Ø (Section VI-E); and derivation of the Bernoulli MPMM filter update when Y k−1 = Ø (Section VI-F).

A. THE BERNOULLI FILTER
The Bernoulli filter [1], [2,Sec. 14.7] is the special case of the multitarget Bayes filter when at most a single target is present-i.e., when where Z k was defined in (2). Note that, as presented in [2, Sec. 14.7], the Bernoulli filter propagates two items, not one: the probability of target existence p k|k = ∫ f ({x k }|Z 1:k )dx k and the target spatial distribution f k|k (x k ) = f ({x k }|Z 1:k )/p k|k . But it is clear that the filter in (109,110) differs from that in [2, Sec. 14.7] only by a change of notation (although the former is significantly simpler in form). A tutorial on the (original) Bernoulli filter can be found in [22].
State estimation using D k|k is as in [2]. The target exists if p k|k > τ for some threshold τ > 1/2 ; and if it exists, its state is the MAP estimate argsup x D k|k (x).
Eqs. (109,110) can be consolidated by substitution into the single-step update equation, (111), as shown at the bottom of the next page, where M Z k (x k−1 ) was defined in (8).

B. THE BERNOULLI MPMM FILTER
This is a special case of the MPMM filter when |X k−1 |, |X k | ≤ 1 for all k ≥ 1. Section V-G described a simplified target-appearance model. This model will allow us to avoid the factor G B in (61) by assuming that G B [h k ] = 1 identically.
Remark 8: Note that this simplified target-appearance model would not be acceptable in the multitarget case, since then the number of targets could never increase.
Given that G B [h k ] = 1, the p.g.fl. (61) of the MPMM transition density reduces to: where either X k−1 = Ø or X k−1 = {x k−1 } for all k ≥ 1. From (44), the p.g.fl. update for the Bernoulli MPMM filter is where the numerator is and the denominator is 2. If Z k−1 = Ø is collected then: These equations are derived in Appendix A. Remark 9: In regard to (129), consider the following special case: κ k (Y k ) = 0 identically (no clutter); p S (x k−1 ) = 1 (target never disappears); and p D (x k ) = 1 (perfect detection); in which case |X k | = |Y k | = 1 for all k ≥ 1. Then Eq. (129) should reduce to (31)-which is indeed the case.

VII. CONCLUSION
The Bernoulli filter is a general solution for tracking a single disappearing and reappearing target, using a sensor whose observations are corrupted by missed detections and a general, known clutter process. The Bernoulli filter presumes restrictive independence assumptions, namely a hidden Markov model (HMM) structure. That is, the current target state depends only on the previous target state; and the measurement that it generates depends only on its current state.
Pieczynski's pairwise Markov model (PMM) relaxes these restrictions. In it, the current target state can additionally depend on the previous measurement; and the current measurement can additionally depend on the previous measurement and the previous target state.
In this paper we: (i) generalized PMMs to the multitarget case (MPMM); (ii) devised a theoretically rigorous formula for the ''standard'' MPMM transition density (see (60,61)); (iii) derived transition models for the elementary MPMM pairs (X k , Y k ) with |X k |, |Y k | ≤ 1 (Sections V-D through V-H); (iv) used them to derive the Bernoulli MPMM filter (an MPMM generalization of the Bernoulli filter, Section VI); and then used it to derive the PMM Bernoulli filter (a generalization of the Bernoulli filter that obeys PMM rather than HMM sensor and target statistics.
Future work will be devoted to generalization of the PMM Bernoulli filter to multiple correlated sensors.

APPENDIX A DERIVATION OF THE BERNOULLI MPMM FILTER
The derivation has two parts: when Y k−1 = Ø (Appendix A.1) and when when Y k−1 = Ø (Appendix A.2).

1) DERIVATION OF THE BERNOULLI MPMM FILTER UPDATE WHEN Y k−1 = Ø
Let us turn to the derivation of the Bernoulli MPMM filter update when Y k−1 = Ø. Let Then from (114), (116), and (117), For fixed h k , abbreviate in which case (134) VOLUME 8, 2020 For W ⊆ Y k , note that and where L y (x) != f (y|x). Thus from the product rule for functional derivatives [2, p. 389], and so substituting g k = 0 and using the fact that we get where and Thus after substitution and collection of like terms we get: δG Consequently, and as claimed, and with normalization factor and where, by (118), First note that, setting h k = 0, and so Taking δ/δY k of both sides with respect to g k and then setting g k = 0 we get: Now note that and so Thus For the first term in this sum, note that where, for W ⊆ Y k , if otherwise (168) VOLUME 8, 2020 and so where For the second term in (166), note that and where for W ⊆ Y k , Thus and so Thus setting g k = 0 in (166) and substituting (173) and (182), we get (129), (183)-(185), as shown at the bottom of the previous page, where

APPENDIX B DERIVATION OF THE PMM BERNOULLI FILTER
The derivation has two parts: when Z k−1 = Ø (Appendix B.1) and when Z k−1 = Ø (Appendix B.2).
2) DERIVATION OF THE PMM BERNOULLI FILTER WHEN Y k−1 = Ø We are to verify (10) assuming that Z k−1 = Ø. On the one hand, (138) can be written as: On the other hand, (139) can be written as (204)-(206), shown at the top of the page, and so Thus adding (203) and (207) we get: Thus the single-step PHD update does result in (10), (210), as shown at the top of the page.