Murat Uney
email: m.uney@ed.ac.uk

Fellow, IEEE Bernard Mulgrew
email: b.mulgrew@ed.ac.uk

Member, IEEE Murat Üney

Member, IEEE Daniel E Clark
email: daniel.clark@telecom-sudparis.eu

Latent Parameter Estimation in Fusion Networks Using Separable Likelihoods

Keywords: sensor networks, hidden Markov models, Markov random fields, pseudo-likelihood, simultaneous localisation and tracking, Monte Carlo algorithms, dynamical Markov random fields

come

Latent parameter estimation in fusion networks using separable likelihoods

I. INTRODUCTION

A wide range of sensing applications including wide area surveillance is underpinned by state space models which are capable of representing a variety of dynamic phenomena such as spatio-temporal (see, e.g., [START_REF] Sarkka | Spatiotemporal learning via infinite-dimensional Bayesian filtering and smoothing: A look at Gaussian process regression through Kalman filtering[END_REF]) processes. In fusion (or, object tracking [START_REF] Vo | Wiley Encyclopedia of Electrical and Electronics Eng., Multitarget Tracking[END_REF]) networks, multi-sensor versions of stochastic state space models, also known as hidden Markov models [START_REF] Cappé | Inference in Hidden Markov Models[END_REF], are used to estimate object trajectories in a surveillance region.

These models, however, are often specified by some latent parameters [START_REF] Cappé | An overview of existing methods and recent advances in sequential Monte Carlo[END_REF] some of which are unknown in practice and need to be estimated based on measurements from the state processes (or, objects). Examples of this problem setting in fusion networks include estimation of noise parameters [START_REF] Singh | Approximate likelihood estimation of static parameters in multi-target models[END_REF], sensor biases [START_REF] Lin | Multisensor multitarget bias estimation for general asynchronous sensors[END_REF], [START_REF] Ristic | Calibration of multi-target tracking algorithms using non-cooperative targets[END_REF] and localisation/calibration of sensors in a GPS denying environment (e.g., in underwater sensing [START_REF] Erol-Kantarci | A survey of architectures and localization techniques for underwater acoustic sensor networks[END_REF]) using point detections from non-cooperative targets [START_REF] Kantas | Distributed maximum likelihood for simultaneous self-localization and tracking in sensor networks[END_REF], [START_REF] Uney | Target aided online sensor localisation in bearing only clusters[END_REF]. Another example is the estimation of the orientations and positions of nodes in a camera network based on feature detections [START_REF] Devarajan | Calibrating distributed camera networks[END_REF].

Such problems fall in the domain of parameter estimation in state space models (see, e.g., [START_REF] Kantas | On particle methods for parameter estimation in state-space models[END_REF] for a review). The parameter likelihood of the multi-sensor problem, however, does not scale well with the number of sensors which specifies the dimensionality of the unknown, or, the length of the measurement window that will be used for estimation. In the presence of multiple objects, the scalability issue is exacerbated by the measurement origin (or, data association) uncertainties that arise. Exact evaluation of the likelihood in this case has combinatorial complexity with the number of sensors [START_REF] Chen | Data association based on optimization in graphical models with application to sensor networks[END_REF], and, in general multiple object models, it is intractable even for a single sensor [START_REF] Jiang | Bayesian tracking and parameter learning for non-linear multiple target tracking models[END_REF]. Estimation using a maximum likelihood (ML) or a Bayesian approach requires repeated evaluation of this likelihood (see, e.g., [START_REF] Kantas | On particle methods for parameter estimation in state-space models[END_REF], [START_REF] Andrieu | Particle markov chain monte carlo methods[END_REF], [START_REF] Ala-Luhtala | An introduction to twisted particle filters and parameter estimation in non-linear statespace models[END_REF]) necessitating the use of efficient approximation strategies.

Intractable or computationally prohibitive likelihoods have motivated a number of lines of work in the statistics literature including likelihood free methods, or, approximate Bayesian computation [START_REF] Marin | Approximate bayesian computational methods[END_REF], and, composite likelihood/pseudolikelihood approaches [START_REF] Varin | An overview of composite likelihood methods[END_REF]. Likelihood free methods can be used for sampling from the parameter posterior in state space models [START_REF] Ehrlich | Gradient free parameter estimation for hidden Markov models with intractable likelihoods[END_REF] including those capable of modelling multiple objects [START_REF] Yildirim | Parameter estimation in hidden Markov models with intractable likelihoods using sequential Monte Carlo[END_REF]. The latter approach is based on developing surrogates to replace the original likelihood, e.g., block based approximations in maximum likelihood [START_REF] Andrieu | On-line parameter estimation in general state-space models[END_REF]. The pseudolikelihood perspective has been useful in networked settings in which constraints on i) the availability of parts of data, and/or, ii) scalability in processing with the number of sources arise. Examples include surrogates built upon local functions for estimation of parametric probability measures (e.g., exponential family distributions) from distributedly stored high dimensional samples [START_REF] Wiesel | Distributed covariance estimation in gaussian graphical models[END_REF]- [START_REF] Mizrahi | Distributed parameter estimation in probabilistic graphical models[END_REF].

It is not straightforward to find such pseudo-likelihoods for parameter estimation in state space models, however, that can resolve these two issues that arise when there are multiple data sources (or, sensors). It is worthwhile to develop and analyse surrogates that provide scalability with the number of sources, and, are suitable to local computations (e.g., local filtering). In [START_REF] Uney | Cooperative sensor localisation in distributed fusion networks by exploiting non-cooperative targets[END_REF], we proposed a pseudo-likelihood which is a product of "dual-term" approximations replacing their intractable exact counterparts. These approximations are separable in that they can be evaluated using single sensor filtering. This feature underpins scalability with the number of sensors. In [START_REF] Uney | A cooperative approach to sensor localisation in distributed fusion networks[END_REF], we have investigated the quality of the dual-term approximation, and, related it to the level of uncertainty in the prediction and estimation of the underlying state process.

In this work, we propose an alternative pseudo-likelihood which is provably a more accurate approximation for parameter estimation in multi-sensor state space models, under typical operating conditions. This approximation is also separable in that it is a scaled product of quadruple terms each of which can be found using single sensor filtering. In order to exploit this quad-term likelihood when there are multiple objects, extra attention should be paid to the handling of the data association uncertainties. We propose to use a hypothesis based parameterisation for the multi-object state space model as detailed in [START_REF] Jiang | Bayesian tracking and parameter learning for non-linear multiple target tracking models[END_REF]& [START_REF] Yildirim | Calibrating the Gaussian multi-target tracking model[END_REF] in order to facilitate the use of the quad-term surrogate in this setting. In the parameterised model, we explicitly point out the combinatorial complexity of exact likelihood evaluation with the number of sensors. Then, we introduce an empirical Bayesian [START_REF] Carlin | Bayes and empirical Bayes methods for data analysis[END_REF] interpretation of local filtering that facilitates the use of separable likelihoods within this model. These modelling aspects detailing the use of separable likelihoods in hypothesis based multi-object models constitute the second contribution of this work.

Separable likelihoods fit well in distributed fusion archictectures in which locally filtered distributions are transmitted in the network, as opposed to sensor measurements [START_REF]Distributed Data Fusion for Network-Centric Operations[END_REF]. Moreover, they facilitate parameter estimation using a message passing computational structure which is desirable in networked problems. Specifically, the proposed likelihood surrogate together with independent parameter priors lead to a pairwise Markov random field (MRF) posterior model. The marginal distributions of this model approximates posterior marginals of the latent parameters to be estimated. We estimate these marginals iteratively using Belief Propagation (BP) [START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF] which consists of successive message passings among neighbouring nodes and updating of local marginals based on these messages. This computational structure lends itself to decentralised estimation, as well as scalable computation at fusion centre.

As an indication of the approximation quality, we consider the Kullback-Leibler divergence (KLD) [START_REF] Cover | Elements of Information Theory[END_REF] of the quad-term likelihood with respect to the actual pairwise likelihood and relate it to the uncertainties in predicting and estimating the underlying state using individual and joint sensor histories. We show that with more accurate local filters the approximation quality improves and the proposed quad-term separable likelihood has an improved error bound compared to the aforementioned dual-term approximation.

We provide a Monte Carlo algorithm for sensor selfcalibration in this framework for linear Gaussian state space (LGSS) models. The algorithm is based on the nonparametric BP approach [START_REF] Sudderth | Nonparametric belief propagation[END_REF] and involves sampling from the updated marginals followed by quad-term likelihood evaluations in the message passing stage. As BP iterations converge to a fixed point, the empirical average of the samples from the marginals constitute (an approximate) minimum mean squared error (MMSE) estimate of the latent parameters. The edge potential are evaluated using the entire measurement history within a selected time period in an offline fashion which is a strategy similar to particle Markov chain Monte Carlo (MCMC) algorithms [START_REF] Andrieu | Particle markov chain monte carlo methods[END_REF]. As such, we differ from [START_REF] Uney | A cooperative approach to sensor localisation in distributed fusion networks[END_REF] in which windowing of measurements are used for enabling online processing.

Preliminary results of the proposed pseudo-likelihood can be found in [START_REF] Uney | Distributed estimation of latent parameters in state space models using separable likelihoods[END_REF]. This article provides a complete account of our solution strategy in multiple object models and is structured as follows: Section II provides the probabilistic model and the problem statement. Then, we detail pairwise pseudolikelihoods in parameterised multi-object models, and, relate this perspective to latent parameter estimation via inference over pairwise MRFs, in Section III. The proposed quad-term node-wise separable likelihood approximation is detailed in Section IV. Section V details the structural and computational properties of the quad-term approximation when the unknowns are respective quantities. Based on these results, we propose a distributed sensor localisation algorithm in linear Gaussian multi-object state space models in Section VI. The efficacy of this algorithm is demonstrated in comparison to the approach in [START_REF] Uney | A cooperative approach to sensor localisation in distributed fusion networks[END_REF], in Section VII. Finally, we conclude in Section VIII.

II. PROBLEM DEFINITION

A. Probabilistic model

Let us consider a set of sensors V " t1, . . . , N u networked over communication links listed by E Ă V ˆV. The graph G " pV, Eq is undirected (i.e., the links are bi-directional), connected, and, might contain cycles.

Next, let us consider a single object with state evolution modelled as a Markov process X k for time index k ě 1. This process is specified by an initial state distribution and a transition density. The state space model with parameters θ is then specified as follows [START_REF] Kantas | On particle methods for parameter estimation in state-space models[END_REF]: The state value x k is a point in the state space X and is generated by the chain

X k |pX 1:k´1 " x 1:k´1 q " πpx k |x k´1 ; θq, X 1 " π b px 1 ; θq, (1)
where .|. denotes conditioning. A measured value z i k P Z i at sensor i P V is generated independently in accordance with the likelihood model

Z i k |pX 1:k " x 1:k , Z i 1:k " z i 1:k q " g i pz i k |x k ; θq (2)
where subscript 1 : k indicates a vector concatenation over time.

In fusion scenarios, there are multiple such objects denoted by a multi-object state

X k fi rX k,1 , . . . , X k,M k s , (3)
that induce measurements according to the above state space model resulting with sensors collecting a multitude of measurements

Z i k fi " Z i k,1 , . . . , Z i k,O i k ı , (4)
where M k is the number of objects and O i k is the number of measurements collected at sensor i at time k. Here, the origin of Z i k,j s are unknown, i.e., the data associations which encode a mapping from these measurement (random) variables to the elements of X k (and, equivalently to the previously collected measurements from the same objects) are not known [START_REF] Vo | Wiley Encyclopedia of Electrical and Electronics Eng., Multitarget Tracking[END_REF].

In the general multi-object tracking model, M k and O i k are random variables with laws determined by probability models regarding how these objects appear in the surveillance region and disappear (which is often referred to as their birth and death, respectively), the law for the false alarms, etc. For the sake of simplicity and ease of presentation in the limited space especially when relating computational complexity to the number of sensors and objects in the following discussion, we assume that all of the objects that exist at time step k " 1 remain in the scene for the time window considered and there are no missed detections and false alarms in sensor measurements which imply that M k " M and O i k " M k , respectively, for some positive integer M with probability one.

In this simplified "closed world" model the multi-object state transition is given by

πpX k |X k´1 ; θq " M ź m"1 πpx k,m |x k´1,m ; θq. ( 5
)
The likelihood of the measurements collected by sensor i is conditioned not only on the multi-object state X k , but, also on a (data association) hypothesis τ i k that encodes the association of measurements to the objects within X k :

l i pZ i k |X k , τ i k ; θq " M ź o"1 g i pz i k,o |x k,τ i k poq ; θq, (6)
and, the prior on τ i k assigns equal probability to all M ! permutations of r1, . . . , M s that τ i k can take, i.e.,

ppτ i k q " 1 M ! . (7)

B. Statement of the problem

We are interested in estimating θ P B using the measurements collected across the network by sensors i P V for a time window of length t. The parameter likelihood of the problem quantifies how well these measurements fit into the state space model with the selected value of the parameter, and, is evaluated via multi-sensor filtering [4, Sec.IV]:

l `Z1 1:t , . . . , Z N 1:t |θ ˘" t ź k"1 p `Z1 k , . . . , Z N k |Z 1 1:k´1 , . . . , Z N 1:k´1 , θ ˘, (8)
where the time updates on the right hand side are given by

p `Z1 k , . . . , Z N k |Z 1 1:k´1 , . . . , Z N 1:k´1 , θ ˘" ÿ τ 1 k ¨¨¨ÿ τ N k ppτ 1 k q ˆ. . . ˆppτ N k q ˆżX M lpZ 1 k , . . . , Z N k |X k , τ 1 k , . . . , τ N k ; θq ˆppX k |Z 1 1:k´1 , . . . , Z N 1:k´1 ; θqdX k , (9)
and, the multi-sensor likelihood inside the integration factorises as

lpZ 1 k , . . . , Z N k |X k , τ 1 k , . . . , τ N k ; θq " ź iPV l i pZ i k |X k , τ i k ; θq. ( 10
)
where the terms in the product are given by [START_REF] Lin | Multisensor multitarget bias estimation for general asynchronous sensors[END_REF].

Here, [START_REF] Erol-Kantarci | A survey of architectures and localization techniques for underwater acoustic sensor networks[END_REF] follows from the chain rule of probabilities. The term in [START_REF] Kantas | Distributed maximum likelihood for simultaneous self-localization and tracking in sensor networks[END_REF] is the contribution at time step k which updates the likelihood of the previous time step and is found using the Markov property that the sensor measurements are mutually independent of the measurement histories, conditioned on the current state for any value of θ. Let us denote this relation by Z j k K KZ j 1:k´1 |X k , θ for i P V (see, e.g., [START_REF] Wainwright | Graphical models, exponential families, and variational inference[END_REF], for this notation). [START_REF] Uney | Target aided online sensor localisation in bearing only clusters[END_REF] follows from that the measurements of different sensors are mutually independent, i.e., Z i k K KZ j k |X k , θ for pi, jq P V ˆV.

This likelihood can be used in a MMSE estimator of θ P B, in principle, for a random variable Θ associated with a prior density ppθq. This estimate is given by the expected value of the posterior distribution ppθ|Z 1 1:t , . . . , Z N 1:t q 9 lpZ 1 1:t , . . . , Z N 1:t |θq ppθq,

θ "

ż B θ ppθ|Z 1 1:t , . . . , Z N 1:t q dθ. ( 12
)
The MMSE estimate can be computed by generating L samples from the posterior distribution in (11) using, for example, MCMC methods [START_REF] Andrieu | Particle markov chain monte carlo methods[END_REF] and using these samples to find a Monte Carlo estimate of the integral in [START_REF] Kantas | On particle methods for parameter estimation in state-space models[END_REF]. In both this approach and maximum likelihood (ML) solutions aiming to maximise [START_REF] Erol-Kantarci | A survey of architectures and localization techniques for underwater acoustic sensor networks[END_REF] with iterative optimisation, repeated evaluations of the likelihood are required. The evaluation of this likelihood is intractable, however, not only because of the pM !q N summations in [START_REF] Kantas | Distributed maximum likelihood for simultaneous self-localization and tracking in sensor networks[END_REF], but, also because of the complexity in finding the integrations involved. The integrands here are i) the multi-sensor likelihood in [START_REF] Uney | Target aided online sensor localisation in bearing only clusters[END_REF], and, ii) the prediction density for X k based on the network's entire measurement history up to time k. In other words, [START_REF] Kantas | Distributed maximum likelihood for simultaneous self-localization and tracking in sensor networks[END_REF] is the scale factor for the posterior density of Bayesian recursions, or, the "centralised" filter given by

ppX k , τ 1:N k |Z 1 1:k , . . . , Z N 1:k ; θq " lpZ 1 k , . . . , Z N k |X k , τ 1:N k ; θq p `Z1 k , . . . , Z N k |Z 1 1:k´1 , . . . , Z N 1:k´1 , θ ˘ˆppτ 1:N k q ˆppX k |Z 1 1:k´1 , . . . , Z j 1:k´1 , θq, (13)
ppX k |Z 1 1:k´1 , . . . , Z N 1:k´1 , θq " ÿ τ 1 k´1 ¨¨¨ÿ τ N k´1 ż X M πpX k |X k´1 , θq ˆppX k´1 , τ 1:N k´1 |Z 1 1:k´1 , . . . , Z N 1:k´1 , θqdX k´1 , (14)
where we denote by τ 1:N k the concatenation of τ i k s and it has pM !q N different configurations. Here, both the prediction [START_REF] Jiang | Bayesian tracking and parameter learning for non-linear multiple target tracking models[END_REF] and update [START_REF] Chen | Data association based on optimization in graphical models with application to sensor networks[END_REF] are OppM !q N q.

In order to address these challenges, multi-object filtering (or, tracking) algorithms often employ two approximations: First, they aim to find the most probable data association hypothesis in [START_REF] Chen | Data association based on optimization in graphical models with application to sensor networks[END_REF] denoted by τ 1:N k , instead of both evaluating this expression for all possible associations and storing them. The benefits of doing so are that i) one can generate tracks (or, object trajectories) as simply marginals of ppX k , τ 1:N k " τ 1:N k | . q for k " 1, ..., t, and, ii) evaluations of the integral in [START_REF] Jiang | Bayesian tracking and parameter learning for non-linear multiple target tracking models[END_REF] in the next time step can be restricted to this value of the hypothesis variable. Equivalently, the posterior distribution in [START_REF] Chen | Data association based on optimization in graphical models with application to sensor networks[END_REF] is factorised as [START_REF] Andrieu | Particle markov chain monte carlo methods[END_REF] where the association variables appear as model parameters in the first term and the second term is similar to a prior distribution on these models with the difference that it is conditioned on the current measurements. At time k ´1, let us select this "empirical prior" as

ppX k , τ 1:N k |Z 1 1:k , . . . , Z N 1:k ; θq " ppX k |Z 1 1:k , . . . , Z N 1:k , θ, τ 1:N k qppτ 1:N k |Z 1 1:k , . . . , Z N 1:k , θq " ppX k |Z 1 1:k , . . . , Z N 1:k , θ, τ 1:N k qppτ 1:N k |Z 1 k , . . . , Z N k , θq
ppτ 1:N k´1 |Z 1 k´1 , . . . , Z N k´1 , θq Ð δ τ 1:N k´1 pτ 1:N k´1 q ( 16
)
where δ is Kronecker's delta function and Ð denotes assignment. The second approximation follows from the first one: Evaluation of the prediction stage in [START_REF] Jiang | Bayesian tracking and parameter learning for non-linear multiple target tracking models[END_REF] reduces to evaluation of the Chapman-Kolmogorov equation for only the most likely value of the association parameters. This approach is often referred to as empirical Bayes [START_REF] Carlin | Bayes and empirical Bayes methods for data analysis[END_REF], and is used to facilitate approximate solutions to otherwise intractable problems.

A similar approximation can be used when evaluating the parameter posterior in [START_REF] Devarajan | Calibrating distributed camera networks[END_REF]. This leads to the following likelihood

l `Z1 1:t , . . . , Z N 1:t |θ, τ 1:N 1:t ˘" t ź k"1 p τ 1:N k `Z1 k , . . . , Z N k |Z 1 1:k´1 , . . . , Z N 1:k´1 , θ ˘, (17)
conditioned on τ1:N 1:t where the factors are defined by

p τ 1:N k `Z1 k , . . . , Z N k |Z 1 1:k´1 , . . . , Z N 1:k´1 , θ ˘fi ż X M lpZ 1 k , . . . , Z N k |X k , τ 1:N k ; θq ˆppX k |Z 1 1:k´1 , . . . , Z N 1:k´1 ; θqdX k . ( 18
)
This likelihood evaluated at τ 1:N 1:k " τ 1:N 1:k replaces the one in [START_REF] Erol-Kantarci | A survey of architectures and localization techniques for underwater acoustic sensor networks[END_REF] when the empirical (model) prior is selected as in [START_REF] Ala-Luhtala | An introduction to twisted particle filters and parameter estimation in non-linear statespace models[END_REF] (see Appendix A for details). We will refer to [START_REF] Marin | Approximate bayesian computational methods[END_REF] as the empirical likelihood. Note that [START_REF] Varin | An overview of composite likelihood methods[END_REF] is the integral term in [START_REF] Kantas | Distributed maximum likelihood for simultaneous self-localization and tracking in sensor networks[END_REF].

The empirical likelihood update term is computationally more convenient, however, alone it is not sufficient for scalability with the number of sensors N : Finding τ 1:N k is equivalently an N `1-dimensional assignment problem which is NP hard even for N " 2 sensors [START_REF] Poore | Some assignment problems arising from multiple target tracking[END_REF] which partly underlies the local filtering paradigm for multi-sensor processing and our interest in compatible solutions. For a moment, let us consider the problem for a single object, i.e., for M " 1. In this case τ i k for i " 1, . . . , N have only one possible configuration (i.e., there is no data association uncertainty). Because the dimensionality of θ is specified by N and (8) will be evaluated for roughly N L samples (when estimating [START_REF] Kantas | On particle methods for parameter estimation in state-space models[END_REF] (see, e.g., [START_REF] Ristic | Beyond the Kalman Filter; particle filters for tracking applications[END_REF])) each of which costing -in the simplest linear Gaussian measurements case (9) 1 -at the least OpN 2 tq, the computational cost will be cubic in the number of sensors which can easily become prohibitive for large N .

The networked setting has additional constraints to take into account: The sensors perform local filtering of their measurements and exchange filtered (track) distributions over G as opposed to transmitting their measurements [START_REF]Distributed Data Fusion for Network-Centric Operations[END_REF]. As a result, the network-wide measurements are not available to evaluate the likelihood of the problem. Instead, local distributions we denote by ppX k , τ j k " τ j k |Z j 1:k q are made available to neighbouring nodes where τ j k is an approximation to the most probable association configuration τ j k found locally, based on only the local sensor measurements at sensor j. There are computationally efficient algorithms for finding such solutions for the single sensor problem (see, e.g., [START_REF] Poore | Some assignment problems arising from multiple target tracking[END_REF] and the references therein). Therefore, a viable solution needs to build upon these densities and local data associations τ i k as opposed to joint multi-sensor filtering in the network.

The problem we address in this work is the design of scalable approximations to [START_REF] Erol-Kantarci | A survey of architectures and localization techniques for underwater acoustic sensor networks[END_REF] for estimating θ in a networked setting based on local filtering results at the nodes. The proposed approach also addresses the aforementioned computational bottleneck at fusion centres in centralised multisensor architectures with a designated node receiving unfiltered sensor measurements.

It is also worth noting that the parameter vector θ P B can be used to represent a wide variety of parameters of the global model some of which can be intrinsic to sensors i P V individually such as parameters pertaining to local noise models. We are particularly interested in a second class of parameters which have dependencies among sensors such as respective parameters such as sensor locations and similar "calibration" parameters. In the former setting, the estimation of local parameters decouple into independent estimation problems which can be solved using a suitable approach (see, e.g., [START_REF] Yildirim | Calibrating the Gaussian multi-target tracking model[END_REF], [START_REF] Kokkala | Combining particle MCMC with Rao-Blackwellized Monte Carlo data association for parameter estimation in multiple target tracking[END_REF], [START_REF] Schlangen | Single-cluster PHD filter methods for joint multi-object filtering and parameter estimation[END_REF]).

In our setting, θ fi rθ 1 , . . . , θ N s where θ i is associated with i P V and its estimation does not decouple and depends on all measurements across the network due to the dependencies of parameters, which, in turn, brings forward the multi-sensor aspects of the problem this work aims to address. Because local filtering is performed, on the other hand, local estimation of data association τ i k is available independent of θ, which we discuss in detail later in Section V.

III. A PSEUDO-LIKELIHOOD AND A PAIRWISE MRF POSTERIOR FOR DECENTRALISED ESTIMATION

Pseudo-likelihoods are constructed from likelihood like functions which are computationally convenient and defined typically over smaller subsets of the data to overcome difficulties posed by intractable likelihoods over the entire set of data (see, e.g., [START_REF] Varin | An overview of composite likelihood methods[END_REF] and the references therein). Let us denote the network-wide data set by Z fi rZ 1 1:t , . . . , Z N 1:t s. A fairly general form for a pseudo-likelihood is given by [START_REF] Varin | An overview of composite likelihood methods[END_REF] lpZ|θq " where S is an index set, ω s is a positive real number, and, Z ds and Z cs are mutually exlusive subsets of Z (for example, Z d1 " Z i 2 and Z c1 " Z j 1 , etc.). These sets can be selected in various ways ensuring that the factors are computationally convenient functions, for example, marginals and/or conditional densities, and, estimates based on lpZ|θq are sensible. Note that ( 8) is also in this form, however, with difficult to evaluate factors.

ź sPS lpZ ds |Z cs , θq ωs (19)
Let us consider θ " rθ 1 , ..., θ N s and a pseudo-likelihood surrogate for the empirical likelihood in ( 17):

lτ 1:N 1:t pZ|θq " ź pi,jqPE l τ i,j 1:t pZ i , Z j |θ i,j q, ( 20
) " ź pi,jqPE t ź k"1 p τ i,j k pZ i k , Z j k |Z i 1:k´1 , Z j k´1 , θ i,j q.( 21
)
Here, E is essentially the set of sensor pairs whose likelihoods we would like to incorporate into the pseudo-likelihood, and, it is convenient to choose them as those that share a communication link, in a networked setting (Section II-A). The pairwise structure above is beneficial to use with the MMSE estimator in [START_REF] Kantas | On particle methods for parameter estimation in state-space models[END_REF]. Note that the MMSE estimate is the concatenation of the expected values of posterior marginals, i.e., ppθ i |Zq for i " 1, . . . , N . These distributions can be found using message passing algorithms over G when the surrogate [START_REF] Yildirim | Parameter estimation in hidden Markov models with intractable likelihoods using sequential Monte Carlo[END_REF] is used in [START_REF] Devarajan | Calibrating distributed camera networks[END_REF] together with independent but arbitrary a priori distributions selected for Θ i s. Specifically, the parameter posterior corresponding to such a selection of prioirs and the pseudo-likelihood ( 20) is a pairwise Markov random field over G " pV, Eq [START_REF] Wainwright | Graphical models, exponential families, and variational inference[END_REF]:

ppθ|Zq 9 ź iPV ψ i pθ i q ź pi,jqPE ψ ij pθ i , θ j q, ( 22
)
ψ i pθ i q " p 0,i pθ i q, ψ ij pθ i , θ j q " l τ i,j 1:t pZ i , Z j |θ i , θ j q, ( 23
)
where the node potential functions (i.e., ψ i s) are the selected priors (e.g., uniform distributions over bounded sets θ i s take values from) and the edge potentials (i.e., ψ ij s) are the pairwise likelihoods for the pairs pi, jqs. This model is illustrated in Fig. 1.

The pairwise MRF model in [START_REF] Wiesel | Distributed covariance estimation in gaussian graphical models[END_REF] allows the computation of posterior marginal ppθ i |Zq through iterative local message passings such as Belief Propagation (BP) [START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF]. In BP, the nodes maintain distributions over their local variables and update them based on messages from their neighbours which summarise the information neighbours have gained on these variables. This is described for all i P V by m ji pθ i q " ż ψ ij pθ i , θ j q ψ j pθ j q ź i 1 Pnepjqzi m i 1 j pθ j q dθ j , (24)

pi pθ i q 9 ψ i pθ i q ź jPnepiq m ji pθ i q.

(25)

In BP iterations, nodes simultaneously send messages to their neighbours using (24) (often using constants as the previously received messages during the first step) and update their local "belief" using [START_REF] Uney | Cooperative sensor localisation in distributed fusion networks by exploiting non-cooperative targets[END_REF]. If G contains no cycles (i.e., G is a tree), pi s are guaranteed to converge to the marginals of [START_REF] Wiesel | Distributed covariance estimation in gaussian graphical models[END_REF], in a finite number of steps [START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF]. For the case in which G contains cycles, iterations of ( 24) and ( 25) are known as loopy BP (LBP). For the case, convergence does not have general guarantees, nevertheless LBP has been been very successful in computing approximate marginals in a distributed fashion, in fusion, self-localisation and tracking problems in sensor networks [START_REF] Ihler | Nonparametric belief propagation for self-localization of sensor networks[END_REF]- [START_REF] Uney | Graphical model-based approaches to target tracking in sensor networks: An overview of some recent work and challenges[END_REF]. In our problem setting, we assume that the models over spanning trees of a loopy G are consistent in that they lead to "similar" marginal parameter distributions, which suggests the existence of LBP fixed points [START_REF] Wainwright | Tree-based reparameterization framework for analysis of sum-product and related algorithms[END_REF] that will be converged when initial beliefs are selected reasonably [START_REF] Yedidia | Constructing free-energy approximations and generalized belief propagation algorithms[END_REF].

IV. QUAD-TERM NODE-WISE SEPARABLE LIKELIHOODS

The pseudo-likelihood introduced in Section III leads to a parameter posterior that admits a pairwise MRF model. This is advantageous in providing a means for decentralised estimation through message passing algorithms in a network. The edge potentials [START_REF] Liu | Distributed parameter estimation via pseudolikelihood[END_REF] of this model, however, are i) conditioned jointly on two sensors' measurements simultaneous access to which is infeasible in a networked setting, and, ii) conditioned on association variables for two sensors and ideally should be evaluated at its most probable configuration τ i,j 1 , . . . , τ i,j t each of which is NP-hard to find, as explained in Section II.

In order to overcome these difficulties, we introduce an approximation which factorises into terms local to nodes, i.e., a node-wise separable approximation. Let us consider the "centralised" pairwise likelihood update term in [START_REF] Andrieu | On-line parameter estimation in general state-space models[END_REF] given some configuration τ i,j k for k " 1, . . . , t, and, drop them from the subscript for the sake of simplicity in notation, as well as the i, j subscript in θ, in the following discussion. This term factorises in alternative ways as follows:

ppZ i k , Z j k |Z i 1:k´1 , Z j 1:k´1 , θq " ppZ i k |Z i 1:k´1 , Z j 1:k , θqppZ j k |Z i 1:k´1 , Z j 1:k´1 , θq (26)
" ppZ j k |Z i 1:k , Z j 1:k´1 , θqppZ i k |Z i 1:k´1 , Z j 1:k´1 , θq (27)
" ´ppZ i k |Z i 1:k´1 , Z j 1:k , θqppZ j k |Z i 1:k´1 , Z j 1:k´1 , θq ¯1{2 ˆ´ppZ j k |Z i 1:k , Z j 1:k´1 , θqppZ i k |Z i 1:k´1 , Z j 1:k´1 , θq ¯1{2 (28)
In the first and second lines above, the chain rule is used. The third equality can be found by taking the geometric mean of the first two expressions. All four factors in Eq.( 28) are conditioned on the measurement histories of both sensors to which one cannot have simulatenous access in a networked setting. We would like to aviod this by leaving out the history of sensor i (sensor j) in the first two (last two) terms of (28), i.e.,

qpZ i k , Z j k |Z i 1:k´1 , Z j 1:k´1 , θq fi 1 κ k pθq ´ppZ i k |Z j 1:k , θqppZ j k |Z j 1:k´1 , θq ¯1{2 ˆ´ppZ j k |Z i 1:k , θqppZ i k |Z i 1:k´1 , θq ¯1{2 (29)
κ k pθq " ż ż ´ppZ i k 1 , Z j k 1 |Z j 1:k´1 , θq ˆppZ i k 1 , Z j k 1 |Z i 1:k´1 , θq ¯1{2 dZ i k 1 dZ j k 1 ( 30
)
where κ k pθq is the normalisation constant that guarantees q to integrate to unity. Note that κ k is a function of the parameters θ.

The appeal of this quadruple term is that its factors depend on single sensor histories. As such, they require filtering of sensor histories of i and j individually enabling the evaluation of their product in a network. This point is discussed later in this section.

A. Approximation quality

We consider the difference between the original centralised update term in [START_REF] Carlin | Bayes and empirical Bayes methods for data analysis[END_REF] and the quad-term approximation introduced in [START_REF]Distributed Data Fusion for Network-Centric Operations[END_REF]. Because these terms are probability densities over sensor measurements, their "divergence" can be quantified using the KLD [START_REF] Cover | Elements of Information Theory[END_REF]:

Proposition 4.1: The KLD between the centralised update and the node-wise separable approximation in [START_REF]Distributed Data Fusion for Network-Centric Operations[END_REF] is bounded by the average of the mutual information (MI) [START_REF] Cover | Elements of Information Theory[END_REF] between the current measurement pair and a single sensor's history conditioned on the history of the other sensor, i.e.,

DpppZ i k , Z j k |Z i 1:k´1 , Z j 1:k´1 , θq||qpZ i k , Z j k |Z i 1:k´1 , Z j 1:k´1 , θqq ď 1 2 IpZ i k , Z j k ; Z i 1:k´1 |Z j 1:k´1 , Θq `1 2 IpZ i k , Z j k ; Z j 1:k´1 |Z i 1:k´1 , Θq. (31)
The proof can be found in Appendix B 2 . The upper bound in [START_REF] Cover | Elements of Information Theory[END_REF] measures the departure of the current pair of measurements, and, one of the sensor histories from a state of conditional independence when they are conditioned on the history of the other sensor. Note that these variables, when conditioned on X k , are conditionally independent, i.e., pZ i k , Z j k qK KZ j 1:k´1 |X k , Θ holds and consequently

IpZ i k , Z j k ; Z i 1:k´1 |X k , Θq " IpZ i k , Z j k ; Z j 1:k´1 |X k , Θq " 0.
Similarly, the average MI term on the right hand side of [START_REF] Cover | Elements of Information Theory[END_REF]

is zero if pZ i k , Z j k qK KZ i 1:k´1 |Z j 1:k´1 , Θ and pZ i k , Z j k qK KZ j 1:k´1 |Z i 1:
k´1 , Θ hold simultaneously. This condition is satisfied, for example, in the case that either of the measurement histories Z i 1:k´1 and Z j 1:k´1 are sufficient statistics for X k (i.e., it can be predicted by both sensors with probability one). This level of accuracy should not be expected as the transition density of state space models introduce some uncertainty. Therefore, it is instructive to relate the KLD in [START_REF] Cover | Elements of Information Theory[END_REF] further to the uncertainty on X k given the sensor histories: Corollary 4.2: The KLD considered in Proposition 4.1 is upper bounded by the weighted sum of uncertainty reductions in the local target prediction and posterior distributions achieved when the other sensor's history is included jointly:

DpppZ i k , Z j k |Z i 1:k´1 , Z j 1:k´1 , θq||qpZ i k , Z j k |Z i 1:k´1 , Z j 1:k´1 , θqq ď 1 2 ˜´HpX k |Z j 1:k´1 , Θq ´HpX k |Z j 1:k´1 , Z i 1:k´1 , Θq HpX k |Z i k´1 , Θq ´HpX k |Z j 1:k´1 , Z i 1:k´1 , Θq ¯1
2 ˜´HpX k |Z j 1:k , Θq ´HpX k |Z j 1:k , Z i 1:k´1 , Θq HpX k |Z i 1:k , Θq ´HpX k |Z i 1:k , Z j 1:k´1 , Θq ¯¸, ( 32
)
where H denotes the Shannon differential entropy [START_REF] Cover | Elements of Information Theory[END_REF].

The proof is provided in Appendix C. Corollary 4.2 relates the approximation quality of the quad-term node-wise separable updates to the uncertainties in the target state prediction and posterior distributions when individual node histories and their combinations are considered. The difference terms on the RHS of (32) quantify the difference in uncertainty between estimating the target state X k using only the local measurements, and, also taking into account the other sensor's measurements. Overall, a better quality of approximation should be expected when the local filtering densities involved concentrate around a single point in the state space.

B. The quad-term pairwise likelihood

The quad-term update in (29) leads to a separable approximate likelihood given by

l ´Zi 1:t , Z j 1:t |θ ¯" t ź k"1 qpZ i k , Z j k |Z i 1:k´1 , Z j 1:k´1 , θq (33)
We refer to this term as the quad-term separable likelihood as it can also be expressed as a (scaled) product of four factors each of which are the products of the four factors of (29) over k. Let us define

r k ij pZ i k , θq fi ppZ i k |Z j 1:k , θq, s k j pZ j k , θq fi ppZ j k |Z j 1:k´1 , θq.
Then, the quad-term update in [START_REF] Uney | Distributed estimation of latent parameters in state space models using separable likelihoods[END_REF] is given by

qpZ i k , Z j k |Z i 1:k´1 , Z j 1:k´1 , θq " 1 κ k pθq ´rk ij pZ i k , θqs k j pZ j k , θq ¯1{2 ´rk ji pZ j k , θqs k i pZ i k , θq ¯1{2 , (34)
where the normalisation factor is given in [START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF], and, equivalently in terms of the four factors above as

κ k pθq " ż ż ´rk ij pZ i k 1 , θqs k j pZ j k 1 , θq ¯1{2 ˆ´r k ji pZ j k 1 , θqs k i pZ i k 1 , θq ¯1{2 dZ i k 1 dZ j k 1 .

Corollary 4.3:

The KLD between the parameter likelihood in [START_REF] Erol-Kantarci | A survey of architectures and localization techniques for underwater acoustic sensor networks[END_REF] and the node-wise separable approximation in [START_REF] Uney | Distributed estimation of latent parameters in state space models using separable likelihoods[END_REF] is bounded by the terms on the right hand sides of ( 31) and [START_REF] Sudderth | Nonparametric belief propagation[END_REF] summed over k " 1, . . . , t as

D ´l ´Zi 1:t , Z j 1:t |θ ¯|| l ´Zi 1:t , Z j 1:t |θ ¯¯" t ÿ k"1 D pp||qq . ( 35
)
Proof. Eq. ( 35) can easily be found after expanding the KLD term explicitly and expressing the logarithm of products involved as sums over logarithms of the factors. Boundedness follows from non-negativity of KLDs and summing both sides of ( 31) and ( 32) over k " 1, ..., t.

As a conclusion, when t is not large -e.g., on the order of tens which is typical in fusion applications-the proposed approximation can be used for parameter estimation via local filtering. Sensors that are more accurate in inferring the underlying state process result with a smaller KLD in [START_REF] Poore | Some assignment problems arising from multiple target tracking[END_REF], which in turn leads to a more favourable estimation performance.

One other approximation based on local filtering distributions was studied in [START_REF] Uney | A cooperative approach to sensor localisation in distributed fusion networks[END_REF] which has the following dual-term product form

upZ i k , Z j k |Z i 1:k´1 , Z j 1:k´1 , θq fi ppZ i k |Z j 1:k´1 , θqppZ j k |Z i 1:k´1 , θq. (36)
In Appendix D, we shown that Dpp||qq ă Dpp||uq, when sensors are equivalent. The entropy bound given in [START_REF] Sudderth | Nonparametric belief propagation[END_REF] is also smaller than that for the dual-term approximation. In other words, the quad-term approximation is more accurate compared to the dual-term approximation, under typical operating conditions.

The scaling factor of the dual term approximation is unity regardless of θ, on the other hand, admitting a significant amount of flexibility in the range of the distributions and likelihoods that can be accommodated in the state space model. For example, the dual-term pseudo-likelihood is used with random finite set variables (RFS) in [START_REF] Uney | A cooperative approach to sensor localisation in distributed fusion networks[END_REF], [START_REF] Üney | Distributed localisation of sensors with partially overlapping field-of-views in fusion networks[END_REF], which, in a sense, have the association variables marginalised out making it possible to avoid multi-dimensional assignment problems in the general multi-object tracking model. For RFS distributions, however, it is not straightforward to compute the scaling factor in [START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF] for the quad-term. In this article, we consider a parametric model instead, which is effectively configured through association variables.

V. QUAD-TERM LIKELIHOOD FOR SENSOR CALIBRATION

PARAMETERS

The results presented so far are fairly general and do not depend on the nature of Θ. When Θ represents respective parameters such as calibration parameters, there are certain simplifications of the expressions involved which provide computational benefits. In particular, parameters such as respective location and bearing angles relate the local coordinate frames of the sensors which collect measurements in their local frame. The local filtering distributions are hence over the space of state vectors in the local frame. A point x k P X (Section II-A) is implicitly in the Earth coordinate frame (ECF), and, associated with its representation in the jth local frame rx k s j through a coordinate transform T with the following properties rx k s j " T px k ; θ j q, (37) rx k s i " T pT ´1prx k s j ; θ j q; θ i q.

As an example, when x k is a location on the Cartesian plane, and, θ j is the position of sensor j, T is given by T px k ; θ j q fi x k ´θj , T pT ´1prx k s j ; θ j q; θ i q " rx k s j `θj ´θi .

For simplicity in notation, we will denote T pT ´1p.; θ j q; θ i q by T θ p.q when semantics is clear from the context 3 .

In this section, it is revealed how local filtering distributions are used in the quad-term update. When θ are respective quantities, these distributions become independent of θ because both the state and the measurement variables are in the same local coordinate frame. In other words, at sensor j

ppX k , τ j k |Z j 1:k , θ j q " pprX k s j , τ j k |Z j 1:k q, (38)
holds for the filtering posterior, for any configuration of τ j k . In the prediction stage of filtering, the multi-object transition kernel in (5) also becomes independent of Θ, so, the Chapman-Kolmogorov equation for finding the prediction density at sensor j (together with the empirical Bayes selection of association priors during iterations as explained in Section II) becomes pprX k s j |Z j 1:k´1 q " ż πprX k s j |rX k´1 s j q ˆpprX k´1 s j , τ j k´1 " τ j k´1 |Z j 1:k´1 qdrX k s j . (39) Note also that the entropy terms in [START_REF] Sudderth | Nonparametric belief propagation[END_REF] that are conditioned on a single sensor's measurements measure the uncertainty of the above densities. Consequently, they also become independent of Θ, i.e., HpX k |Z j 1:k´1 , Θq equals to HpX k |Z j 1:k´1 q for example, which highlights its relevance to the local prediction accuracy.

A. The quad-term update for calibration

Now, let us expand the quad-term time updates in [START_REF] Wainwright | Graphical models, exponential families, and variational inference[END_REF], and, explicitly show the aforementioned simplifications. We start with s k j which is the scale factor of the local Bayesian filter at sensor j:

s k j pZ j k , θq " p τ j k pZ j k |Z j 1:k´1 , θq, " ż l j pZ j k |X k , τ j k , θqppX k |Z j 1:k´1 , θqdX k , " ż l j pZ j k |rX k s j , τ j k qpprX k s j |Z j 1:k´1 qdrX k s j , (40)
where in the last line independence from θ is asserted. Because s k j does not depend on θ, we denote s k j pZ j k , θq by s k j pZ j k q in the rest of the article.

Next, let us consider r k ij which has terms in different coordinate frames:

r k ij pZ i k , θq " p τ i,j k pZ i k |Z j 1:k , θq " ż l i pZ i k |X k , τ i k , θqppX k , τ j k |Z j 1:k , θqdX k " ż l i `Zi k |T pX k ; θ i q, τ i k ˘p ´T ´1pX k ; θ j q, τ j k |Z j 1:k ¯dX k " ż l i `Zi k |T θ prX k s j q, τ i k ˘p ´rX k s j , τ j k |Z j 1:k ¯drX k s j . (41)
In the third line above, the coordinate transformations are substituted explicitly. The last line follows from that the filtering distribution [START_REF] Schlangen | Single-cluster PHD filter methods for joint multi-object filtering and parameter estimation[END_REF] is in the jth local frame.

B. Evaluation of the quad-term update based on single sensor filtering distributions

Here, we discuss the evaluation of the quad-term likelihood given local filtering distributions and single sensor association configurations which we denote for sensor j by pprX k s j , τ j k " τ j k |Z j 1:k q and τ j k , respectively, as explained in Section II-B.

Instead of considering evaluation for the most probable association hypothesis τ i,j k which is infeasible to find, we propose to use the local results τ i,j k fi pτ i k , τ j k q as a reasonable approximation to this configuration and substitute them in ( 40)- [START_REF] Uney | Graphical model-based approaches to target tracking in sensor networks: An overview of some recent work and challenges[END_REF]. These approximations can be found regardless of θ as discussed earlier in this section by using one of the well studied algorithms in the literature [START_REF] Vo | Wiley Encyclopedia of Electrical and Electronics Eng., Multitarget Tracking[END_REF] such as solving a 2 ´D association problem at each time step [START_REF] Poore | Some assignment problems arising from multiple target tracking[END_REF] to find τ j k . We detail this approach for a linear Gaussian state space model in Section VI.

Given these local results and their exchange over the network, one can consider an in-network computation scheme for evaluating (40)- [START_REF] Uney | Graphical model-based approaches to target tracking in sensor networks: An overview of some recent work and challenges[END_REF]. Specifically, these terms (and, the other factors of the quad-term update which are obtained by replacing i and j in these expressions) can be found at the sensor platform where the measurements to be substituted for evaluation are stored, i.e., sensors j and i, respectively, for s k j and r k ij . Substitution of the measurement histories on the conditioning side will have been carried out by local filtering. More explicitly, s k j in (40) (or, s k i ) becomes a product of similar terms when τ k j is substituted in [START_REF] Cetin | Distributed fusion in sensor networks: A graphical models perspective[END_REF], and, its computation is carried out during the local filtering of sensor j's (or, sensor i's) measurements using

s k j pZ j k q " M ź o"1 s k j,o pz j k,o q (42)
s k j,o pz j k,o q fi ż g j pz j k,o |x 1 k qp m 1 px 1 k |Z j 1:k´1 qdx 1 k
where m 1 " τ j k poq and the density in the integral is the m 1 th marginal of the local prediction density, the mth of which is given by

p m px 1 |Z j 1:k´1 q fi ż ppX k " " x k,1 , . . . , x k,m´1 , x 1 , x k,m`1 , x k,m`2 , . . . , x k,M s ,τ j k " τ j k |Z j 1:k´1 q dx k,1 . . . dx k,m´1 dx k,m`1 . . . dx k,M .
The term r k ij in (41) (or, r k ji ) is also computed based on these local filtering distributions. The integration in the RHS of [START_REF] Uney | Graphical model-based approaches to target tracking in sensor networks: An overview of some recent work and challenges[END_REF] implicitly assumes that the ordering of individual objects in the local multi-object vectors are the same. In a networked setting, however, this is not necessarily the case and the identities of the fields in the state vector may differ [START_REF] Guibas | The identity management problem-a short survey[END_REF]. In order to tackle with this unknown correspondance, we introduce an additional permutation random variable γ k (see, e.g., [START_REF] Huang | Fourier theoretic probabilistic inference over permutations[END_REF]) for relating the fields of a multi-object vector X k as ordered locally at sensor i and j, such that the mth field of the state vector at sensor i refers to the same object in the γ k pmqth field of the state vector at sensor j. For example, rx k,m s i " T θ prx k,γ k pmq s j q.

Suppose that an estimate γk of this quantity is provided. After substituting in [START_REF] Uney | Graphical model-based approaches to target tracking in sensor networks: An overview of some recent work and challenges[END_REF] together with τ i,j k one obtains

r k ij pZ i k , θq " M ź o"1 r k ij,o pz i k,o , θq, (43)
r k ij,o pz i k,o , θq fi ż g i `zi k,o |T θ px 1 k q ˘pm 1 px 1 k |Z j 1:k qdx 1 k
where m 1 " γk pτ i k poqq and the density inside the integral is the m 1 th marginal of the filtering distribution local to sensor j. In Appendix E, we show that (43) replaces the likelihood for γ k when m 1 " γ k pτ i k poqq and an ML estimate γk can be found in a way similar to solving the data association problem in local filtering, which is detailed later in Section VI-A.

Finally, the scale factor [START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF] is computed. This involves finding the measurement distributions in (30) using the prediction distribution in [START_REF] Ihler | Nonparametric belief propagation for self-localization of sensor networks[END_REF] for both sensors i and j together with their likelihoods. This leads to the following two decomposition: The first term in [START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF] is found as

ppZ i k 1 , Z j k 1 |Z i 1:k´1 , θq " M ź o"1 p o pz i 1 k,o , z j 1 k,ξpoq |Z i 1:k´1 , θq p o pz i , z j |Z i k´1 , θq fi ż g i pz i |x 1 k qp m 1 px 1 k |Z i 1:k´1 qdx 1 k ˆż g j pz j |T θ ´1px 1 k qqp m 1 px 1 k |Z i 1:k qdx 1 k ( 44
)
where ξpoq " τ j´1 k ˝τ i k poq maps the oth measurement at sensor i to the corresponding one in sensor j, and, m 1 " τ i k poq in the second line.

The second term in [START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF] is found as

ppZ i k 1 , Z j k 1 |Z j 1:k´1 , θq " M ź o"1 p o pz i 1 k,o , z j 1 k,ξpoq |Z j k´1 , θq p o pz i , z j |Z j k´1 , θq fi ż g j pz j |x 1 k qp m 1 px 1 k |Z j 1:k´1 qdx 1 k ˆż g i pz i |T θ px 1 k qqp m 1 px 1 k |Z j 1:k q dx 1 k ( 45
)
where m 1 " γk ˝τ i k poq in the last line is the object that corresponds to the oth measurement at sensor i.

Consequently, the scale factor is found as

κ k pθq " M ź o"1 κ k,o pθq (46)
κ k,o pθq fi ż `po pz i , z j |Z i 1:k´1 , θq ˆpo pz i , z j |Z j 1:k´1 , θq ¯1{2 dz i dz j , using the densities found in ( 44) and ( 45).

The expressions above describe the evaluation of the quadterm likelihood in terms of single sensor filtering distributions that can be obtained using any filtering algorithm with individual measurement histories. The use of the local filtering distributions provides scalability with the number of sensors for parameter estimation in the state space model in Section II.

These computations can be distributed in the pair pi, jq as follows: Both sensors i and j perform local filtering and exchange the resulting posterior densities at every step, as well as s i and s j , respectively, found using [START_REF] Wainwright | Tree-based reparameterization framework for analysis of sum-product and related algorithms[END_REF]. Based on the received densities, sensors i and j evaluate r ij and r ji , respectively, using [START_REF] Yedidia | Constructing free-energy approximations and generalized belief propagation algorithms[END_REF]. As part of the filtering process, they realise the Chapman-Kolmogorov equation in [START_REF] Ihler | Nonparametric belief propagation for self-localization of sensor networks[END_REF] with their local posterior, as well as the remote posterior recently received. These densities are then used in ( 44), [START_REF] Guibas | The identity management problem-a short survey[END_REF], and, [START_REF] Huang | Fourier theoretic probabilistic inference over permutations[END_REF] to compute the scale factor for the next step. The scale factor hence found in the previous step, therefore, is substituted in [START_REF] Wainwright | Graphical models, exponential families, and variational inference[END_REF] together with the four terms computed. This is repeated for k " 1, . . . , t, and, the quad-term likelihood [START_REF] Uney | Distributed estimation of latent parameters in state space models using separable likelihoods[END_REF] is computed, as a result.

VI. A MONTE CARLO LBP ALGORITHM FOR SENSOR CALIBRATION IN LINEAR GAUSSIAN STATE SPACE MODELS

In this section, we consider a linear Gaussian state space (LGSS) model within the probabilistic graphical model in Fig. 1, and, specify an algorithm for estimation of θs that combines the quad-term calibration likelihood evaluation detailed in Section V with BP message passing on the resulting pairwise MRF model (Section III). This algorithm uses Monte Carlo methods for realising BP [START_REF] Sudderth | Nonparametric belief propagation[END_REF] and facilitates scalability by building upon single sensor filtering as required by the quad-term approximation. As a result, an efficient inference scheme over the model in Fig. 1 is achieved.

The state space model we consider is specified by a linear state transition with process noise that is additive and Gaussian, and, linear measurements with independent Gaussian measurement noise, i.e.,

πpx k |x k´1 q " N px k ; Fx k´1 , Qq (47)
g j pz j k |x k ; θ j q " N pz j k ; H j rx k s j , R j q (48)

for j " 1, . . . , N , where N p.; µ, Pq is a multi-dimensional Gaussian density with mean vector µ and covariance matrix P.

Here, x k is the concatenation of position and velocity (on a 2 ´D Euclidean plane, without loss of generality). The matrices F and Q model motion with unknown acceleration (equivalently, manouevres), and, are selected as

F " « I, ∆T ˆI 0, I ff , Q " σ 2
« q 1 I, q 2 I q 2 I, q 3 I ff where I and 0 are the 2 ˆ2 identity and zero matrices, respectively. ∆T is the time difference between consecutive steps. Q is positive definite and parameterised with σ 2 , and, 0 ă q 1 ă q 2 ă q 3 ă 1 specifying the magnitude of the uncertainty, and, contributions of higher order terms, respectively 4 . In the measurement model, R j is the measurement noise covariance, and, H j is the observation matrix which we assume forms an observable pair with F (e.g., H j " rI, 0s).

A. Local single sensor filtering

We now focus on filtering and provide explicit formulae that adopts the recursions in ( 13)-( 18) for a single sensor in the LGSS model. We use the empirical Bayes approach explained in Section II-B for scaling with time under data association uncertainties. This approach corresponds to the single frame data association solution in multi-object tracking [START_REF] Poore | Some assignment problems arising from multiple target tracking[END_REF].

First, let us consider the prediction stage at time k in which we are given a filtering density evaluated at the most likely data association hypothesis τ j k´1 of the previous step 5 . The latter is a product of its marginals ppX k´1 , τ j k´1 " τ j k´1 |Z j 1:k´1 q " M ź m"1 p m px k´1,m , τ j k´1 " τ j k´1 |Z j 1:k´1 q where p m px k´1,m , τ j k´1 " τ j k´1 |Z j 1:k´1 q " ppx k´1,m |z j,m 1:k´1 q. Here, z i,m 1:k´1 denotes the measurements induced by object m from step 1 to k ´1, i.e.,

z j,m 1:k´1 fi ˆzj k´1,ρ j k´1 pmq , z j k´2,ρ j k´2 pmq , . . . , z j 1,ρ j 1 pmq ˙,
where ρ is the inverse of τ , i.e., ρ˝τ is the identity permutation. Consequently, the mth marginal of the posterior at k ´1 is a Gaussian density that can be obtained equivalently by Kalman filtering [START_REF] Ristic | Beyond the Kalman Filter; particle filters for tracking applications[END_REF] over z j,m 1:k´1 , i.e., p m px k´1,m |z j,m 1:k´1 q " N px k´1,m ; xj k´1,m , P j k´1,m q, (49) with the mean and covariance matrices over time specifying the mth "track."

Hence, the prediction density in (39) (Section V-A) evaluated at τ j k´1 " τ j k´1 for the state transition in ( 47) is given by

ppX k |Z j 1:k´1 q " M ź m"1
N px k,m ; xj k|k´1,m , P j k|k´1,m q (50) xj k|k´1,m " Fx j k´1,m , P j k|k´1,m " FP j k´1,m F T `Q, where the last two lines are Kalman prediction equations with p.q T denoting matrix transpose.

Next, let us consider the update stage in which we use the prediction density [START_REF] Silverman | Density Estimation for Statistics and Data Analysis[END_REF] with M measurements concatenated in Z j k and the measurement likelihood. This likelihood is found by substituting [START_REF] Bertsekas | The auction algorithm: A distributed relaxation method for the assignment problem[END_REF] in [START_REF] Lin | Multisensor multitarget bias estimation for general asynchronous sensors[END_REF]. First, we use this term within the likelihood for the association variable τ j k which -as it is mutually independent from τ j 1 , . . . , τ j k´1 -is given by

l j pZ j 1:k |τ j k q 9 p τ j k pZ j k |Z j 1:k´1 q (51)
"

ż lpZ j k |X k , τ j k qppX k |Z j 1:k´1 qdX k " M ź o"1 ż g j pz j k,o |x 1 k qp τ j k poq px 1 k |Z j 1:k´1 qdx 1 k " M ź o"1 ż g j pz j k,o |x 1 k qppx 1 k |z j,τ j k poq 1:k´1 qdx 1 k .
The first line above follows from that the joint distribution of Z j 1:k´1 and τ j k is independent of the latter. The prior distribution for τ j k is non-informative as given in [START_REF] Ristic | Calibration of multi-target tracking algorithms using non-cooperative targets[END_REF], so, the ML estimate using (51) coincides with the MAP estimate and it is given by τ j k " arg max

τ j k PSM l j pZ j 1:k |τ j k q (52)
where S M is the set of M -permutations. An equivalent problem is found by taking the logarithm of the objective function in the combinatorial optimisation problem above as follows:

τ j k " arg max τ j k PSM M ÿ o"1 c `o, m " τ j k poq ˘(53) cpo, mq fi log ż g j pz j k,o |x 1 k qppx 1 k |z j,m 1:k´1 qdx 1 k ,
for o, m " 1, . . . , M . This cost for the LGSS model is explicitly found using the prediction distribution [START_REF] Silverman | Density Estimation for Statistics and Data Analysis[END_REF] and the measurements within the KF innovations [START_REF] Ristic | Beyond the Kalman Filter; particle filters for tracking applications[END_REF] as 53) is a 2 ´D assignment problem which can be solved in polynomial time with M (despite that the search space S M has a factorial size) using one of the well known solvers [START_REF] Burkard | Assignment Problems[END_REF] including the auction algorithm [START_REF] Bertsekas | The auction algorithm: A distributed relaxation method for the assignment problem[END_REF]. This algorithm operates over a matrix of costs obtained by C " rcpo, mqs to iteratively find the M pairs corresponding to the best permutation τ j k in the ML problem [START_REF] Vo | The gaussian mixture probability hypothesis density filter[END_REF]. Here, com-putation of the M 2 cost matrix usually has the predominant computational time.

cpo, mq " log N pz j k,o ; ẑj k,m , S j k,m q (54) ẑj k,m " H j xj k|k´1,m , S j k,m " R j `Hj P j k|k´1,m H T j . The optimisation in (
Next, we consider the state distribution update (see [START_REF] Andrieu | Particle markov chain monte carlo methods[END_REF] in Section II-B) and assert the empirical (model) prior in [START_REF] Ala-Luhtala | An introduction to twisted particle filters and parameter estimation in non-linear statespace models[END_REF]. As a result, the filtering density at k becomes a product of its marginals each of which is a Gaussian as in [START_REF] Casella | Monte Carlo Statistical Methods[END_REF] found by the KF update [START_REF] Ristic | Beyond the Kalman Filter; particle filters for tracking applications[END_REF], i.e.,

ppX k , τ j k " τ j k |Z j 1:k q " M ź m"1 N px k,m ; xk,m , P k,m q (55)
where,

xk,m " xk|k´1,m `Kk,m pz j k,ρ j k pmq ´Hj xk|k´1,m q P k,m " pI ´Kk,m H j q P j k|k´1,m K k,m fi P j k|k´1,m H T j S j k ´1 ,m .
Note that, because the posterior density is now non-zero only for τ j k " τ j k , the marginalisation over τ j k (see, e.g., ( 14)) in the following prediction stage reduces to [START_REF] Ihler | Nonparametric belief propagation for self-localization of sensor networks[END_REF] and [START_REF] Silverman | Density Estimation for Statistics and Data Analysis[END_REF].

B. Evaluation of the calibration quad-term in the LGSS Model

Let us consider the evaluation of s k j given by ( 42) using the formulae for the LGSS model introduced in Section VI-A. By comparison with the cost term in ( 53) and (54), it can easily be seen that s k j is the exponential of the association cost for τ j k , i.e.,

s k j pZ j k q " exp M ÿ o"1 c `o, τ j k poq ˘. (56)
Next, let us consider [START_REF] Yedidia | Constructing free-energy approximations and generalized belief propagation algorithms[END_REF] for evaluating r k ij . Evaluation of this term involves finding the object identity correspondance γk by solving a 2-D assignment as explained in Appendix E. In the LGSS model, the assignment cost matrix D " rdpo, mqs is found as

dpo, mq " log N pz i k,o ; ẑi k,m , S i k,m q (57) ẑi k,m " H i T θ px j k,m q, S i k,m " R i `Hi T θ pP j k,m qH T i .
for o, m " 1, . . . , M . Here, the second order statistics P j k,m is also transformed by applying any rotations involved in T θ to its eigenvectors. The best assignment which here encodes γk is found using the auction algorithm [START_REF] Bertsekas | The auction algorithm: A distributed relaxation method for the assignment problem[END_REF] as well, similar to the assignment in the Bayesian filtering update (Sec. VI-A). Using this estimate, the quad-term factor is computed using

r k ij pZ i k , θq " exp M ÿ o"1 d ´o, γk `τ i k poq ˘¯. (58)
In order to evaluate the other factors of the quad-term update, i.e., s i k and r k ji , similar computations are used. It suffices to replace i in the subscripts/superscripts of the expressions above with j, and, vice versa.

Finally, let us consider the scale factor in [START_REF] Huang | Fourier theoretic probabilistic inference over permutations[END_REF]. Let us use the notation introduced in the previous section for expressing the densities inside the integration. Starting with [START_REF] Üney | Distributed localisation of sensors with partially overlapping field-of-views in fusion networks[END_REF], one obtains

p o pz i , z j |Z i 1:k´1 , θq " ppz i , z j |z i,τ i k´1 poq 1:k´1 , θq (59)
" N prz i , z j s T ; µ 1 , Σ 1 q µ 1 " « ẑi k,m H j T ´1 θ px i k,m q ff Σ 1 " « S i k,m 0 0 R j `Hj T ´1 θ pP i k,m qH T j ff
where m " τ i k´1 poq, and, ẑi k,m and S i k,m are computed as in (54) with i substituted in place of j.

The second density in [START_REF] Huang | Fourier theoretic probabilistic inference over permutations[END_REF] conditioned on sensor j's history, i.e., [START_REF] Guibas | The identity management problem-a short survey[END_REF], is similarly found as

p o pz i , z j |Z j 1:k´1 , θq " N prz i , z j s T ; µ 2 , Σ 2 q µ 2 " « H i T θ px j k,m q ẑj k,m ff Σ 2 " « R i `Hi T θ pP j k,m qH T i 0 0 S j k,m
ff where m " γk´1 ˝τ i k´1 poq and, ẑj k,m and S j k,m are given in (54).

Using the densities above and integration rules for Gaussians, the oth term of the scale factor in [START_REF] Huang | Fourier theoretic probabilistic inference over permutations[END_REF] is found as

κ k,o pθq " `ˇΣ ´1 1 ˇˇˇˇΣ ´1 2 ˇˇ˘1 {4 ˇˇΣ ´1 1 `Σ´1 2 2 ˇˇ1 {2 exp " ´1 4 `µT 1 Σ ´1 1 µ 1 `µT 2 Σ ´1 2 µ 2 1 4 `Σ´1 1 µ 1 `Σ´1 2 µ 2 ˘T `Σ´1 1 `Σ´1 2 ˘´1 ˆ`Σ ´1 1 µ 1 `Σ´1 2 µ 2 ˘( . (60
) In a distributed setting, the scale factor expressions above are computed both at sensors i and j, for which the prediction stage given in ( 50) is carried out for both the local posterior and the posterior recevied from the other sensor at time k ´1.

C. Sampling from the calibration marginals using nonparametric BP

In this section, we introduce particle based representations and Monte Carlo computations [START_REF] Casella | Monte Carlo Statistical Methods[END_REF] for the realisation of (loopy) BP message passings. Note that Sections VI-A and VI-B specify the evaluation of edge potentials given in [START_REF] Uney | Distributed estimation of latent parameters in state space models using separable likelihoods[END_REF] for pτ i 1:t , τ j 1:t q " pτ i 1:t , τ j 1:t q. For sampling from the marginal parameter posteriors, we adopt the approach detailed in [26, Sec.VI] for carrying out LBP belief update and messaging in [START_REF] Uney | Cooperative sensor localisation in distributed fusion networks by exploiting non-cooperative targets[END_REF] and [START_REF] Mizrahi | Distributed parameter estimation in probabilistic graphical models[END_REF], respectively. Given L equally weighted samples from pi pθ i q, i.e.,

θ plq i " pi pθ i q, (61)
for l " 1, . . . , L, the edge potentials are evaluated to obtain

ψ ij pθ plq i , θ plq j q " l ´Zi 1:t , Z j 1:t |θ " pθ plq i , θ plq j q ¯. ( 62
)
Consider the BP message from node j to i in [START_REF] Mizrahi | Distributed parameter estimation in probabilistic graphical models[END_REF]. Suppose that independent identically distributed (i.i.d.) samples from the (scaled) product of the jth local belief and the incoming messages from all neighbours except i are given, i.e., θplq j " pj pθ j q ź i 1 Pnepjq{i m i 1 j pθ j q for l " 1, ..., L.

These samples are used with kernel approximations in order to represent the message from node j to i (scaled to one), in the NBP approach [START_REF] Sudderth | Nonparametric belief propagation[END_REF]. We use Gaussian kernels leading to the approximation given by

mji pθ i q " L ÿ l"1 ω plq ji N pθ i ; θ plq ji , Λ ji q, (64)
θ plq ji " T pT ´1p θplq j ; θ plq j q; θ plq i q, ω plq ji " ψ i,j pθ plq i , θ plq j q ř L l 1 "1 ψ i,j pθ pl 1 q i , θ pl 1 q j q
, where the kernel weights are the normalised edge potentials. Λ ji is related to a bandwidth parameter that can be found using Kernel Density Estimation (KDE) techniques. In particular, we use the rule-of-thumb method in [START_REF] Silverman | Density Estimation for Statistics and Data Analysis[END_REF] and find

Λ ji " ˆ4 p2d `1qL ˙2{pd`4q Ĉji , Ĉji " ÿ l 1 ÿ l ω pl 1 q ji ω plq ji pθ pl 1 q ji ´m ji qpθ plq ji ´m ji q T , mji " L ÿ l"1 ω plq ji θ plq ji
where mji and Ĉji are the empirical mean and covariance of the samples, respectively, and d is the dimensionality of θ ji s. Given these messages, let us consider sampling from the updated marginal in [START_REF] Uney | Cooperative sensor localisation in distributed fusion networks by exploiting non-cooperative targets[END_REF]. We use the weighted bootstrap (also known as sampling/importance resampling) [START_REF] Smith | Bayes statistics without tears: A sampling-resampling perspective[END_REF] with samples generated from the (scaled) product of Gaussian densities with mean and covariance found as the empirical mean and covariance of the particle sets, respectively. In other words, given mji and Ĉji as above, we generate θ plq i " f pθ i q, l " 1, . . . , L, f pθ i q 9 N pθ i ; mi , Ĉi q ź jPnepiq N pθ i ; mji , Ĉji q.

The particle weights for these samples to represent the updated marginal is given by

ω plq i " ωplq i { L ÿ l 1 "1 ωpl 1 q i ωplq i " ´p0,i pθ plq i q ź jPnepiq mji pθ plq i q ¯{f pθ plq i q
where p 0,i is the prior density selected for θ i (and, the node potential in [START_REF] Wiesel | Distributed covariance estimation in gaussian graphical models[END_REF]). Thus, the local calibration marginal is estimated by

Pi pdθ i q " L ÿ l"1 ω plq i δ θ plq i pdθ i q. ( 65
)
Algorithm 1 Pseudo-code for estimation of θ using the quadterm separable likelihood within Belief Propagation.

1: for all j P V do Ź Local filtering 2: for k " 1, . . . , t do

3:

Find ppX k , τ j k " τ j k |Z j 1:k q in (55) as described in Section VI-A 4:

Find s k j pZ j k q in (56) 5: end for 6: end for 7: for all j P V do Ź Sample from priors 8: Sample θ plq i " p 0,i pθ i q for l " 1, . . . , L as in (61) 9: end for As the final step of the bootstrap, tθ plq i , ω plq i u M l"1 is resampled (with replacement) leading to equally weighted particles from pi pθ i q, i.e., tθ plq i u L l"1 . We follow similar bootstrap steps in order to generate the samples in (63).

After nodes iterate the BP computations described above for S times, each node estimates its location by finding the empirical mean of tθ plq i u L l"1 . These steps are summarised in Algorithm 1.

VII. EXAMPLE: SELF-LOCALISATION IN LGSS MODELS

In this example, we demonstrate the quad-term node-wise separable likelihood in sensor self-localisation. The LGSS model given by ( 47) and ( 48) is used with process noise parameters selected as σ " 0.5, q 1 " 1{4, q 2 " q 3 " 1{2, q 4 " 1. The measurement model for sensor i is given by H i " rI, 0s and R i " σ 2 n I with σ n " 10 modelling noisy position measurements in the local coordinate frame.

Let us consider the multi-object multi-sensor scenario depicted in Fig. 2. 16 sensors observe 4 objects moving with data association uncertainties. The locations of the sensors are to be estimated with respect to sensor 1 which is selected as the origin of the network coordinate system. Therefore θ " rθ 1 , . . . , θ 16 s with the prior distribution for θ 1 selected as Dirac's delta, i.e., p 0,1 pθ 1 q " δpθ 1 q. For the other nodes, the localisation prior, i.e., p 0,i pθ i q for i " 2, . . . , 16, is a uniform distribution over the sensing region.

We use Algorithm 1 specified in Section VI for estimating θ. The MRF model we consider is specified by the pairwise graph G in Fig. 2 (blue edges). We use L " 100 points in (61) to represent the local belief densities. We select the sensor data time window length as t " 10 (starting at time 21 until 30 in the scenario in Fig. 2). We follow the steps in Algorithm 1 for S " 16 iterations.

A typical run is illustrated in Fig. 3. Here, the scatter plot of particles from marginal posteriors are given over iterations. Note that the network coordinate system is established by sensor 1 through its informative prior, and, LBP emanates this information towards the outer nodes while learning the edge potentials using node-wise separable likelihood evaluations ψ i,j pθ plq i , θ plq j q given in ( 23) and [START_REF] Uney | Distributed estimation of latent parameters in state space models using separable likelihoods[END_REF]. For performance assesment, first, we consider the mean squared error for θ output by our algorithm, as we have built our discussion on MMSE estimators in [START_REF] Kantas | On particle methods for parameter estimation in state-space models[END_REF]. We find this value empirically by taking the average of the squared norm of estimation errors over 100 Monte Carlo simulations. In Fig. 4, we present a semi-log plot of this quantity over iterations (blue line). Note that, convergence occurs in less than ten iterations which is a favourable feature. We compare this algorithm with the RFS based dual-term pseudo-likelihood proposed in [START_REF] Uney | A cooperative approach to sensor localisation in distributed fusion networks[END_REF]. When evaluating this term, we use a Poisson multi-object model output by using the Gaussian mixture probability hypothesis density (GM-PHD) filter [START_REF] Vo | The gaussian mixture probability hypothesis density filter[END_REF] with the LGSS model. The averaged MSE performance for the case that dual-term likelihoods are used as edge potentials is depicted with the green dashed line in Fig. 4. The quad-term approximation is seen to provide faster convergence with both pseudo-likelihoods leading to an on par accuracy in the steady regime, in this example. The edge update time for the quadterm update averages to 0.601 per edge per particle compared to 1.312 for the dual-term update demonstrating its relative efficiency. Note that, the MSE is a network-wide term and the local error norms are smaller. The localisation miss-distance averaged over sensors is given in Fig. 5. The average error (ȏ ne standard deviation) in the final step is 2.60 ˘0.70m with a maximum value of 4.96 which is less than 0.5% of the edge distances of 1000m. These results demonstrate that the proposed scheme is capable of providing self-localisation with favourable accuracy and small error margins.

VIII. CONCLUSIONS

In this work, we have addressed the prohibitive complexity of latent parameter estimation in state space models when there are many sensors collecting measurements. We proposed a pseudo-likelihood, namely the quad-term node-wise separable likelihood, as an accurate surrogate to the actual likelihood which is extremely costly to evaluate. The separable structure of this quad-term approximation makes it possible to evaluate it using local filtering operations, hence, scale with the number of sensors.

In order to use the proposed approximation in the case of multiple objects, we employed a parameterised multiobject state space model in which different configurations of the parameter specify different hypothesis of object-tomeasurement and object-to-object associations. Specifically, we introduced an empirical Bayesian perspective for evaluating separable likelihoods in this model, using only local Bayesian filtering. This approach substitutes the estimates of the hypothesis variables when evaluating the likelihood of the latent parameter, and, decouples these two inference tasks if the hypothesis variables can be estimated locally. Therefore, it can be extended to general hypothesis variables that capture, for example, variable number of objects M k , less than one probability of detection, i.e., P D ă 1, and, measurements with false alarm, by adapting the assignment problems involved accordingly (see, for example [START_REF] Poore | Some assignment problems arising from multiple target tracking[END_REF]).

The associated posterior distribution is a MRF over which distributed inference is possible using message passing algorithms such as LBP. We specified a particle message passing algorithm for sampling from latent parameter marginals for a linear Gaussian state space model with multiple objects. This algorithm is demonstrated in simulations for sensor selflocalisation using point measurements from non-cooperative objects in a potentially GPS denying environment. It is possible to estimate unkown orientation angles as well, by appropriately defining the transform in [START_REF] Kokkala | Combining particle MCMC with Rao-Blackwellized Monte Carlo data association for parameter estimation in multiple target tracking[END_REF] and taking into account the rotations implied by θ i s in the LBP steps of Algorithm 1 (additional details can be found in [START_REF] Üney | Enabling self-configuration of fusion networks via scalable opportunistic sensor calibration[END_REF]).

APPENDIX

A. Empirical Bayes parameter likelihood

The empirical Bayes parameter posterior follows from the decomposition of the posterior in (11) using the chain rule of probabilities as ppθ|Z 1 1:t , . . . , Z N 1:t q " (66) ÿ 1:t , . . . , Z N 1:t , τ 1:N 1:t qppτ 1:N 1:t |Z 1 1:t , . . . , Z N 1:t q.

The first term inside the summations is the parameter posterior conditioned on the association variables and the second term is similar to a prior with respect to the first term. The fact that this term is conditioned on the measurements makes an empirical selection possible as discussed in Section II-B. Let us use a similar empirical prior selection approach as used in [START_REF] Ala-Luhtala | An introduction to twisted particle filters and parameter estimation in non-linear statespace models[END_REF], i.e., ppτ [START_REF] Cover | Elements of Information Theory[END_REF] suggest that ǫ ě 0 regardless of the problem setting (as conditioning reduces entropy). Under normal sensing conditions it is reasonable to expect that ǫ ą 0 holds as sensor measurements are highly informative on all the variables at time k. For sensors with similar sensing capabilities, it is also reasonable to expect that their measurement history are interchangeable. In other words, HpZ i k |Z j 1:k´1 , Θq `HpZ j k |Z i 1:k´1 , Θq " HpZ j k |Z j 1:k´1 , Θq `HpZ i k |Z i 1:k´1 , Θq (81) holds in which case (80) is nonzero and Dpp||uq ą Dpp||qq. This condition can be relaxed for a difference of ǫ between the RHS and LHS of (81).

Comparison in terms of entropy upper bounds is more straightforward. Let us consider [START_REF] Sudderth | Nonparametric belief propagation[END_REF] and Corollary 4.2 in [START_REF] Uney | A cooperative approach to sensor localisation in distributed fusion networks[END_REF] which we repeat here for convenience:

Dpp||uq ď H ù ´Hú , (82)
H ù " HpX k |Z j 1:k´1 , Θq `HpX k |Z i 1:k´1 , Θq, H ú " HpX k |Z i 1:k´1 , Z j 1:k´1 , Θq maxtHpX k |Z i 1:k´1 , Z j 1:k , Θq, HpX k |Z j 1:k´1 , Z i 1:k , Θqu. For sensors of identical capabilities, both terms in the maximisation should be identical as the accuracy of the state estimate should not differ for using either of the sensors' current measurement in addition to the histories of both. As a result, we can replace H ú with

H ú " HpX k |Z i 1:k´1 , Z j 1:k´1 , Θq1 2 ´HpX k |Z i 1:k´1 , Z j 1:k , Θq `HpX k |Z j 1:k´1 , Z i 1:k , Θq ¯.
Now, note that the sum of the positive weighted terms in the RHS of (32) (let us denote by H q ) is smaller than H ù owing to that conditioning reduces entropy [START_REF] Cover | Elements of Information Theory[END_REF]. The sum of the negative weighted terms (let us denote by H q ), for the case ,equals to H ú , leading to H ù ´Hú ´pH q ´Hq q ą 0, which indicates that the entropy bound of the quad-term update is smaller than that for the dual-term update.

E. Local ML estimate of object correspondances

The semantic of the object correspondance γ k is slightly different from that of the random permutation variables used in identity management [START_REF] Guibas | The identity management problem-a short survey[END_REF] in a way closer to data association because local (track) identities are synonymous with the measurements they are associated with instead of signal features etc. Therefore, in our model, γ k is uniquely defined when given τ i k and τ j k , and, γ 1 , . . . , γ k are mutually independent. The ML estimate of γ k using the likelihood for data set pZ i k , Z j 1:k q is also conditioned on the association configuration pτ i k , τ i k q and θ is given by γk " arg max

γ k PSM log lpZ i k , Z j 1:k |γ k , τ i k " τ i k , τ j k " τ j k , θq (83)
where S M is the set of M -permutations. This likelihood can be decomposed using the chain rule of probability as

lpZ i k ,
where the second equality follows from the independence of sensor j's measurements up to time k ´1 from the association variables at time k. The first term on the RHS is easily identified as (43) evaluated for γ k when m 1 " γ k pτ i k poqq. This term, when subsituted in (83), leads to a 2-D assignment problem given by γk " arg max

γ k PSM M ÿ o"1 dpo, m 1 " γ k pτ i k poqqq
dpo, m 1 q fi log r k ij,o pz i k,o , θq, for o, m 1 " 1, . . . , M , where r k ij,o is given in [START_REF] Yedidia | Constructing free-energy approximations and generalized belief propagation algorithms[END_REF]. After finding the M 2 costs above, this problem can be solved using the auction algorithm [START_REF] Bertsekas | The auction algorithm: A distributed relaxation method for the assignment problem[END_REF] in polynomial time with M . This algorithm finds the M pairs corresponding to the best permutation γk . We use a similar approach for the data association problem in Bayesian filtering in Section VI-A.

Fig. 1 .

 1 Fig. 1. A multi-sensor state space -or, hidden Markov-model (black dashed box on the right representing a chain over k) and a Markov Random field model of the parameter posterior (the blue edges on the left).

Fig. 2 .

 2 Fig. 2. Example scenario: 16 Sensors collect measurements from 4 objects (T1-T4) with association uncertainties. Initial positions of the objects are denoted by black squares. Trajectories for 60 time steps are depicted. The blue lines depict the edges of the MRF model used for estimation.

Fig. 3 .

 3 Fig. 3. Node beliefs in LBP iterations: Marginal posterior estimates after iteration 4 (upper-left), 6 (upper-right), 8 (lower-left), and, 10 (lower-right).

Fig. 4 .

 4 Fig. 4. Log-normalised error margin versus the iteration number n.

Fig. 5 .

 5 Fig. 5. Localisation miss-distance averaged over nodes versus the iteration number. 100 Monte Carlo runs displayed with the boxes centered at the median (red). Edges (blue) indicate the 25th and 75th percentiles.

  Now, let us compare (79) and (78): The negative weighted terms are equal, so, the difference of the positive weighted terms are considered: ∆u `´∆q `" HpZ j k |Z i 1:k´1 , Θq ´HpZ j k |Z j 1:k´1 , Θq (80)

	`HpZ j k |Z i 1:k , Θq `HpZ i k |Z j 1:k , Θq ı
	`HpZ i k |Z j 1:k´1 , Θq ´HpZ i k |Z i 1:k´1 , Θq	`ǫ
	ǫ fi HpZ j k |Z i 1:k´1 , Θq ´HpZ j k |Z i 1:k , Θq
	`HpZ i k |Z j 1:k´1 , Θq ´HpZ i k |Z j 1:k , Θq
	Properties of differential entropy

1:N 1:t |Z 1 1:t , . . . , Z N 1:t q " t ź k"1 ppτ 1:N k |Z 1 k , . . . , Z N k q ppτ 1:N k |Z 1 k , . . . , Z N k q Ð δ τ 1:N k´1 pτ 1:N k´1 q. (

67

)

After substituting from (67) in (66), the parameter posterior is found as

Specifically, for linear Gaussian measurements with no data association uncertainty, the marginal parameter likelihood involves computation of the innovation covariance for the so called group-sensor measurements in joint multi-sensor filtering.

Note that the results presented in this section are valid for any selection of τ i,j 1:t as they relate random variables which are conditioned on the data association. The divergences and bounds, nevertheless, are more relevant for τ i,j 1:t " τ i,j 1:t .

Note that, when the calibration parameters also include orientation angles, T θ p.q involves rotation matrices accordingly.

It can easily be shown that this state transition model is invariant under the selected coordinate frame for x k s.

Note that the empirical prior on τ j k´1 is selected as in[START_REF] Ala-Luhtala | An introduction to twisted particle filters and parameter estimation in non-linear statespace models[END_REF] leading to this posterior be identically zero for all values of τ j k´1 other than τ j k´1 .

ACKNOWLEDGMENT

The authors would like to thank Dr. Sinan Yildirim for discussions on the general multi-object tracking model in [14] and [27], and, Prof. Arnaud Doucet for discussions on pseudo-likelihood methods in state space models.

This work was supported by the Engineering and Physical Sciences Research Council (EPSRC) Grant number EP/K014277/1 and the MOD University Defence Research Collaboration (UDRC) in Signal Processing. Murat Üney and Bernard Mulgrew are with the Institute for

ppθ|Z 1 1:t , . . . , Z N 1:t q " ppθ|Z 1 1:t , . . . , Z N 1:t , τ 1:N 1:t " τ 1:N 1:t q 9 lpZ 1 1:t , . . . , Z N 1:t |θ, τ 1:N 1:t " τ 1:N 1:t qppθq where the likelihood in the last line is given by [START_REF] Marin | Approximate bayesian computational methods[END_REF] and [START_REF] Varin | An overview of composite likelihood methods[END_REF] in Section II-B.

B. Proof of Proposition 4.1

Proof. Let us expand the KLD term in [START_REF] Cover | Elements of Information Theory[END_REF] by substituting its arguments given in [START_REF] Kantas | Distributed maximum likelihood for simultaneous self-localization and tracking in sensor networks[END_REF], [START_REF]Distributed Data Fusion for Network-Centric Operations[END_REF] and [START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF]:

In the equations above, κ k pθq is a normalisation constant given by [START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF]. Eq.( 68) is obtained after multiplying both the numerator and the denominator of the quotient inside the logarithm by ppZ i 1:k´1 |Z j 1:k´1 , θqppZ j 1:k´1 |Z i 1:k´1 , θq and a rearrangement of the terms. The definition of MI [START_REF] Cover | Elements of Information Theory[END_REF] results with the first two terms in Eq.(69). The last term is the expectation of the normalisation constant over the joint distribution of the sensor histories Z i 1:k´1 and Z j 1:k´1 , and, Θ. Let us now consider the normalisation constant:

" 1.

The inequality (70) follows from Hölder's Inequality. Consequently, the last term in (69) is non-positive, and, (31) is obtained.

C. Proof of Corollary 4.2

Proof. We apply the chain rule of information to the MI terms on the RHS of (31) leading to

and,

The MI terms on the RHSs of the equations above are for random variables which form Markov chains with the current the state variable X k . Consider the (conditional) chains

k´1 , Θ for the RHS of Eq.( 71). The Data Processing Inequality [START_REF] Cover | Elements of Information Theory[END_REF] applied to these terms lead to

`HpX k |Z j k , Z j 1:k´1 , Θq ´HpX k |Z j k , Z j 1:k´1 , Z i 1:k´1 , Θq A similar break down of Eq.( 72) results with 73) and ( 74) into (31) results with [START_REF] Sudderth | Nonparametric belief propagation[END_REF].

D. Comparison of the quad-term and dual-term updates

Let us compare the KLDs of the quad-term and dual-term updates. The KLD of the dual-term in ( 36) is given by [START_REF] Uney | A cooperative approach to sensor localisation in distributed fusion networks[END_REF]

The MI terms on the RHS of the above equation can be expanded using that IpA; B|Cq " HpA|Cq ´HpA|C, Bq [START_REF] Cover | Elements of Information Theory[END_REF]. The third term can be expanded in alternative ways as follows:

"

77) After decomposing the first two terms on the RHS of (75) similarly and adding to the average of ( 76