Asma Rabaoui
email: asma.rabaoui@ec-lille.fr

Nicolas Viandier

Juliette Marais

Emmanuel Duflos

DPMs for the Density Estimation in a Dynamic Nonlinear Modeling: Application to GPS Positioning in Urban Canyons

Keywords: Global Positioning System (GPS), navigation, urban canyon, pseudorange errors, density estimation, nonparametric Bayesian Methods, sequential Monte Carlo methods, Rao-Blackwellized particle filter

In Global Positioning System (GPS), classical localization algorithms assume, when the signal is received from the satellite in line-of-sight (LOS) environment, that the pseudorange error distribution is Gaussian. Such assumption is in some way very restrictive since a random error in the pseudorange measure with an unknown distribution form is always induced in constrained environments especially in urban canyons due to multipath/masking effects. In order to ensure high accuracy positioning a good estimation of the observation error in these cases is required. To address this, an attractive flexible Bayesian nonparametric noise model based on Dirichlet process mixtures (DPM) is introduced. Since the considered positioning problem involves elements of non-Gaussianity and non-linearity and besides it should be processed on-line, the suitability of the proposed modeling scheme in a joint state/parameter estimation problem is handled by an efficient particle filter called Rao-Blackwellized particle filter (RBPF).

Our approach is illustrated on a data analysis task dealing with joint estimation of vehicles positions and pseudorange errors in a GNSS based localization context where the GPS information may be inaccurate because of multipath/masking effects.

I. INTRODUCTION

GPS is a radio navigation system which relies on radio-frequency signals emitted by a constellation of satellites.

Consequently, it can be easily distorted by the measuring system and propagation impairments (atmospheric layers, multipath, diffraction and mask phenomena). Today, Global Navigation Satellite Systems (GNSS) have penetrated the transport field through a variety of applications such as monitoring of containers, fleet management or car navigation. These applications do not necessarily request high availability, integrity and accuracy of the positioning system. However, for safety applications (such as guidance of autonomous vehicles), performances require to be much more stringent.

Unfortunately, most of all these transport applications are mainly used in dense urban environments. Even with modern GPS receivers, high positioning accuracy is only achieved in LOS conditions when the signal is received directly without any reflexion. Since the signals delivered by sensors are easily blurred because of hard external conditions in urban canyons, they may deliver totally erroneous measurements, especially when the signal reaches the receiver after interaction with several objects/obstacles. In urban areas, some obstacles (buildings, trees, etc) can be modeled using a simulator and thus some of the errors can be predicted or modeled. However, some other obstacles (cars, pedestrians, etc) can appear suddenly and can induce a random error in the pseudorange measure.

In order to ensure high accuracy positioning, a good estimation of the observation error in these cases is required.

The incoming signal is in general the sum of the direct signal and several delayed replica. However, non-LOS signals, which are signals received after reflections on the surrounding obstacles, frequently occur in dense environments. In these environments, most of the time no direct signals can reach the receiver and thus localization accuracy is severely degraded due to high propagation time delays.

In our paper, we consider positioning accuracy degradation problem when we don not receive the direct signal but just reflected replica. In this case, the receiver tracks only reflected signals. Such phenomena make the pseudorange error distribution becomes a non-centered Gaussian distribution because of additional errors on pseudorange estimation. As a consequence, classical localization methods assuming that state and observation noises are zero mean Gaussian are not efficient anymore and can severely impair positioning accuracy.

A. Previous Work

Literature focusing on techniques for localization performance enhancement in constrained environments is abundant. The most spread rely on multi-sensors approaches [START_REF] Giremus | A particle filter approach for joint detection/estimation of multipath effects on gps measurements[END_REF], [START_REF] Caron | Particle filtering for multisensor data fusion with switching observation models. Application to land vehicle positioning[END_REF], which goals are to compensate the lack of GNSS performances by adding other sensors (odometer, inertial measurement unit, etc). However multi-sensors methods show some drawbacks and in particular the system cost and complexity.

In [START_REF] Caron | Particle filtering for multisensor data fusion with switching observation models. Application to land vehicle positioning[END_REF], [START_REF] Nahimana | Reception state estimation of gnss satellites in urban environment using particle filtering[END_REF], [START_REF] Viandier | GNSS performance enhancement in urban environment based on pseudo-range error model[END_REF], the proposed approaches allow to adapt the error model in the filtering process to the reception condition of each satellite signal. This modeling permits to handle the overall reception process which switches between some preselected measurement models. The switching process permits to define three reception states (LOS, non-LOS and no reception cases). The principle idea of such modeling scheme is to define an indicator variable which governs the behavior of the statistical model and allow for switching from one model to another.

An interesting way of defining the statistical structure of this indicator variable has been developed in [START_REF] Caron | Particle filtering for multisensor data fusion with switching observation models. Application to land vehicle positioning[END_REF].

Switching state-space models, also called jump state-space models, have been widely studied in the literature.

In [START_REF] Bar-Shalom | Multitarget-multisensor tracking: principles and techniques[END_REF], [START_REF] Tugnait | Adaptive estimation and identification for discrete systems with markov jump parameters[END_REF] a Jump Markov Linear System (JMLS) has been proposed for the case where the considered state-space model is conditionally linear and Gaussian. The estimation problem was handled using finite mixture modeling (Interacting Multiple Model (IMM) [START_REF] Blom | The interacting multiple model algorithm for systems with markovian switching coefficients[END_REF], [START_REF] Mazor | Interacting multiple model methods in target tracking: a survey[END_REF] and Generalized Pseudo-Bayes (GPB) [START_REF] Bar-Shalom | Multitarget-multisensor tracking: principles and techniques[END_REF] algorithms). For non-linear and non-Gaussian cases, efficient sequential Monte Carlo algorithms have been applied [START_REF] Andrieu | Efficient particle filtering for jump markov systems. application to time-varying autoregressions[END_REF], [START_REF] Doucet | Recursive state estimation for multiple switching models with unknown transition probabilities[END_REF], [START_REF] Jilkov | Online bayesian estimation of transition probabilities for markovian jump systems[END_REF], [START_REF] Doucet | Particle filters for state estimation of jump markov linear systems[END_REF].

In [START_REF] Nahimana | Reception state estimation of gnss satellites in urban environment using particle filtering[END_REF] the proposed Jump Markov System (JMS) shows quite good results in terms of accuracy and stability in the context of GPS positioning. These acceptable performances are owed to the adaptation of the error model according to the reception state. However, the proposed approach for error estimation using the Expectation-Maximization (EM) algorithm has some limits especially related to slow convergence. Moreover, the non-LOS case is modeled by a finite Gaussian mixture and the number of Gaussian components is fixed. Consequently, the mixture model can not represent the true measurement error.

In this paper, in order to limit costs, we have chosen to work only with GNSS signals, no additive sensors will be used. For a better accuracy, we propose a new statistical filtering method based on a better definition (and use) of the observation noise for each satellite signal. In order to ensure high accuracy positioning a good estimation of the observation error in such cases is required.

B. The Motivating Problem

In this paper, the state evolution equation governing the position of a land vehicle is a dynamic non-linear model and the measurement equation (observation model) is nonlinear as well with additive observation errors. Within this setting, a joint estimation problem of hidden states and noise probability densities will be considered. Often, the observations which are the GPS measurements arrive sequentially in time and one is interested in performing the joint estimation on-line.

Note that in a dynamic nonlinear model with additive observation errors, it is usual to assume that errors are normally distributed or can be approximated by a finite Gaussian mixture. This can cause problems when there are for example outlying errors leading to many modes in the density distribution form whose variability over time induces poor inference about the model parameters. As we are considering in this paper the case where multipath propagation affects severely GPS signals qualities, it is our intention therefore to model the observation errors using a highly flexible family of density functions. It is important to notice that mixture modeling is considered as a successful and widely used density estimation method, capable of representing the phenomena that underlie many real-world datasets. However, the main difficulty in mixture analysis is how to choose the number of mixture components. Model selection methods in general treat the number of components as an unknown constant and set its value based on the observed data. Such an approach lacks flexibility, since in practice we often need to model the possibility that new observations come from as yet unseen components.

We propose to use a flexible density modeling based on a Bayesian nonparametric approach involving an infinite mixture model. The Dirichlet Process Mixture (DPM) model, studied in Bayesian nonparametrics [START_REF] Escobar | Bayesian density estimation and inference using mixtures[END_REF], [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF], is the most used approach for nonparametric density estimation. The construction of probability measures on the space of distribution functions has received considerable attention in the Bayesian nonparametrics literature [START_REF] Muller | Nonparametric Bayesian data analysis[END_REF]. A rich literature exists on its theoretical properties, and it has been used in a variety of problems [START_REF] Gelfand | Bayesian nonparametric Spatial Modeling with Dirichlet process mixture mixing[END_REF], [START_REF] Ishwaran | Approximate Dirichlet process computing in finite normal mixtures: smoothing and prior information[END_REF], [START_REF] Dorazio | On selecting a prior for the precision parameter of Dirichlet process mixture models[END_REF]. Despite this activity, only a few work has been conducted on how to tune the hyperparameters of a nonparametric mixture model. Most of the time, in real-world application contexts, the observed data can help in defining a prior for some of the hyperparameters in the considered Bayesian model. Computing every unknown hyperparameter with pure sampling is not the most effective approach in all cases since this can be practically difficult and time consuming.

Our objective is to estimate jointly the hidden states {x t ; t ∈ N} and the error density parameters that pertain to the observation noise vectors {w t ; t ∈ N}. These density parameters are the latent variables {θ t ; t ∈ N} to be estimated using a nonparametric Bayesian approach. For on-line Bayesian filtering, the particle filter permits to compute the posterior distribution p(x t , θ t |z 1:t ) recursively over time. In this case, we have to take into consideration the fact that computing everything with pure sampling can be computationally expensive. One solution consists in considering the augmented hidden state vector {Y t = (x t , θ t ); t ∈ N}. By the partition of the state-space Y t into two sub-spaces drawn by θ t and x t (for all t ≥ 1), improved formulations of the particle filter can be applied.

Rao-Blackwellization [START_REF] Doucet | Sequential Monte Carlo methods in practice[END_REF], [START_REF] Ristic | Beyond the Kalman Filter[END_REF] is one way of improving the efficiency of the particle filter. In the case where it is possible to evaluate some of the filtering equations analytically and the others with Monte Carlo sampling, Rao-Blackwellization will lead to estimators with less variance than what could be obtained with pure Monte Carlo sampling [START_REF] Casella | Rao-Blackwellisation of sampling schemes[END_REF] The idea is to partition the state vector so that one component of the partition is a conditionally nonlinear Gaussian state-space model; for this component one can work out the solution analytically and use the Extended Kalman filter to compute p(x 1:t |θ 1:t , z 1:t ). The particle filter is then used only for the nonlinear non-Gaussian portion of the state-space to compute p(θ 1:t |z 1:t ). In this way, we can say that the hard part of the problem will be reduced into the computation of the posterior distribution p(θ 1:t |z 1:t ).

C. Contributions and paper organization

This paper considers the problem of density estimation from a Bayesian nonparametric viewpoint, using a hierarchical mixture model. A typical choice is to assume the unknown density as a mixture of a parametric family with a discrete random probability as a mixing distribution. This random probability is induced by the Dirichlet Process and the hierarchical model is the Dirichlet Process Mixture model. By using DPM, it is important to mention that no matter what additive errors distributions are involved we are confident that our family of densities will be able to capture the right shape and hence statistical inference for the parameters of interest will be improved and reliable. A drawback of using a model based on the Dirichlet process is that it is infinite dimensional and therefore inference will be complicated. However, recent innovations in sampling algorithms within infinite dimensional frameworks has made considerable progress in recent years to such an extent that it is now possible to perform exact inference without the need to set up arbitrary approximations.

The flexibility of DPM models grows when we assume one more step in the hierarchy, i.e. when the the DPM hyperparameters are random. Our aim using this family of densities is first to be able to capture the right shape of the noise probability density functions (pdfs) and then to improve the statistical inference for the parameters of interest. Therefore, we will show how to tune the DPM hyperparameters in a flexible way. The posterior distribution some hyperparameters will be estimated as part of the Gibbs sampler, and some other parameters will be chosen carefully in a data-adaptive way for a better fitting of the data distribution shape.

To sum up, this paper proposes several contributions. The first concern is about the modeling of the observation noises using DPMs. Then, we focus on the suitability of this family of densities in the estimation problem handled by an improved particle filter called Rao-Blackwellized particle filter. The use of Rao-Blackwellization permits to improve the efficiency of the filtering process since we do not need to compute everything with pure sampling which reduces significantly the computational cost. An another contribution is about the application context and the interesting validation schemes. The efficiency of the proposed approaches is demonstrated by applying a validation step involving simulated GNSS signals.

The remainder of this paper is organized as follows. Section II presents the pseudorange error modeling. Section III is dedicated to the Bayesian modeling of the problem. Section IV is about the nonparametric density estimation problem. We show how to tune the DPM hyperparameters in a flexible way for a better fitting of the data distribution shape. Section V discusses the use of an efficient Particle filter to perform optimal estimation. In section VI the efficiency of the proposed approaches is demonstrated by conducting validation experiments involving simulated GNSS signals. Section VII concludes this paper with a summary and discussion.

II. PSEUDORANGE ERROR MODELING

In multi-sensors based systems, each sensor transmits a signal (an information) to the receiver (or antenna). In this paper, the GNSS (Global Navigation Satellite System) constellation [START_REF] Kaplan | understanding GPS: Priciples and Applications[END_REF] is considered as a sensor network.

Consequently, each GNSS satellite is considered as a sensor [START_REF] Rabaoui | Studies on dpm for the density estimation of pseudorange noises and evaluations on real data[END_REF]. These sensors work according to three different operation modes which are :

• Failure mode, when the sensor provide no information.

• Degraded mode, when the information delivered by the sensor is not accurate or is noisy.

• Normal mode, when the information delivered by the sensor can be considered as accurate.

The failure mode appears when the antenna cannot receive a signal from a satellite because of local masks. In the normal mode, the satellite signal reaches the antenna in LOS. Finally, the degraded mode occurs when the signal reaches the antenna after one or more reflections.

The navigation signal transmitted by each satellite includes a precise time at which the signal was transmitted.

The distance or range from a receiver to each satellite may be determined using this time of transmission which is included in each navigation signal. By noting the time at which the signal was received at the receiver, a propagation time delay can be calculated. This time delay when multiplied by the speed of propagation of the signal will yield a pseudorange from the transmitting satellite to the receiver.

The range is called a pseudorange because the receiver clock may not be precisely synchronized to GPS time and because propagation through the atmosphere introduces delays into the navigation signal propagation times. These result, respectively, in a clock bias and an atmospheric bias. Clock biases may be as large as several milliseconds.

Using information from at least four satellites, the position of a receiver with respect to the center of the Earth can be determined using triangulation techniques. The pseudorange measure deduced from the signal propagation time is expressed as follow :

ρ s t = d s t -c(δt u + δt s ) + I s t + T s t + m s t + w t (1) where
• ρ s t is the pseudorange between the satellite s and the receiver at time t; • d s t = p tp s t is the true satellite-receiver distance; p t = [x t , y t , z t ] T is the three position coordinates of the vehicle in the reference system of coordinates East, North, Up (ENU) 1 and p s t = [x s t , y s t , z s t ] T is a vector containing the position coordinates of the satellite s;

• δt u is the receiver clock offset with respect to the GPS reference time 2 and c = 3 × 10 8 m/s is the celerity (the speed of light);

• δt s is the satellite clock offset;

• I s t and T s t are respectively the ionospheric and the tropospheric errors; • w t is the receiver noise considered as a white Gaussian random variable with variance σ 2 is general equal to unity;

• m s t is the error due to the signal reflections in case of multipath propagations.

The atmospheric propagation errors and the satellites clock bias can be corrected by correction models [START_REF] Kaplan | understanding GPS: Priciples and Applications[END_REF].

Consequently, after corrections, (1) can be rewritten as (2):

ρ s t = d s t -cδt u + m s t + w t (2)
It follows that the pseudorange error can be defined in function of two different error sources as :

ǫ s t = m s t + w t
where ǫ s t is the sum of the error due to the signal reflections between the satellite s and the receiver and the receiver noise. According to the reception conditions, ǫ s t can switch between different observation models as following :

LOS : m s t = 0, ǫ s t ∼ N (0, σ 2 ) non-LOS : m s t = 0, ǫ s t ≁ N (0, σ 2 ) (3)
In the case when a signal is received in LOS, the pseudorange error distribution is considered white-Gaussian.

In the other case, the pseudorange error distribution is unknown. In constrained environments, the density form can change abruptly or evolve slowly during a long period of time according to the obstacle nature. Moreover, several moving obstacles (vehicles, pedestrians, etc) can induce random errors. Therefore, to estimate accurately the pseudorange error density, a flexible model was proposed in this paper. It should be noted that our aim is to estimate at each time t the position p t and the observation noise density m s t from the set of collected measurements {ρ s1 t , ρ s2 t , . . . , ρ sl t } where l is the number of visible satellites at time t.

III. BAYESIAN MODELING

In this section we are interested in designing a dynamic model for the parameters to be estimated. Our ultimate objective is to estimate accurately the state of a moving vehicle. We design by state the vector containing vehicle positions p t = [x t , y t , z t ] T and velocities ṗ = [x t , y t , z t ] T . Since we are interested in estimating jointly the hidden state x t = [p t , ṗt ] and the pseudorange noise density, we will consider an augmented hidden state vector containing both hidden state parameters and observation noise parameters. The observation error ǫ t is distributed according to 1 This Cartesian coordinate system is far more intuitive and practical than ECEF or Geodetic coordinates [START_REF] Kaplan | understanding GPS: Priciples and Applications[END_REF]. The local ENU coordinates are formed from a plane tangent to the Earth's surface fixed to a specific location and hence it is sometimes known as a "Local Tangent"

or "local geodetic" plane. By convention the east axis is labeled x, the north y and the up z.

2 GPS has a reference time called GPS time maintained on the earth at the US naval observatory. The satellite clocks, though very accurate, are different pieces of equipment and so are not, in general, exactly synchronized with each other. Thus, they have different offsets with respect to GPS time.

an unknown probability density function F whose parameters denoted θ t is be estimated. Let E t = [x t , θ t ] stands for the vector containing all the parameters to be estimated.

In the following, a state-space approach providing a general framework for describing the dynamic state estimation problem at hand is provided. Here, we will use the state-space formulation to introduce the filtering algorithms.

Statistical Bayesian filtering aims at computing the posterior probability density function of a state vector x t ∈ X from sequentially obtained sensor measurements (pseudoranges) ρ t ∈ Z. In our context, it relies on the following model :

x t = f(x t-1 ) + v t (4)
ρ t = h(x t , ǫ t ) (5)
ǫ t ∼ F (.) (6)
where f is the state transition function, h is the observation function, they are assumed to be nonlinear, and v t and ǫ t are respectively the state and the measurement noise vectors. Note that in this work we will be interested in estimating ǫ t and not v t . We will assume that v t is distributed according to a known density function with fixed parameters.

We are considering a real-world data analysis tasks which involves estimating unknown quantities from some given observations. In most of such applications, prior knowledge about the phenomenon being modeled is available.

This knowledge allows us to formulate Bayesian models, that is prior distributions for the unknown quantities and likelihood functions relating these quantities to the observations. Within this setting, all inference on the unknown quantities is based on the posterior distribution obtained from the Bayes' theorem. Often, the observations arrive sequentially in time and one is interested in performing inference on-line. It is therefore necessary to update the posterior distribution as data become available. It is worth noting that computational simplicity in the form of not having to store all the data might be an additional motivating factor for sequential methods.

In this paper, we will focus on SMC integration methods [START_REF] Doucet | Sequential Monte Carlo methods in practice[END_REF] which are a set of simulation-based techniques that provide a convenient approach to computing the posterior distributions in our nonlinear and non-Gaussian case and lead most of the time to accurate filtering results. In the context of positioning with GPS, it is typically practical to process data on-line as it arrives from the visible satellites at each time t, both from the point of view of storage costs as well as for rapid adaptation to changing signal characteristics. In this paper the objective is to study a dynamic nonlinear model with additive observation errors. The inferential problem is to estimate sequentially the hidden states driving the observations along with the starting value. At the same time, our focus is to model the additive observation errors using a highly flexible family of density functions based on an infinite mixture model. In the following, in order to introduce the statistical Bayesian filtering, we first review some theoretical aspects about nonlinear dynamic state-space modeling. Then, we will provide some typical ways dealing with density estimation of errors via mixture based models.

A. Nonlinear dynamic state-space model and statistical Bayesian filtering

The unobserved signals {x t ; t ∈ N} are modeled as a Markov Process of initial distribution p(x 0 ) and transition equation p(x t |x t-1 ) provided by the motion model ( 4). The observations {ρ t ; t ∈ N ⋆ }, are assumed to be conditionally independent given the process {x t ; t ∈ N} and the marginal distribution p(ρ t |x t ) deduced from (5). To sum up, the model is described by

       p(x 0 ) p(x t |x t-1 ), t ≥ 1 p(ρ t |x t ), t ≥ 1
We denote by x 0:t {x 0 , . . . , x t } and ρ 1:t {ρ 1 , . . . , ρ t }, respectively, the states and the observations up to time t. Our aim is to estimate recursively in time the posterior distribution p(x 0:t |ρ 1:t ) and its associated features (including the marginal distribution p(x t |ρ 1:t ), known as the filtering distribution). The posterior distribution can be calculated in two steps theoretically: prediction and update. In the prediction step, we integrate the state distribution from the previous state using the system model. The update operation modifies the prediction distribution making use of the latest observation. Given the initial distribution p(x 0 ), transition distribution p(x t |x t-1 ), and the likelihood distribution p(ρ t |x t ), the objective of the filtering is to estimate the state at time t given the observations up to time t. From Bayesian perspective, when an observation ρ t becomes available at time t, we can obtain this posterior distribution p(x 0:t |ρ 1:t ) via the Bayes rule. The expression of this distribution requires the evaluation of complex high-dimensional integrals. Since we are in a nonlinear case, the multi-dimensional integrals are intractable and some approximate methods must be used. More details about the used filtering approach will be given in section V.

B. Bayesian analysis of mixtures for density estimation

In this section, a Bayesian nonparametric approach based on a mixture model for the density estimation in a nonlinear dynamic model will be introduced. Later, we will study the inference procedure when the nonparametric component is applied to the additive observation errors {ǫ t ; t ≥ 1}. More precisely, we will show that by using a Dirichlet process prior [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF], [START_REF] Escobar | Bayesian density estimation and inference using mixtures[END_REF], [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF] for the observation noise density F , we can derive a Bayesian nonparametric mixture model. Before we discuss the Dirichlet process prior and the nonparametric Bayes' analysis, let us first consider some typical ways to deal with a mixture model.

1) General

Framework: Suppose we have observed data ρ 1 , . . . , ρ t , such that given (µ 1 , Σ 1 ), . . . , (µ t , Σ t ), the ρ i 's are independent, and the ρ i has a normal distribution with mean µ i and variance Σ i . Let θ i = (µ i , Σ i ). Suppose a random variable ρ comes from any one from k components, each of them having a different probability density function. If ρ comes from i-th component with probability π i and pdf f (ρ|θ i ), then the pdf of an observation ρ can be represented in the mixture form

f (ρ|φ) = k i=1 π i f (ρ|θ i ), (7)
where φ = {k, π 1 , . . . , π k , θ 1 , . . . , θ k }, k is the number of components, the π i 's are the mixture weights which must be non-negative and sum to one and the θ i 's denote the vector of all unknown parameters associated with ith component. In some applications it is desired to categorize or classify each observation into one of the k populations. under a Bayesian approach this can be done by computing a classification probability, the probability that an observation comes a certain population. Alternatively, one may be interested in the elements of φ, which include both the mixing weights π 1 , . . . , π k , and the density parameters θ 1 , . . . , θ k . The functional form of each distribution is usually assumed to be known. In some cases the number of components k may not be known, and is estimated from the data, as in cluster analysis [START_REF] Fowlkes | A method for comparing two hierarchical clusterings[END_REF]. Two kinds of Bayesian analysis of mixtures, parametric and nonparametric, have been developed. Parametric Bayesian analysis incorporates prior information about interesting parameters in the inference. Sources of prior information include beliefs about parameters before any data are observed. On the other hand, Bayesian nonparametric mixture modeling has been developed recently. In this paper, our concern is about the nonparametric mixture models.

2) Nonparametric Bayesian Analysis of Mixture Distributions: Suppose data ρ i (i = 1, . . . , t), are distributed as f (ρ i |θ i ), independently over i. The mixture structure simply imposes the constraint that, for some positive integer k, there exist k distinct numbers θ * = (θ * 1 , . . . , θ * k ) such that, for each i = 1, . . . , t, θ i = θ * j for some j = 1, . . . , k. One way of generating such a mixture is to use a Dirichlet process for the prior distribution of θ = (θ 1 , . . . , θ k ). An important result of the Dirichlet process is that the probability of θ i = θ j for i = j is positive. Such Bayesian analyses have been of limited use due to the computational difficulties. However, these difficulties can now be overcome using Gibbs sampling techniques [START_REF] Escobar | Bayesian density estimation and inference using mixtures[END_REF]. This permits the extension of standard Bayesian parametric analyses to Bayesian nonparametric analyses, by linking the Gibbs sampling computations for Dirichlet process priors to those for the standard models.

IV. BAYESIAN NONPARAMETRIC DENSITY ESTIMATION

In [START_REF] Viandier | Gnss pseudorange error density tracking using dirichlet process mixture[END_REF] it was shown that an important characteristic of the observation error is the non-stationarity. In the context of dynamic models, the assumption of stationarity is wrong. In fact, when reflections occur; and by considering a long period of observation time, the pseudorange error distribution can be easily approximated by a finite gaussian mixture model (Fig. 1). However, by considering a short period of time, pseudorange errors can be modeled by only one gaussian distribution whose parameters (mean and variance) evolve over time. This drawback can be limited by the use of DPM. DPM models generate an infinite Gaussian mixture. In this mixture model the Dirichlet process prior is used to specify latent patterns of heterogeneity, particularly when the distribution of latent parameters is thought to be multimodal. It will shown through this paper that this method is flexible and well suited for on-line applications. In traditional application of mixture models, namely, the estimation of a (univariate) density function, the purpose is to determine the appropriate number of mixture components as well as to estimate their corresponding parameters.

More generally, the problem of density estimation arises when data are considered to be sampled from a certain distribution but this distribution is itself assumed unknown. Such problems (and their multiple variations) arise in nearly all branches of Statistics. Nonparametric Bayesian methods are specially well suited for density estimation.

References can be found, among others in [START_REF] Fox | Hierarchical dirichlet processes for tracking manoeuvering targets[END_REF], [START_REF] Teh | Hierarchical Bayesian nonparametric models with applications[END_REF], [START_REF] Griffiths | Infinite latent feature models and the indian buffet process[END_REF], [START_REF] Muller | Nonparametric Bayesian data analysis[END_REF]. The focus of our work will be on mixtures of normals. As it will be explained in this paper, the main motivation behind this choice is the application context.

A. Nonparametric Density estimation

In mixture models the Dirichlet Process is used to specify latent patterns of heterogeneity, particularly when the distribution of latent parameters is thought to be clustered. However, neither the number nor the location of clusters needs to be specified a priori. Therefore, DP mixture models are considered very attractive because DP mixtures allow for uncertainty in the choice of parametric forms and in the number of mixing components. The DP mixture model has found widespread application in recent statistical research [START_REF] Ishwaran | Approximate Dirichlet process computing in finite normal mixtures: smoothing and prior information[END_REF], [START_REF] Brown | A Bayesian semiparametric joint hierarchical model for longitudinal and survival data[END_REF], [START_REF] Iorio | An ANOVA model for dependent random measures[END_REF], [START_REF] Muller | A method for combining inference across related nonparametric Bayesian models[END_REF]. Analyses based on DP mixtures are sensitive to the choice of the Dirichlet process parameters; the need to treat these parameters carefully is discussed in [START_REF] Escobar | Bayesian density estimation and inference using mixtures[END_REF].

Let consider a pdf F and ρ 1 , . . . ρ t a set of vectors statistically distributed according to F , ρ k ∼ F (.) for 1 ≤ k ≤ t. We will consider a nonparametric model allowing to estimate F as follows:

F (ρ) = θ f (ρ|θ)dG(θ) (8)
where θ ∈ Θ is called the latent variable (also called cluster), f (.|θ) is the mixed pdf and G is the mixing distribution.

Under a Bayesian framework, it is assumed that G is a Random Probability Measure (RPM) distributed according to a prior distribution. In this paper, we have selected as a RPM the Dirichlet Process (DP) prior.

B. Dirichlet Processes

In [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF] the Dirichlet Process (DP) has been introduced as a probability measure on probability measures. The DP is parameterized by a base distribution G 0 on a (measurable) space and a positive scaling parameter α. Let's suppose now that we draw a random measure G from a DP, and independently draw random variables θ k from G :

θ k |G ∼ G (9)
where G ∼ DP (G 0 , α). Marginalizing out the random measure G, the joint distribution of {θ 1 , . . . , θ t } follows a Polya urn scheme [START_REF] Blackwell | Ferguson distributions via Polya urn schemes[END_REF]. In another representation of the DP called the stick-breaking scheme, G is introduced explicitly as an infinite sum of atomic measures [START_REF] Stthuraman | A constructive definition of Dirichlet priors[END_REF]. The realizations of a DP are expressed as follow :

G = ∞ k=1 π k δ θk (10) with θ k ∼ G 0 , π k = β k k-1 j=1 (1 -β j ) and β k ∼ B(1, α)
where B denotes the Beta distribution and δ θk denotes the Dirac delta measure located in θ k . The underlying random measure G is then discrete with probability one. Using [START_REF] Mazor | Interacting multiple model methods in target tracking: a survey[END_REF], it comes that the following flexible prior model is adopted for the unknown distribution F

F (y) = ∞ k=1 π k f (y|θ k ) (11)
More details about other important properties of Dirichlet processes can be found in [START_REF] Antoniak | Mixtures of dirichet processes with applications to bayesian nonparametric problems[END_REF], [START_REF] Blackwell | Ferguson distributions via Polya urn schemes[END_REF], [START_REF] Stthuraman | A constructive definition of Dirichlet priors[END_REF], [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF].

C. Dirichlet Process Mixture (DPM)

DP mixture models form a very rich class of Bayesian nonparametric models.They emerge by employing a DP prior for the mixing distribution in a mixture of a parametric family of distributions. DP mixtures have dominated the Bayesian nonparametric literature after the machinery for their fitting, using Markov chain Monte Carlo (MCMC) methods, was developed following the work in [START_REF] Escobar | Bayesian density estimation and inference using mixtures[END_REF].

DPM models [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF] have become increasingly popular for modeling when conventional parametric models would impose unreasonably stiff constraints on the distributional assumptions. Recall that in the finite mixture model, each data point is drawn from one of fixed distributions. To allow the number of mixture components to grow with the data, we move to the general mixture model setting. This is best understood with the hierarchical graphical model in [START_REF] Escobar | Bayesian density estimation and inference using mixtures[END_REF]. The core of the DPM model can basically be thought of as a simple Bayes model given by the likelihood

y k ∼ f (y k |θ k ) and prior θ k ∼ G(θ k ),
with added uncertainty about the prior distribution G :

y k |θ k ∼ f (y k |θ k ), k = 1, . . . , t θ k |G ∼ G G ∼ DP (G 0 , α)
Through this hierarchy, a marginalization step is often used in actual implementations. [START_REF] Blackwell | Ferguson distributions via Polya urn schemes[END_REF] proved that by integrating over G the joint distribution of θ = (θ 1 , • • • , θ t ) may be factored into a product of successive conditional distributions of the following form :

θ 1 |G 0 ∼ G 0 (12)
θ k |θ 1 , • • • , θ k-1 , α, G 0 ∼ 1 k -1 + α k-1 j=1 δ(θ j ) + α k -1 + α G 0 , for 2 ≤ k ≤ t
This factorization implies that θ is randomly partitioned into classes of distinct values such that the elements of θ within a class share the same value. Therefore, a random draw from the Dirichlet process DP (α, G 0 ) may be computed as follows :

θ t+1 |θ 1:t , G = θ * c , with probability n t,c /(α + t) G 0 , with probability α/(α + t) (13)
where

θ * c (c = 1, 2, • • • , k c ) denotes the c th of k c distinct values in θ 1:t = (θ 1 , • • • , θ t
) and n t,c denotes the number of elements in θ 1:t that equals θ * c . In other words, a random draw from G either equals one of the previous draws or is drawn independently from the base probability measure G 0 . This is due to the discreteness of the random measure G.

D. Estimation of DP Hyperparameters

In a hierarchical model, a hyperprior can be placed on the parameters of the DP by simply putting prior distributions on α and the parameters of G 0 . The data is then used to calculate the posterior distribution. We can say that the DP adds another layer in the hierarchical model. The importance of treating the Dirichlet Process hyperparameters carefully is discussed in [START_REF] Escobar | Bayesian density estimation and inference using mixtures[END_REF]. In applications, some authors have chosen a fixed G 0 , such as a normal distribution with large variance relative to the data [START_REF] Ishwaran | Approximate Dirichlet process computing in finite normal mixtures: smoothing and prior information[END_REF] or a Beta(a, b), where a and b are stipulated without further discussion [START_REF] Liu | Nonparametric hierarchical Bayes via sequential imputations[END_REF]. Other authors have used a parametric family of priors for G 0 , such as uniform-inverse Wishart [START_REF] Maceachern | Estimating mixture of Dirichlet process models[END_REF] or normal-inverse gamma [START_REF] Maceachern | Estimating mixture of Dirichlet process models[END_REF], [START_REF] Ghahramani | The infinite gaussian mixture model[END_REF], with fixed hyperprior distributions on the parameters of the prior family. The hyperprior distributions are usually chosen by sensitivity analysis, or else they are chosen to be diffuse with respect to the observed data.

In this paper, Concerning the centering distribution G 0 , priors on G 0 hyperparameters will be specified using a data-adaptive way without applying a pure sampling approach. Most of the time, in real-world application contexts, the observed data can help in defining a prior for some of the hyperparameters in the considered model. Computing every unknown hyperparameter with pure sampling is not the most effective approach in all cases since this can be practically difficult and time consuming.

Since α is an important parameter is this problem, one may want to put a prior on α and then use the data to calculate a posterior distribution. In understanding what values of α to use, the relationship between α and the expected number of clusters in the data might be considered. In [START_REF] Escobar | Bayesian density estimation and inference using mixtures[END_REF] was developed a clever Gibbs sampler wherein draws from the conditional posterior distribution of α are computed by drawing successive samples from relatively familiar distributions (Beta and Gamma).

E. Prior and posterior distributions for α

In [START_REF] West | Hyperparameter estimation in dirichlet process mixture models[END_REF], the author shows how, with respect to a flexible class of prior distributions for parameter α, the posterior may be represented in a simple conditional form that is easily simulated. As a result, for a better fitting of DPM models, inference about this parameter may be developed with existing Gibbs sampling algorithms [START_REF] Escobar | Bayesian density estimation and inference using mixtures[END_REF].

We have mentioned in this paper that a key feature of the Dirichlet is its discreteness, which in our context implies that the pairs (µ k , Σ k ), (k = 1, . . . , t), concentrate on a set of some l ≤ t distinct pairs. To sample the precision parameter α, we first determine the prior distribution for l, the number of normal components in the mixture (clusters). At each stage of the simulation analysis, a specific value of l is simulated from the posterior for l (together with sampled values of the means and variances of the normal components) which also depends critically on this hyperparameter α. In [START_REF] Escobar | Bayesian density estimation and inference using mixtures[END_REF], [START_REF] West | Hyperparameter estimation in dirichlet process mixture models[END_REF] it was shown how, based on a specific but flexible family of prior distributions for α, the parameter vector θ may be augmented to allow for simulation of the full joint posterior now including α.

As in [START_REF] West | Hyperparameter estimation in dirichlet process mixture models[END_REF], a Gamma prior for α will be used. Suppose α ∼ G(a, b), a gamma prior with shape a > 0 and scale b > 0. In this case, p(l|α, t) may be expressed as a mixture of two Gamma posteriors, and the conditional distribution of the mixing parameter given α and l is a simple Beta.

(α|x, l) ∼ π x G(a + l, b -log(x)) + (1 -π x )G(a + l -1, b -log(x)) (14)
with weights π x defined by

π x (1 -π x ) = (a + l -1) t(b -log(x)) ( 15
)
where 0 < x < 1, l > 1 and (x|α, l) ∼ β(α + 1, t). More details about this sampling can be found in [START_REF] West | Hyperparameter estimation in dirichlet process mixture models[END_REF]. It is important to mention that α should be sampled at each Gibbs iteration stage in the simulation process that allow to generate the θ k according to [START_REF] Dorazio | On selecting a prior for the precision parameter of Dirichlet process mixture models[END_REF]. The current sampled values of l and α permit to compute a new value of α.

In [START_REF] Antoniak | Mixtures of dirichet processes with applications to bayesian nonparametric problems[END_REF] the expected number of components sampled from a Dirichlet process is given by

E[l|α, t] = t l=1 l p(l|α, t) (16) = α t l=1 1 α + l -1 ≈ α ln( t + α α )
Thus, although l → ∞ with probability 1 as t → ∞, the number of components increases approximately logarithmically with the number of observations.

F. The G 0 parameters DP mixtures have received a considerable interest in the Bayesian nonparametric literature. Despite this activity, little has been written on choosing a centering distribution G 0 for these models. In this work, the latent sources of heterogeneity to be specified through the G 0 (µ, Σ) parameters are considered using a data-adaptive way. G 0 is a normal distribution with unknown mean µ and covariance Σ. The couple (µ, Σ) can be chosen to be distributed according to a normal-scaled inverse Wishart distribution, denoted as :

(µ, Σ) ∼ N W -1 (λ, ν, Σ 0 , dof ) (17)
such that :

µ|Σ, λ, ν ∼ N (λ, Σ/ν) Σ|Σ 0 , dof ∼ W -1 (Σ 0 , dof )
where λ is the mean of the Normal distribution, ν is a scale parameter, Σ 0 and dof are respectively the inverse scale matrix and the degree of freedom for the inverse Wishart distribution W -1 (see [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF] for more details about these distributions).

The same base distribution G 0 should not be applied in all navigation situations (LOS, non-LOS and blocked).

Here, we propose to adapt the priors on the G 0 parameters at each time t. More precisely, we propose to use an additional parameter r s t defined at each time t for each satellite s. This parameter adds one more step in the DPM hierarchical model as shown in Fig. 2. The parameter r s t is an indicator of the propagation state for each satellite s at each time t. It will be used to update the parameters of G 0 . In practice, when a LOS reception state is detected the hyperparameters of the base distribution G 0 will be sampled from the normal-scaled inverse Wishart distribution whose parameters λ, ν, Σ 0 and dof should be chosen carefully. However, when a non-LOS reception state is detected, it will be interesting to adapt the distribution of the hyperparameters as following :

if r s t = 0 (LOS) then, (µ, Σ) ∼ N W -1 (λ 1 , ν 1 , Σ 01 , dof 1 ) if r s t = 1 (non-LOS) then, (µ, Σ) ∼ N W -1 (λ 2 , ν 2 , Σ 02 , dof 2 ), (18)
with Σ 02 > Σ 01 , dof 2 > dof 1 , k 2 < k 1 and λ 2 = λ 1 + res s t . The additive term res s t is the estimation of the pseudorange error computed at each time t as following : where ρ s t is the corrected pseudorange measurement and p t|t-1 = x t|t-1 , y t|t-1 , z t|t-1 T is the predicted position and b t|t-1 = cδt u is the predicted receiver clock bias computed at time t.

res s t = ρ s t -p t|t-1 -p s t -b t|t-1 (19)

V. RAO-BLACKWELLIZED PARTICLE FILTERING FOR DYNAMIC BAYESIAN MODELS

A. Sequential Monte Carlo Filtering

Unscented Kalman filters (UKF) and Extended Kalman filters (EKF) which are widely used as optimal solution in case of white Gaussian systems are not suitable for treating the models which exhibit non centered and non Gaussian distributions. Since, as we know, real data can be very complex. Except in a few simple cases, there is no closed-form solution to this problem. It is therefore necessary to adopt numerical techniques in order to compute reasonable approximations. Sequential Monte Carlo methods (particle filters) are powerful tools that allow us to accomplish this goal. This filtering method is able to cope with noises of any distribution and provide a convenient and attractive approach to computing the posterior distributions.

The particle filter approximates the posterior probability density function via the discrete weighted sum:

p(x t |ρ 1:t ) ≈ N i=1 w (i) t δ(x t -x (i) t ) (20)
The individual weights w i , are computed by applying the principle of importance sampling. Since it is difficult or impossible to directly sample p(x t |ρ 1:t ), particles are simulated according to a proposal distribution :

x (i) t ∼ q(x t |x (i) 1:t-1 , ρ 1:t ), for i = 1, . . . , N
Then, to correct for the discrepancy between the proposal and the target distribution, importance weights will be assigned to the generated particles :

w t = w(i) t / N k=1 w(k) t where wt = w (i) t-1 p(ρ t |x (i) 1:t-1 , ρ 1:t-1 )p(x (i) t |x (i) 1:t-1 ) q(x (i) t |x (i) 1:t-1 , ρ 1:t ) (21)
According to [START_REF] Casella | Rao-Blackwellisation of sampling schemes[END_REF], the set of particles is updated and reweighted using a recursive version of importance sampling.

The most likely particles yield high importance weights. An additional resampling procedure is used for replacing particles with large weights and removing those with small weights. The variance introduced by the resampling procedure can be reduced by proper choice of the resampling method [START_REF] Doucet | Sequential Monte Carlo methods in practice[END_REF].

B. Improving the Particle Filtering by Rao-Blackwellisation

A method to increase the efficiency of sampling techniques is to reduce the size of the state space by marginalizing out if possible some of the variables analytically; this is called Rao-Blackwellisation [START_REF] Casella | Rao-Blackwellisation of sampling schemes[END_REF]. In the following, a Rao-Blackwellised particle filtering (RBPF) [START_REF] Doucet | On sequential monte carlo sampling methods for bayesian filtering[END_REF] will be introduced in order to improve filtering.

As showed previously, we can compute the posterior distribution on-line by applying Bayes' rule sequentially.

However, in general, the required integrals cannot be computed in closed form. Particle filtering therefore approximates the posterior using sequential importance sampling. Now, if we partition the state-space into two sub-spaces, drawn by θ t and x t (for all t ≥ 1). Then, by the chain rule of probability, we can write :

p(x 0:t , θ 1:t |ρ 1:t ) = p(x 0:t |θ 1:t , ρ 1:t )p(θ 1:t |ρ 1:t ) (22)
If we can update p(x 0:t |θ 1:t , ρ 1:t ) analytically and efficiently, then we only need to sample p(θ 1:t |ρ 1:t ) using the particle filter. Since we are now sampling in a smaller space, in general we will need far fewer particles to reach the same accuracy as standard PF. This is the key idea behind RBPF.

The following is a generalization of the Rao-Blackwellized particle filter to DPM as presented in [START_REF] Doucet | Sequential Monte Carlo methods in practice[END_REF] and applied in [START_REF] Caron | Bayesian inference for dynamic models with dirichlet process mixtures[END_REF]. At time t, we use the following empirical distribution to approximate p(x t |ρ 1:t ) through a set of N particles θ

1:t , . . . , θ (N )

1:t , P N (x t |ρ 1:t ) = N i=1 w (i) t p(x t |θ (i) 1:t , ρ 1:t ) (23)
with

p(x t |θ (i) 1:t , ρ 1:t ) = N (x t ; x t|t (θ (i) 1:t ), Σ t|t (θ (i) 1:t )) (24)
By using an Extended or an Unscented KF, for each particle i we can compute recursively the terms x t|t (θ

(i) 1:t ) and Σ t|t (θ (i) 1:t ).
For more details about how to build the algorithm refer to [START_REF] Caron | Bayesian inference for dynamic models with dirichlet process mixtures[END_REF], [START_REF] Caron | Particle filtering for multisensor data fusion with switching observation models. Application to land vehicle positioning[END_REF]. In RBPF, each particle maintains not just a sample from p(θ 1:t |ρ 1:t ) but also a parametric representation of the distribution p(x t |θ (i) 1:t , ρ 1:t ) (the parametric representation in our case are a mean vector and a covariance matrix). For our application, the θ t samples are updated as in standard PF, and then the x t distributions are updated using an EKF, conditional on θ t . The overall algorithm is shown in Algorithm 1.

Following what we presented above, the optimal importance distribution is q(θ t |θ (i)

1:t-1 , ρ 1:t ) = p(θ t |θ (i) 1:t-1 , ρ 1:t ).
Without deriving approximations of this importance distribution, no efficient importance distributions can be used directly since the associated importance weights will be computationally intractable. In this paper, for the sake of simplicity, the importance distribution is chosen to be equal to the Polya urn distribution p(θ t |θ (i) 1:t-1 ) as in [START_REF] Dorazio | On selecting a prior for the precision parameter of Dirichlet process mixture models[END_REF]. This approximation yields to the simplification of the weight expression [START_REF] Casella | Rao-Blackwellisation of sampling schemes[END_REF].

Finally, the minimum mean squared error (MMSE) estimate and posterior covariance matrix of x t are given by

x t|t = N i=1 w (i) t x (i) t|t (25) Σ t|t = N i=1 w (i) t [Σ (i) t|t + ( x (i) t|t -x t|t )( x (i) t|t -x t|t ) T ] (26)
Algorithm 1: Computing the posterior distribution p(x 1:t , θ 1:t |ρ 1:t ) using a Rao-Blackwellized Particle Filter

At time t = 0,

Step 0: initialization

• For i = 1, ..., N , sample (x (i) 0 , Σ (i) 0|0 ) ∼ p 0 (x 0 ) and set t = 1 • ω (i) 0 = 1 N At time t ≥ 1,
Step 1: Importance Sampling step

• For i = 1, . . . , N -Sample θ (i) t ∼ p(θ t |θ (i) 1:t-1 ) -Evaluate the importance weights ω (i) t ∝ ω (i) t-1 p(z t |θ (i) 1:t-1 , θ (i) t , ρ 1:t-1 )
Step 2: Selection step

• Multiply/Discard particles ( θ (i) t ; i = 1, . . . , N ) with respect to high/low normalized importance weights ω (i) t to obtain N equally weighted particles.

Step 3: Apply an Extended Kalman Filter to compute p(x t | θ

(i) t , z 1:t ) • Compute (x t|t-1 (θ (i) 1:t ), Σ t|t-1 (θ (i) 1:t ), x t|t (θ (i) 1:t ), Σ t|t (θ (i) 1:t ) =EKF(x t-1|t-1 (θ (i) 1:t-1 ), Σ t-1|t-1 (θ (i) 1:t-1 ), x t|t (θ (i) 1:t ), Σ t-1|t-1 (θ (i) 1:t-1 , θ (i)
t , ρ t ))

VI. APPLICATIONS TO GNSS POSITIONING IN URBAN CANYON

A. Error Density Modeling using DPM

For density estimation, the heterogeneity in the distribution of the latent parameters is a problematic issue in Bayesian nonparametric modeling. We know that in real applications, the sources of this heterogeneity cannot be completely known. We will prove through some experiments that the Dirichlet process prior may provide a satisfactory model in case where the sources of heterogeneity in the observation errors density are unobserved.

The DP prior is robust to errors in model specification and allows heterogeneity in the patterns distributions to be specified in a data-adaptive way.

For example, the predictive distribution of densities based on our Dirichlet process prior can be multimodal (NLOS case) owing to its natural potential for identifying clusters. In contrast, standard parametric prior distributions of densities, such as the normal distribution, cannot capture multimodality because of their shape restrictions. But, if the actual distribution of density can be approximated by a unimodal distribution (LOS case), the Dirichlet process prior also adapts to this situation. In fact, the parametric model based on finite Gaussian mixtures is a limiting case of our modeling scheme when applying Dirichlet process prior with α → ∞. Therefore, we can say that the Dirichlet process prior is a natural extension of fully parametric models.

In our application of on-line pseudorange error density estimation, the latent variables θ i are the mean and the variance of each Gaussian distribution included in the infinite mixture, i.e. θ i = (µ i , σ i ), for i ≥ 1. To compute a vehicle position at each step of the filtering process, N particles θ 1 i , • • • , θ N i are computed. The problem with GNSS applications is that signal propagation can be considered in different reception modes. In LOS reception, the pseudorange error noise follows a zero-mean Gaussian distribution with a standard deviation σ inferior or equal to 1 (actually it depends on the used receiver). In NLOS reception, the pseudorange error can be higher than 100 m. In the proposed DPM model, the α will be estimated by sampling as detailed in section IV-D. The same base distribution G 0 will not be applied in all navigation situations. Here, we propose to adapt α and G 0 each time using the proposed approaches in section IV-D.

σ t,s ∼ W -1 (Σ 1 , dof 1 ) * If NLOS reception, set the following parameters of G 0 µ t,s ∼ N (µ 0 + res s t , σ t,s /k 2 ) σ t,s ∼ W -1 (Σ 2 , dof 2 )
Algorithm 3: Sampling the α For each iteration i of the Gibbs sampler (i = 1, . . . , N ), do

• Set the parameter α to some predetermined value

• Consider n t the number of measurements provided by the satellite s at time t

• Compute the number of components in the mixture (clusters) at time t denoted k t , (k t ≤ n t )

• Set a prior for α : α ∼ G(a, b) (a, b > 0)

• Compute a conditional distribution of the mixing parameter denoted x using the Beta function as following (x|α, k t ) ∼ β(α + 1, n t )

• Use x to compute weights π x as in [START_REF] Muller | Nonparametric Bayesian data analysis[END_REF] • p(α|k t , n t ) is a mixture of two gamma posteriors computed as in ( 14)

B. Experimental set-up description

The Ergospace software simulates electromagnetic signal propagation in 3D realistic constrained environments.

This software is developed to be a support for the development of GNSS applications. The deterministic method of Ray Tracing determines possible paths of received rays (example given in Fig. 3). Propagation models characterize interactions between the signal and the studied environment. Receiver behavior (mobile or other) and performance are also taken into account in calculations.

The software provides output data files in Matlab format that we process in order to compute the position. In this section, simulations have been performed with a 3D scene of the town of Rouen, in France. The simulated route corresponds to an existing bus line. In these simulations, the mobile ran at a constant speed of 30 km/h, during 10 minutes and is covering a distance of 8236 meters. GNSS signals from both GPS and Galileo constellations can be received with a maximum of 3 reflections. In most urban transport applications, localization should be provided with a good accuracy. In our experiments, we have fixed the accuracy threshold to 3 meters [START_REF] Viandier | Positioning urban buses: GNSS performances[END_REF].

C. Results and discussion

In this part we will present some results that pertain to the choosing of the DPM model parameters. Moreover, we designed experiments involving some comparison with some commonly used estimation and filtering methods to show the robustness, flexibility and merits of our approach. Note that in the following, the DPM simulations are based on the implementation of the DP mixture model described in the previous subsection, applied to the case of conjugate models with normal structure. The precision parameter α may be calibrated automatically following the procedure proposed in Algorithm 3. The parameters of the centering distribution G 0 are simply adapted to the data characteristics as shown is Algorithm 2.

In order to highlight the effectiveness of the proposed algorithm based on RBPF, we have conducted preliminary experiments using EKF, PF and RBPF algorithms. Using these methods, localization performances were compared along the same trajectory during the same period of time. In the EKF algorithm, pseudorange error distributions are assumed to be centered Gaussian. In the PF algorithm, we considered a finite Gaussian mixture. About 500 particles are used in the above algorithms. Accuracy has been calculated with reference to an exact position given by the simulator (Fig. 4).

The use of RBPF algorithm is fully justified by the results presented below (Fig. 5 and Fig. 6), as it yields consistently lower error rates and a high localization accuracy. It appears that using DPM models can improve the accuracy and the availability of the localization even in very constrained environments.

Fig. 6 shows the evolution over time of the positioning error. The proposed algorithm has proven to be more accurate than the other algorithms. Besides, another key point permitting to conclude that those methods are more suitable than the others is the algorithm stability when abrupt propagation changes occur during mobile dynamic evolution (Fig. 5). The GM and the DPM are used in a sequential Monte Carlo filtering approach. Table I shows the localization performances for the two models in terms of mean error and availability of the navigation solution. An indicator of the availability of the navigation solution is a threshold over the error values which is set equal to 3 meters.

The parameters of a DP include the precision measure α and the base probability measure G 0 . In problems where α is unknown, estimation results can be extremely sensitive to the choice of prior assumed for α [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF]. In most applications of DPM models, the precision parameter α is not known and must be well chosen. In this paper, some empirical tests were conducted to show the influence of α values on the estimation performance.

In fact, a random draw from DP (α, G 0 ) either equals one of the previous draws or is drawn independently from the base probability measure G 0 . The parameter α obviously plays an important role in the distribution of θ. In equation [START_REF] Dorazio | On selecting a prior for the precision parameter of Dirichlet process mixture models[END_REF], note that the probability that θ i differs from all previously drawn parameter values is proportional to α; therefore, higher values of α lead to a higher probability of many unique values (i.e., a higher number of classes relative to the sample size n). Conversely, lower values of α lead to a greater chance of clustered distributions with and their corresponding reception states.

In preliminary experiments, we considered a range of values in order to explore the extend to which our results are sensitive to different choices of α. While the overall shape of the distribution remains the same across values of α, higher values of α appear to produce an increased number of localized features and even indicate additional modes for a subset of the observations. The estimated means and variances of the model parameters are, however, quantitatively very similar. Furthermore, we find little gain from using large values of α. In our DPM model, we used a gamma distribution as a prior for the α. The major advantage of this choice is that with this prior the conditional distributions are easy to sample by applying the approach of section IV-E. Fig. 7 shows that estimating the α leads to a good positionning accuracy compared to the result obtained when we use a fixed α.

The resulting posterior distribution over α based on a gamma prior shown in Fig. 8 is considerably variable. are mostly received in LOS and those of satellite 13 both in LOS and NLOS). It seems to be extremely hard to have such informations on data in a real-world application. Therefore, the parameters a et b of the gamma prior were adjusted carefully. For example, we noticed that small values of a and b lead to nearly similar values of the α probability density, which means a lack of variability in the distribution of θ.

VII. CONCLUSION

Our work deals with the problem of improving localization accuracy in urban canyons when signals delivered by the GNSS system are received most of the time after many reflexions. The errors induced in the pseudorange measurements have an unknown density form. Thus, in this paper we have discussed a Bayesian nonparametric model based on Dirichlet process mixtures that permits robust adaptive inference to be carried out on data with some unknown specifications. The DP prior is extremely flexible, it allows to model a huge variety of distributional forms. In the considered application, its effectiveness was proved, since in our case little is known about priors to be considered for the errors distribution. We illustrated through experiments that a good choice of the DPM hyperparameters leads to a better density estimation. Since the proposed density estimation solution was strongly motivated by the complexity of the considered data, we proved that the proposed approaches outperform standard sequential methods, especially in a nonlinear dynamic context in which joint estimation of states and observation noise distributions is needed.

Many interesting questions remains to be addressed. One of these is whether non-stationarity can be involved in the nonparametric model by considering time-varying Dirichlet process mixtures. At each time step it can be interesting to consider that the unknown distribution follows a DPM model. Some interesting and simple approaches, such as utilized in [START_REF] Caron | Generalized polya urn for time-varying dirichlet process mixtures[END_REF], need to be investigated. A second is to develop more accurate and efficient approximations of the importance distribution q(θ t |θ (i) 1:t-1 , z 1:t ) in the proposed RBPF.

Fig. 1 .

 1 Fig. 1. The pseudorange errors can be approximated by a finite Gaussian mixture model where the number of components in the mixture should be specified in advance and the choice of the analysis window length should be done carefully by considering data evolution aspects.

Fig. 2 .

 2 Fig. 2. DPM hierarchical model adaptation using the indicator r s t of the signal reception state.

Algorithm 2 :

 2 Tuning the hyperparameters of the DP For t = 1, . . . , n • For s = 1, . . . , N t s (N t s : Number of visible satellites at time t) -For each iteration i of the Gibbs sampler (i = 1, . . . , N ), do * Sample α according to Algorithm 3 * If LOS reception, set the following parameters of G 0 µ t,s ∼ N (µ 0 , σ t,s /k 1 )

Fig. 3 .Fig. 4 .

 34 Fig. 3. An example illustrating the vehicle trajectory on a 3D realistic constrained environment (city of Rouen in France) using the Ergospace software simulator.

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig.5. Localization accuracy obtained by applying various filtering algorithms : EKF, PF (with GM for the density estimation) and the proposed approach based on RBPF (with DPM for the density estimation)

Fig. 8 .

 8 Fig. 8. Posterior distribution over the dispersion parameter α when a Gamma prior is used.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Posterior distribution of the number of clusters k

Fig. 11 .

 11 Fig. 11. Locations of the clusters in the DPM model for different satellites (satellites 18, 1 and 13)

Fig. 9

 9 Fig.9illustrates the posterior distribution of k. The posterior clearly indicates that k << n, which suggests that the distribution of the latent variables is clustered. Practically, in most situations a considerable uncertainty exists about the probable value of the number of model clusters k. In our context, nor the prior mean neither the prior variance of the k are known. As shown is[START_REF] Gelfand | Bayesian nonparametric Spatial Modeling with Dirichlet process mixture mixing[END_REF], the number of clusters increases approximately logarithmically with the number of observations. This is illustrated in Fig.10which shows how the prior over the number of components changes as a function of n. Fig.11shows clusters and theirs locations for two satellites (signals from satellite 18

TABLE I COMPARISON

 I OF POSITIONING ERRORS USING DPM AND GM

		Mean error Min error Max error %error<3m
	DPM	2.01	0.0014	8.284	89%
	GM	6.14	0.0019	25.49	51%
	fewer unique values in θ.