Scalable Bayesian Distributional Regression for GNSS Multipath Error Modeling with Uncertainty Quantification

Hayden Dorahy; Chris Hide; Olivier Julien; Ilaria Martini,; Ian Sheret

doi:10.33012/navi.755

Abstract

Global navigation satellite system (GNSS) positioning relies on accurate stochastic models for measurement errors to provide reliable position estimates and bounds. This paper presents a novel approach for modeling errors in GNSS measurements with Bayesian distributional regression, using variational inference to scale model complexity and data set size beyond typical computational limits. This methodology is applied to an automotive data set for GNSS pseudorange (multipath) errors, targeting a Student’s t-distributed model to realistically characterize heavy-tailed data. The distribution is regressed on signal quality indicators from the GNSS receiver to segment different environmental effects on the measurements. Bayesian penalized tensor product splines are used to model nonlinear relationships based on the signal quality indicators. Detailed analyses of goodness-of-fit diagnostics show that the model is able to fit well to the data, and model uncertainty is quantified such that users may be aware of and compensate for inaccuracies in modeling.

Keywords

1 INTRODUCTION

Global navigation satellite system (GNSS) positioning is an inverse problem in which a set of noisy measurements (e.g., satellite pseudoranges from a GNSS receiver) is used to estimate a state (at least the receiver’s position and clock offset). As with any inverse problem, the quality of the estimate depends heavily on the realism of both of the following items:

The functional model, used to describe the relationship between the measurements and state.
The stochastic model, a characterization of the random noise in the system or measurements.

When a GNSS position estimate is given, also provided are a set of position bounds defining a confidence region that reportedly contains the true position (with at least a specified confidence level, e.g., 95%). In actuality, the position bounds returned by typical automotive and consumer GNSS positioning solutions are often unrepresentative of the true system, owing to mismodeling.

In recent years, for safety-critical applications such as autonomous vehicles and drones, requirements for GNSSs have become much more demanding, compared with the needs of aviation in the past. A focus has been placed on not only high-integrity (high confidence level) but also high-precision (small confidence region) position bounds from GNSSs (Wörner et al., 2016), requiring more advanced algorithms and models.

Trust in a system requires trust in its supporting models. In this work, we present a novel method and case study for fitting realistic stochastic models for GNSS measurement errors using distributional regression, which can effectively scale to large amounts of data, with a focus on detailed distributional analysis and quantifying model trust.

1.1 Background

A GNSS receiver uses the time difference between the receiver’s clock at signal reception and the satellite’s clock at signal transmission (scaled by the speed of light) to measure the receiver–satellite distance, referred to as the pseudorange (Teunissen & Montenbruck, 2017). One of the most significant and challenging-to-model error sources in GNSS measurements comes from the effect of local multipath – the interference of direct GNSS signals with those reflected off nearby surfaces. To build a realistic stochastic model for these errors, it is clear that the effect of the environment must be considered. For example, a receiver might observe more severe signal reflections off buildings in a dense urban setting compared with an open field, resulting in larger error magnitudes.

When using GNSS measurements to estimate a receiver’s position (and position bounds), the uncertainty of the errors in the measurements must be specified by a probability distribution. In high-integrity, safety-critical GNSS applications that target reliable position bounds, such as aviation, stochastic models are often developed using “cumulative distribution function (CDF) overbounding” methods (DeCleene, 2000; Rife et al., 2006), which model heavy-tailed measurement errors with probability distributions that overestimate the tail probability density, in order to provide a conservative estimate of the position bounds (Khanafseh et al., 2018). When limited to the use of normal distributions, as often implemented in practice, these models rarely provide an accurate characterization of the stochastic properties of the measurement error. Owing to widespread use of the Kalman filter, standard positioning algorithms are typically restricted to normally distributed error models, which poorly describe heavy-tailed data, although they can benefit in accuracy from model conditioning (Groves & Jiang, 2013).

In contrast, some positioning algorithms use and benefit from the flexibility of non-normal probability distributions, to more closely match the real world and improve estimates. Some notable examples include factor graph navigation solvers (Taylor & Gross, 2024), navigation filters for robust positioning in the presence of outliers (Huang et al., 2016; Roth et al., 2013), and novel GNSS algorithms targeting automotive safety, such as the single-epoch position bound (Bryant et al., 2021), Kalman integrated protection level (Navarro, 2016; Welte et al., 2017), and Student’s t-filters for GNSS integrity (Al Hage et al., 2020, 2023), which use Student’s t-distribution to model heavy-tailed data. We primarily contribute to methods targeted at these applications; however, the methodology presented in this work can also generalize to approximating normally distributed models for standard positioning applications.

1.2 Distributional Regression

In distributional regression, the conditional probability distribution of a response variable¹ is estimated as a function of covariates² (Kneib et al., 2023). By using a distributional regression technique, we can model the distribution of the measurement error based on a set of variables related to the effect of the environment (or other relevant phenomena) on measurements.

Numerous distributional regression approaches have been developed, notably generalized additive models for location, scale, and shape (GAMLSS) (Rigby & Stasinopoulos, 2005). These models use a specified family of distributions to describe the response, can model the distribution parameters as smooth functions of covariates, and are generally fitted using maximum (penalized) likelihood estimation (Stasinopoulos et al., 2017). Comprehensive reviews of distributional regression techniques have been given by Kneib et al. (2023) and Klein (2024).

Owing to the use of explicit probability distributions for most GNSS positioning algorithms, techniques such as quantile or expectile regression (Koenker & Bassett Jr, 1978; Waltrup et al., 2015) are not considered.

We opt for Bayesian distributional regression because of its flexibility in terms of both model specification and existing software, in addition to natural interpretations of model uncertainty (trust) in Bayesian inference. This class of models can be computationally expensive to estimate when used with large data sets. To achieve significantly improved scaling with data set size, we use stochastic variational inference (SVI; Hoffman et al. (2013)) in place of Markov chain Monte Carlo (MCMC) methods for model estimation. This step is implemented and graphics processing unit (GPU)-accelerated with Python using Pyro, a well-known, powerful deep probabilistic modeling framework (Bingham et al., 2018).

In this work, we consider the modeling of pseudorange measurement errors caused by local multipath as a case study; however, similar methods can be applied to other error sources (e.g., those observed in GNSS carrier-phase measurements, atmospheric delay corrections, etc.). In the following sections, the case study data set is introduced, and the modeling methodology is detailed.

2 DATA

The data set we consider here consists of pseudorange residuals (the response variable) observed by a moving vehicle, in urban and open-sky conditions around central Europe, over > 70 h. Because of the dependence of measurement error characteristics on both GNSS constellation and frequency band, we consider only the 1.3 million pseudorange observations from Global Positioning System (GPS) satellites on the L1 frequency in this case study; however, the modeling methodology discussed herein can be readily applied to any other set of measurements.

Pseudorange residuals contain any error in the pseudorange measurement that is not part of the functional model (e.g., not containing terms such as satellite/receiver clock bias, atmospheric delays, etc.) and primarily consist of local multipath effects. These residuals are estimated offline in large batches using a factor graph estimation framework (Taylor & Gross, 2024). See Appendix A for additional information about the data set’s hardware configuration and estimation process.

To make the pseudorange residuals more amenable to modeling attempts, the original GNSS observations have undergone the following pre-processing steps (in the receiver’s positioning engine):

Application of a strict selection process to discard poor-quality signals, based on removal of the following:
- signals with C/N₀ < 30 dBHz
- all signals obtained when the vehicle is close to stationary (moving at a speed below 1 m/s)
- very-high-variability signals that contribute to inconsistent solutions over a window in time, as measured by window consistency (see Appendix B)
Antenna calibration, to remove systematic biases caused by the antenna platform
Phase smoothing³ of pseudoranges over a symmetric five-epoch (at 1 Hz) window

These exact pre-processing steps are not strictly required for the modeling approach presented here; however, for the characterization of GNSS pseudorange residuals in urban environments using consumer-grade hardware, these steps help to significantly reduce skewness and tail-heaviness in the residual distribution and are thus applied in this case study. Although it is possible to model skewed data (see Section 3 for more details), [not only will skewed models] increase the computational complexity of the modeling process, but, to the authors’ knowledge, skewed distributions have not been considered in measurement error models for GNSS positioning applications in the literature or commercially. However, such considerations are an interesting topic for future work. Furthermore, discarding measurements obtained when the vehicle is static heavily reduces the impact of long-lived, time-correlated biases due to multipath; similarly, the (frequency-, azimuth-, and elevation-dependent) systematic biases due to antenna group delay, which can cause long-term effects, are removed by antenna calibration (Bryant et al., 2020). These additional processes reduce time correlation and spatial (azimuth and elevation) correlation between residuals, allowing for better observability of the underlying statistics and distribution of the data (Zieba, 2010), although phase smoothing does reintroduce some correlation between consecutive measurements.

The time series of residuals from each satellite have similar statistical properties overall (e.g., mean: –0.007 ± 0.009, standard deviation: 0.252 ± 0.016, and autocorrelation), as expected from careful handling of satellite-dependent effects during pre-processing and estimation; hence, we group all GPS L1 observations together into a single data set. It is challenging to rigorously assess stationarity and ergodicity of the data, especially when considering the additional effects of the signal environment. The time series pass standard stationarity tests (i.e., augmented Dickey–Fuller test, Kwiatkowski–Phillips–Schmidt–Shin test), and the time correlation within sequences of pseudorange residuals (periods of continuous phase lock) from a single satellite is relatively low, with their autocorrelation decaying exponentially to below 0.1 after 8 s, on average. Additionally, the overall distribution of pseudorange residuals remains similar over time: the 0.1, 0.5, and 0.9 quantiles of the pseudorange residuals vary with standard deviations of 0.03, 0.01, and 0.03 m, respectively, when split into 50 bins stratified by time.

GNSS measurement error characteristics can be affected by changes to the entire data processing chain, including the antenna platform, receiver signal processing, and any pre-processing applied to the raw GNSS observations. Consequently, regardless of the modeling approach, any positioning algorithm that uses data-driven models should maintain the same or similar data-generating process used in the modeling stage.

To account for the effects of the varied signal environment on the measurements, we model the distribution of the pseudorange residual as a function of two covariates (which may also be referred to here as measurement quality indicators): the signal’s carrier-to-noise-density ratio (C/N₀) and its consistency over time (window consistency), computed in the receiver’s positioning engine. C/N₀ is measured in dBHz, and the window consistency is positive-valued and measured in meters (a lower value indicates improved signal consistency). The window consistency of a signal is computed by using the root mean square of an estimate of a signal’s measurement residuals over a symmetric five-epoch window, obtained by solving for differences in receiver position, clock, and biases over the window, and by using delta-pseudorange/phase measurements and nonlinear (robustified) weighted least squares. See Appendix B for more details.

Figure 1 shows the density of observations over the two-dimensional (2D) covariate space. These covariates were selected because they have the strongest relationship with the distribution of pseudorange residuals – and because using too many covariates increases the covariate space dimensionality and sparsity of data, reducing model interpretability and increasing the difficulty of conditional distribution estimation. Satellite elevation and signal-level metrics, such as correlation function deformation, have been compared against C/N₀ and window consistency using visual diagnostics (e.g., quantile regression curves), but show weaker relationships with the distribution of pseudorange residuals. In the future, it would be interesting to investigate the use and development of more informative quality indicators and to expand the number of quality indicators used in distributional regression, which may require more detailed analysis for covariate selection based on cross-validation or an information criterion (O’Neill & Burke, 2023).

FIGURE 1

Density of observations as a function of measurement quality indicators

3 METHODOLOGY

To specify a Bayesian distributional regression model, we must first choose the distribution of the response variable. Standard commercial GNSS positioning algorithms almost solely use the normal distribution for stochastic models, owing to its useful properties and inexpensive computation. Although Student’s t-distribution presents a significant increase in complexity when used with positioning algorithms, this distribution has also seen some utility in state-of-the-art GNSS applications because of its increased flexibility in modeling heavier-tailed data.

The distribution of pseudorange residuals, conditional on the covariates, can be analyzed by inspecting observations within small regions of covariate space. Analyses using distribution plots show that the data are heavy-tailed, and quantile–quantile plots show that the data can be reasonably well modeled by Student’s t-distribution (see Figure 2 for details on an example). Distributions such as Johnson’s S_U and sinh-arcsinh distributions (Jones & Pewsey, 2009) with the ability to model skewness in addition to tail weight are also potential candidates; however, these distributions are not considered here because of their further increase in complexity beyond Student’s t-distribution for positioning applications.

Comparison of fitted Student’s t-distribution and normal distribution to the pseudorange residual conditional distribution (for C/N0 ∈ [47, 47.5] dBHz and window consistency ∈ [0.24, 0.25] m). (a) Empirical histogram and (fitted) theoretical PDF, (b) Folded cumulative distribution (log10-scale y-axis), (c) Detrended q-q plot (Student’s t-model), (d) Detrended q-q plot (normal model). The fitted Student’s t-distribution has 10 degrees of freedom. A folded CDF (shown in (b)) is a CDF that has been reflected across the line of cumulative density = 0.5. With a log-scale vertical axis, this plot is useful for inspecting the empirical distribution tails, which closely match the theoretical Student’s t-distribution. A quantile–quantile (q–q) plot compares the quantiles of a sample to theoretical quantiles of a given distribution, typically plotted diagonally. The detrended q–q plots in (c) and (d) show sample-theoretical quantiles (hence, quantile deviations) on the vertical axis, for improved resolution (Buuren & Fredriks, 2001), along with confidence intervals that represent the region that a quantile deviation would fall between if the sample came from the theoretical distribution (with a specified probability). Confidence intervals are computed using the distribution of an order statistic (David & Nagaraja, 2004, Chapter 2).

FIGURE 2

Comparison of fitted Student’s t-distribution and normal distribution to the pseudorange residual conditional distribution (for C/N₀ ∈ [47, 47.5] dBHz and window consistency ∈ [0.24, 0.25] m). (a) Empirical histogram and (fitted) theoretical PDF, (b) Folded cumulative distribution (log₁₀-scale y-axis), (c) Detrended q-q plot (Student’s t-model), (d) Detrended q-q plot (normal model).

The fitted Student’s t-distribution has 10 degrees of freedom. A folded CDF (shown in (b)) is a CDF that has been reflected across the line of cumulative density = 0.5. With a log-scale vertical axis, this plot is useful for inspecting the empirical distribution tails, which closely match the theoretical Student’s t-distribution. A quantile–quantile (q–q) plot compares the quantiles of a sample to theoretical quantiles of a given distribution, typically plotted diagonally. The detrended q–q plots in (c) and (d) show sample-theoretical quantiles (hence, quantile deviations) on the vertical axis, for improved resolution (Buuren & Fredriks, 2001), along with confidence intervals that represent the region that a quantile deviation would fall between if the sample came from the theoretical distribution (with a specified probability). Confidence intervals are computed using the distribution of an order statistic (David & Nagaraja, 2004, Chapter 2).

Note that we refer to the three-parameter Student’s t-distribution, which parameterizes not only its degrees of freedom (v), but also location (µ) and scale (σ). In the GAMLSS literature, this parameterization is referred to as the Student t-distribution family and is described by the following probability density function (PDF):

$\begin{matrix} f (x; μ, σ, ν) & = \frac{Γ (\frac{ν + 1}{2})}{Γ (\frac{ν}{2}) σ \sqrt{π ν}} {(1 + \frac{{(x - μ)}^{2}}{σ^{2} ν})}^{- \frac{ν + 1}{2}} \end{matrix}$ 1

where Γ(·) is the gamma function. Based on assumptions on continuity of the response’s conditional distribution, the distributional parameters (µ, σ, v) are modeled as smooth functions of the two covariates. We represent the pseudorange residual, carrier-to-noise-density ratio, and window consistency observations as y_i, x_1i, and x_2i, respectively, for i = 1,⋯, n signals and define the distributional regression model as follows (using notation similar to that of the GAMLSS literature):

$\begin{matrix} y_{i} & \sim Student (μ_{i}, σ_{i}, ν_{i}) \end{matrix}$ 2

$\begin{matrix} μ & = η_{1} \\ \log (σ) & = η_{2} \\ \log (ν) & = η_{3} \end{matrix}$ 3

Here, we use link functions to relate the parameters to the predictors η_k, for k = 1, 2, 3 (using an identity link for μ and log links to constrain σ and v as positive).

We use a class of smoothers⁴ referred to as Bayesian penalized tensor product splines to model the continuous relationship between the predictors η_k and covariates x₁, x₂. These smoothers provide advantages over the simpler choice of polynomials when modeling a regression surface, as they provide flexibility with a reduced risk of nonlocal influence of data on the predictor (Magee, 1998).

3.1 Bayesian Penalized Tensor Product Splines

A Bayesian penalized tensor prsssoduct spline is the result of a combination of different techniques for constructing a smoother.

The foundation of this smoother is a B-spline, as described by Eilers and Marx (2010). A B-spline is a smoothing function of one variable that consists of a sum of B-spline basis functions a_i with weights α_i, for i = 1, …, I:

$\begin{matrix} s (x_{1}) & = \sum_{i = 1}^{I} α_{i} a_{i} (x_{1}) \end{matrix}$ 4

Each a_i is a degree-3 B-spline basis function (with a formulation identical to that of Equation (10) by Eilers and Marx (2010)), centered on each of the I spline knots that evenly span the range of the variable x₁. See Figure 3(a) for a visualization of a B-spline and its basis functions (where I = 20).

FIGURE 3

Illustrations of spline basis functions. (a) Example B-spline and (weighted) basis functions. The 20 weighted basis functions that comprise this B-spline are shown in blue. Each individual curve is a basis function a_i(x₁) multiplied by its weight α_i. The locations of the first three curves are marked (a₁, …, a₃), and the weight α₂ (corresponding to the curve height) of the second curve is marked, (b) Example of a single basis function for a tensor product B-spline

We can generalize the spline for multivariate regression using a “tensor product” construction (Wood, 2006) to create a new basis b_ij, with weights β_ij, from the marginal bases a_1i and a_2j for the two covariates:

$\begin{matrix} s (x_{1}, x_{2}) & = \sum_{i = 1}^{I} \sum_{j = 1}^{J} β_{i j} a_{1 i} (x_{1}) a_{2 j} (x_{2}) \\ = \sum_{i = 1}^{I} \sum_{j = 1}^{J} β_{i j} b_{i j} (x_{1}, x_{2}) \end{matrix}$ 5

Figure 3(b) presents a visualization of an example b_ij. The predictors are modeled as tensor product B-splines and can be written in matrix form:

$\begin{matrix} η_{1} & = Z_{1} β_{1} \\ η_{2} & = Z_{2} β_{2} \\ η_{3} & = Z_{3} β_{3} \end{matrix}$ 6

Here, each Z_k is a basis matrix depending on thess covariates x₁ and x₂, and β_k is a vector containing the basis function weights that we wish to infer.

Eilers and Marx (2010) recommended using many equally spaced spline knots and incorporating a roughness penalty into the regression to tune the smoothness of the resulting function. Incorporating a roughness penalty introduces the assumption of smoothness to the model (which we expect to be realistic for continuous covariates). Notably, this step assists in estimating the difficult-to-observe tail heaviness parameter, v, by allowing more information from nearby datapoints in covariate space to contribute to local parameter estimates.

Because the tensor product spline is constructed from marginal bases, the knots are located on a regular 2D grid across covariate space. If adjacent knots on the grid are considered “neighbors,” the neighboring relationships can be represented on a lattice graph (see Figure 4).

FIGURE 4

Lattice graph 𝒢 = (𝒱, ℰ) with vertices 𝒱, corresponding to tensor product spline knots, and edges ℰ, corresponding to neighbor relationships

The plot is aligned with Cartesian axes to demonstrate the construction of the tensor product spline basis from marginal bases along x₁ and x₂.

For our Bayesian regression approach (where the weights β are considered as latent random variables), we use a stochastic analogue of a roughness penalty, based on a random walk prior proposed for Bayesian P-splines (Lang & Brezger, 2004). It is possible to generalize this concept to the multivariate case for tensor product splines by using a Gaussian Markov random field (GMRF). In essence, a GMRF prior constrains the differences between weights of basis functions with neighboring knots to influence the smoothness of our regression spline.

GMRFs are often applied in spatial statistics, as they are computationally efficient to use in modeling large data sets compared with traditional Gaussian random fields because of their sparsity. Here, we use the GMRF representation for a Matern Gaussian random field described by Lindgren et al. (2011) through the solution of a stochastic partial differential equation (SPDE). This approach allows us to use the flexible and interpretable Matérn field with the computational advantages of a GMRF:

$\begin{matrix} {(k^{2} - Δ)}^{\frac{α}{2}} τ β (u) & = W (u) \\ u \in ℝ^{d}, α & = ν + \frac{d}{2} \end{matrix}$ 7

Here, ∆ is the continuous Laplace operator, and 𝒲(u) is spatial white noise with zero mean and unit variance. The solution to this SPDE is a Matérn Gaussian random field, with parameters τ, κ, v >0. τ affects the precision, κ the spatial scale, and v the smoothness (Mejia et al., 2019). For our 2D problem (d = 2), we set v = 1 and κ = 0.1 and consider τ² a hyperparameter with a (conjugate) gamma-distributed prior. The prior for β is then a multivariate normal distribution:

$\begin{matrix} β | τ & \sim N (0, Q^{- 1}) \\ Q & = τ^{2} {(I + R)}^{⊤} (I + R) \end{matrix}$ 8

where I is the identity matrix and R is the (sparse) Laplacian matrix of the graph representing neighboring knot relations (e.g., Figure 4) (Sidén, 2020). To aid in intuition, the dependence between spline weights in Equation (5), using this GMRF prior, can be written as follows:

$\begin{matrix} (κ^{2} + 4) β_{i, j} & - β_{i - 1, j} - β_{i + 1, j} \\ \begin{matrix} - β_{i, j - 1} - β_{i, j + 1} & \sim N (0, τ^{- 2}) \end{matrix} \end{matrix}$ 9

Note that the lattice graph of knots should extend beyond the range of our data in covariate space, to satisfy GMRF and spline assumptions. For this application, we use 169 knots for each tensor product spline that models a distributional parameter (i.e., 13 knots for a single marginal spline basis).

3.2 Model and Inference

Using this spline formulation for the predictors η_k, for k = 1, 2, 3, we can assemble a Bayesian hierarchical model:

$\begin{matrix} y_{i} | \dots & \sim Student (μ_{i}, σ_{i}, ν_{i}) \\ β_{k} | τ_{k} & \sim N (m_{k}, Q_{k}^{- 1}) \\ τ_{k}^{2} & \sim Gamma (a_{k}, b_{k}) \end{matrix}$ 10

where m_k, a_k, and b_k are fixed, Q_k has a dependency on τ_k as shown in Equation (8), and y_i| … represents the response variable given all of the model’s latent random variables (β₁, β₃,…, τ₁,…, τ₃). The gamma distribution uses scale/shape parameterization, with scale a_k = 2 and shape $b_{k} = \frac{s_{k}}{2}$ , where s_k is equal to 0.05, 0.3, 1 for k = 1, 2, 3, respectively. Using this parameterization, the expected value of the gamma distribution is s_k, which might be interpreted as the “expected variation” of the predictors η_k across the GMRF. From testing on our data set, we have seen that the model estimates are relatively insensitive to these values: when these values vary by ±50%⁵, the model estimates can change by 0.1% – 2.5% for σ and 0.5% – 7% for v (or, respectively, 0.7% and 2.9%, on average), depending on the density of observations (high density corresponds to less sensitivity; see Figure 1 for densities). The model estimate for μ can change by 0.0003 – 0.01 m (or 0.003 m on average).

We can make Bayesian inferences on the latent random variables β_k, τ_k given pseudorange residuals y_i with their posterior distribution. Using Bayes’ theorem, the PDF of the posterior can be written as follows:

$\begin{matrix} p (θ | y) & = \frac{p (y, θ)}{p (y)} = \frac{p (y | θ) p (θ)}{p (y)} \end{matrix}$ 11

where:

$\begin{matrix} p (y) & = \int p (y, θ) d θ \end{matrix}$ 12

Here, the latent random variables (β₁, τ₁, …) in the model are collected into a single vector θ for ease of notation.

To summarize, given our pseudorange residual data, we wish to infer the distribution of θ (which encodes the relationship between covariates and the distribution of the pseudorange residuals), denoted by the posterior p(θ|y). Our model specifies the conditional distribution of the residuals p(y|θ) with Student’s t-distribution, as well as the prior information p(θ) avaisslable for the latent random variables (corresponding to smoothness assumptions of the spline).

Monte Carlo methods such as MCMC are typically used to provide estimates for integrals of the posterior (for example, to compute the median/expected value of θ|y) because, for many models, including ours, these integrals are unavailable analytically. For our model and data set size, this approach is computationally infeasible. Instead, we use variational inference, a well-known method for approximate Bayesian inference (Blei et al., 2017).

3.3 Variational Inference

The idea behind variational inference is to approximate the posterior with a more tractable variational distribution. This approximation is achieved by finding a probability distribution, with density q_ϕ(θ), within a variational family⁶ parameterized by ϕ that best approximates the posterior. The measure of similarity between the distributions is measured using the Kullback-Leibler (KL) divergence, resulting in an objective for the optimal variational distribution q_ϕ*(θ), where the following holds:

$\begin{matrix} ϕ^{*} & = \arg \min_{ϕ} KL (q_{ϕ} (θ) | | p (θ | y)) \end{matrix}$ 13

Because of the often analytically intractable p(y) term (referred to as evidence in the Bayesian context) contained in the posterior, the KL divergence is reformulated as follows:

$\begin{matrix} KL (q_{ϕ} (θ) | | p (θ | y)) & = E_{q_{ϕ} (θ)} [\log q_{ϕ} (θ) - \log p (θ | y)] \\ = E_{q_{ϕ} (θ)} [\log q_{ϕ} (θ) - \log p (y, θ) + \log p (y)] \\ = \log p (y) - E_{q_{ϕ} (θ)} [\log p (y, θ) - \log q_{ϕ} (θ)] \\ = \log p (y) - ELBO (ϕ) \end{matrix}$ 14

where 𝔼_{qϕ (θ)}[·] is the expectation with respect to q_ϕ (θ). The optimization problem is then reframed as the maximization of the evidence lower bound (ELBO), which is equivalent to Equation (13) because the KL divergence is non-negative:

$\begin{matrix} ϕ^{*} & = \arg \max_{ϕ} ELBO (ϕ) \\ ELBO (ϕ) & = E_{q_{ϕ} (θ)} [\log p (y, θ) - \log q_{ϕ} (θ)] \end{matrix}$ 15

This optimization-based technique for Bayesian inference scales to large data sets or complex models significantly better than sampling methods (Blei et al., 2017).

Many popular, off-the-shelf probabilistic programming tools are available for performing variational inference on Bayesian models, including Stan (Kucukelbir et al., 2015), PyMC (Oriol et al., 2023), Edward (Tran et al., 2016), and Pyro (Bingham et al., 2018). We use Pyro for its flexibility, control, and familiarity (with Python/PyTorch).

There are several different methods for maximizing the ELBO, depending on both the model and variational family. Pyro’s SVI algorithm uses a form of stochastic gradient ascent to maximize the ELBO, using unbiased Monte Carlo estimates of the ELBO gradient (Pyro, 2017).

Additionally, we take advantage of subsampling in Pyro’s SVI, which involves using subsamples of the data set to estimate the likelihood and gradients. This approach can significantly reduce the memory footprint and enhance computational efficiency when working with large data sets, alongside improved global convergence properties typically associated with stochastic gradient algorithms (Bottou, 1991).

3.4 Variational Family

A variational family 𝒟 is a set of PDFs parameterized by ϕ, over our latent random variables θ, from which we can select a variational distribution q_ϕ(θ) ~ 𝒟. We use a set of multivariate normal distributions, where ϕ includes means and covariances. In particular, we use a set of “low-rank” multivariate normal distributions, with factored covariance matrices:

$\begin{matrix} Σ & = W W^{⊤} + D \end{matrix}$ 16

where W is a matrix with significantly fewer columns than rows and D is a diagonal matrix. The use of this low-rank covariance factorization is a compromise between the inflexibility of a diagonal covariance and the quadratic scaling of the number of free parameters in a full-rank covariance with the number of latent random variables. In addition, the factorization enables reduced inverse and determinant computation, which, with the use of matrix identities, scales with the number of columns of W. Because we expect to capture the correlation between the components of the penalized tensor product splines described in Section 3.1 and the relationship between the σ and ν parameters in Student’s t-distribution, some representation of multivariate dependence in the posterior is necessary.

3.5 GPU Acceleration

Because this type of optimization problem incorporates many single-instruction multiple-data operations, we use a GPU to accelerate computation. This further (and vastly) improves the scaling of computation time with data set size (see Table 1 for a rough comparison).

View this table:

TABLE 1 Benchmarks of 100 SVI Iterations for a CPU and GPU with Different Numbers of Observations, from a Synthetic Data Set (Tested on an AWS EC2 g5.2xlarge Instance)

4 SYNTHETIC DATA STUDY

To test the methodology and understand what we can expect from the model, we examine results from a synthetic data set with properties similar to those of the target data set. The synthetic data set consists of 1 million Student’s t-distributed random observations, y_i, with the distributional parameters (μ, σ, and v) varying as a function of two covariates x₁ and x₂. The values of the covariates x₁ and x₂ are generated from a multivariate normal distribution with mean $[\begin{matrix} 45 \\ 0.4 \end{matrix}]$ and covariance $[\begin{matrix} 5 & 0 \\ 0 & 0.1 \end{matrix}]$ , to emulate the characteristic of observations being clustered in a single region of covariate space. Then, the response variable y_i is generated from a Student’s t-distribution conditional on x₁ and x₂:

$\begin{matrix} y_{i} & \sim Student (μ_{i}, σ_{i}, ν_{i}) \\ μ & = 0.3 z_{1} \\ σ & = 0.1 \exp (0.5 + 0.2 (z_{2} - z_{1}) + 0.01 (z_{1}^{2} + z_{2}^{2})) \\ ν & = 10 + \exp (2.2 - z_{1}^{2} - 0.4 z_{2}^{2}) \end{matrix}$ 17

where $z_{1} = \frac{1}{10} (x_{1} - 45)$ and z₂ = 5(x₂ – 0.4). These values were chosen to emulate the experimental data; therefore, y_i, x₁, and x₂ may be considered with the same units described in Sections 2 and 3 (meters, dBHz, and meters, respectively).

In Figure 5, the true distributional parameters are visualized as surfaces, with the fitted model’s credible intervals for the parameters overlaid. We can see that the model’s credible interval surfaces encapsulate and tightly match the true values of the μ and σ parameters, even toward the edges of the surfaces, where the observations are very sparse. However, this tight matching is not observed for v, where the intervals are larger and less closely follow the shape of the true surface.

FIGURE 5

Synthetic data: (a–c) Surface plots representing the relationship between the value of Student’s t-distributional parameters μ, σ, v, on the vertical axis and the covariates x₁ and x₂ are shown. The multicolored surface represents the true values of the conditional distribution parameters for the synthetically generated data, and the surrounding transparent light-blue surfaces represent 99% credible intervals (CIs) from the fitted model. (d–f) To aid in interpretation of the surface plots, we also show a 2D slice visualizing the relationship between the parameters and x₁ only, for a fixed value of x₂ = 0.4. Note that the surface is drawn only in the region of covariate space where the majority (99.9%) of observations lie.

To give an indication of the scales shown in Figure 5, the median widths of the 99% credible intervals (where we assign more weight to regions of covariate space with more frequent observations) are 0.005, 0.003, and 3.1, for μ, σ, and v, respectively. Because v represents the degrees of freedom, a larger estimate uncertainty is expected from the difficulty in observing the conditional distribution’s tail. When generating different random realizations of the synthetic data set, the 99% credible interval surfaces generally encapsulate the true surface over the entire covariate space; however, in some rare circumstances, the credible interval surface does not bound the truth in small regions (see Figure 6), as expected to some degree owing to the randomness in the data.

FIGURE 6

Synthetic data: (a) Surface plot demonstrating a rare example of the v parameter credible intervals intersecting with the true surface in a small region (x₁ ≈ 48, x₂ ≈ 0.5); (b) the corresponding slice for x₂ = 0.5

The ability of the model to infer v with reasonable accuracy, even in sparse regions of covariate space, comes from the model construction and smoothness assumptions (Section 3.1). Improvements on these estimates for limited data would require stricter assumptions, based on prior knowledge of the structure of the relationship between the covariates and distribution parameters.

5 RESULTS

Here, we analyze the model fitted to the target data set described in Section 2 with the aid of diagnostic plots. To obtain final estimates of the model parameters, we use the median of the (approximated) Bayesian posterior.

In Figure 7, the distribution parameter estimates from the fitted model are visualized as surfaces dependent on the covariates x₁ and x₂ (carrier-to-noise-density ratio in dBHz and window consistency in meters, respectively). We observe that the parameters have clear relationships with both the carrier-to-noise-density ratio and the window consistency. The scale parameter shows a smooth variation, expectedly increasing with both signal noise and inconsistency. The relationship between the degrees of freedom, v, and the covariates is less smooth, and we observe wider credible intervals.

FIGURE 7

(a–c) Surface plots representing the relationship between the value of fitted Student’s t-distributional parameters μ, σ, v and the covariates (carrier-to-noise-density ratio and window consistency) are shown. The multicolored surface represents the estimate (corresponding to the posterior median) of the conditional distribution parameter, and the surrounding transparent light-blue surfaces represent 99% credible intervals of the estimate. (d–f) To aid in interpretation of the surface plots, we also show a 2D slice visualizing the relationship between the parameters and C/N₀ only, for a fixed value of window consistency = 0.45. Note that the surface is drawn only in the region of covariate space where the majority (99.9%) of observations lie.

The location of the conditional distribution, μ, has a pronounced dependence on the carrier-to-noise-density ratio. These higher-frequency effects may be driven by receiver characteristics or antenna group delay biases (Caizzone et al., 2022; Kaufmann & Julien, 2021) and are successfully captured by the model.

The model’s quantified uncertainty for the conditional distribution parameters gives an indication of the accuracy we might expect when using these estimates. The median widths of the 99% credible intervals (where we assign more weight to regions of covariate space with more frequent observations) are 0.005, 0.004, and 1.1, for μ, σ, and v, respectively.

In Figure 8, the goodness of fit of the conditional distribution is assessed with detrended q–q plots within several small regions of covariate space (see Figure 2(c) for notes on interpretation of a q–q plot). It is clear that the fitted model adapts well to the varying distribution of pseudorange residuals across covariate space, although the data are not always perfectly described by Student’s t-distribution (which cannot model some slight skewness that is visible in the data).

FIGURE 8

Detrended q–q plots of data within nine bins in covariate space, compared with the distribution specified by the fitted model at the coordinates of each bin center (shown at the top of each subfigure)

The bin centers were chosen using combinations of 0.1, 0.5, and 0.9 quantiles for x₁ and x₁, and the bin widths were chosen such that there are approximately 5000 samples within each bin. To assess the goodness of fit of the conditional distribution, we examine how closely the points in the q–q plot fall between the confidence intervals.

To more succinctly summarize the goodness of fit of the model, we use normalized quantile residuals (Stasinopoulos et al., 2017), which we briefly describe here. Suppose that $F (y; \hat{μ}, \hat{σ}, \hat{v})$ is the model CDF. Then, if the model is correct, the quantile residuals û_i = F(y_i) for observations y_i (i = 1,…, n) have a uniform distribution between zero and one (the probability integral transform). Normalized quantile residuals r_i = Φ^–1(u_i), where Φ^–1 is the inverse CDF of the standard normal distribution, have a standard normal distribution for a perfect model.

In Figure 9(a), we can see that the normalized quantile residuals have a distribution close to the standard normal. Further, Figure 9(b) presents a visual summary of the goodness of fit of the model across covariate space, with stratified normalized quantile residuals. The Shapiro-Wilk normality test (Shapiro & Wilk, 1965) is used to color a region of covariate space, based on the normalized quantile residuals observed in the corresponding rectangular bins. The p-value associated with the test statistic is used to assess the normality of the normalized quantile residuals. To construct this figure, because the sample size affects the results of the Shapiro-Wilk test, a subsample of 500 observations is used within each bin to limit the effect of differing observation frequency over covariate space (Figure 1). To reduce the variance of the test results, the p-values are recomputed 40 times using new subsamples (with replacement), and the median p-value is used to color each bin in the figure.

FIGURE 9

Summary plots indicating the model’s goodness of fit using normalized quantile residuals. (a) Histogram of the normalized quantile residuals compared with a PDF of the standard normal distribution, (b) Goodness of fit of the model across covariate space, measured using the Shapiro–Wilk normality test on normalized quantile residuals. The p-value represents the probability that the sample of normalized quantile residuals came from a normally distributed population (i.e., purple gives evidence toward a low goodness of fit), (c) Skewness of the pseudorange residuals over covariate space.

The model achieves a good fit over most of covariate space; however, the smaller p-values shown in Figure 9(b) (indicating a lower probability that the normalized quantile residuals are normally distributed) give evidence toward a poorer fit of the model in the outer regions of covariate space. If, in Figure 9(b), we define a bin with a p-value⁷ greater than 0.05 as a “good fit,” then 98% of the data are within a good-fit region of covariate space (yellow), 1% of data are in a poor-fit region (dark purple), and the remaining 1% of data occupy the region of covariate space where observations are too sparse to analyze (uncolored).

Based on deeper analysis using quantile-quantile plots and a similar stratified figure (Figure 9(c)), the regions of poor fit can be primarily attributed to some skewness that cannot be modeled by the symmetric Student’s t-distribution. A simple solution to limit the effects of this mismodeling in a positioning system is to discard the 1% of measurements observed with quality indicators within these poorly fitting regions, alongside the additional 1% considered as unmodeled owing to sparsity. Alternatively, this skewness may be reduced in the data-generating process (e.g., with changes to the receiver’s positioning engine) or accounted for (as discussed earlier in Section 3) by using a more complex conditional distribution.

Overall, the flexibility from the tensor product splines in the model specification allows the Student’s t-distribution parameters to shift adequately to fit the conditional distribution of the data, providing a good basis for a local multipath error model in the majority of environments observed in the automotive data set. For use of these models in safety-critical applications, further verification based on large data sets and theory for relevant goodness-of-fit tests directly relating to integrity requirements will likely need to be developed in future work.

6 MODEL UNCERTAINTY IN GNSS

Here, we proceed beyond a basic interpretation of model trust to describe the potential use of uncertainty quantification in stochastic models as supplementary information for GNSSs.

GNSSs provide position estimates and position bounds based on measurements and assumed models. Unfortunately, it is well known that all models are wrong; that is, a statistical or dynamical model will not perfectly describe a real system’s data-generating process. Consequently, although it is clear that this is a second-order effect beyond accurate estimation of models, quantifying and compensating for the impact of the “wrongness” of these models in GNSSs (see Figure 10) may be beneficial for robust position estimation and required for safety-critical systems.

FIGURE 10

Example illustration of the sensitivity of position bounds arising from model uncertainty within a positioning system

An uncertain model used in a positioning algorithm will result in uncertainty in the algorithm’s outputs.

As described in this paper, the Bayesian framework enables a quantification of uncertainty in fitted stochastic models where, by using the variational inference approach taken in Section 3.4, measurement error model parameters ψ can be represented as jointly multivariate normally distributed random variables:

$ψ \sim N (ψ^{*}, Σ_{ψ})$ 18

where ψ* is the best estimate of the model parameters (corresponding to the approximate posterior mean/median) and ∑_ψ is their covariance describing uncertainty in the estimate.

In a typical GNSS, only the best estimate of the measurement error model parameters, ψ*, would be used. However, the additional information of the uncertainty in measurement error models is an interesting candidate to be used. This additional information could be used to provide a more complete stochastic model for more informed position estimates and bounds – i.e., treating ψ as a set of random variables within a suitably flexible positioning system (e.g., factor graph optimization or related Bayesian filtering methods), conceptually similar to Bayesian hierarchical modeling (Gelman et al., 2014). This additional information could also be used to assess the impact of this uncertainty on a system using a form of sensitivity analysis.

With sensitivity analysis, a system’s outcomes as a response to model uncertainty can be analyzed. This approach can be formalized by writing the system f as a function of its model parameters ψ, with system outputs z:

$\begin{matrix} z & = f (ψ) \end{matrix}$ 19

Using the information on parameter uncertainty (ψ are normally distributed; see Equation (18)), we can conduct a local sensitivity analysis using the gradient of f. This analysis can be achieved by simply using error propagation (assuming local linearity of f): the uncertainty in the system output is given by J∑_ψJ^⊤, where J is the Jacobian of f evaluated at ψ*. More advanced global sensitivity analysis methods also benefit from the measure of uncertainty gained from Bayesian methods (Saltelli et al., 2008).

However, a recognizable shortcoming of sensitivity analysis based on parametric models is that the parameterization of the stochastic models is fixed and assumed correct – for example, in the case study presented here, the conditional distribution of the pseudorange residuals and its parametric relationship with the covariates is assumed. Nevertheless, if the basic principles of the modeling methods and analysis described in this paper are developed further in the future, perhaps to improve upon the assessment of model goodness of fit in relation to GNSS integrity requirements and strict validation of models, such uncertainty quantification may be beneficial for safety-critical systems (as the quality and knowledge of uncertainty in models used for these systems are paramount) – for example, to measure the sensitivity of reported position bounds for a given confidence level.

7 CONCLUSION

In this paper, we have presented an approach for modeling GNSS measurement errors using distributional regression and quantified the model uncertainty with Bayesian inference. A large automotive data set consisting of GNSS pseudorange measurement residuals was used to demonstrate the proposed model and estimation framework.

By using variational inference for model estimation, we were able to scale the data set size and model complexity such that nonlinear relationships between distributional parameters and measurement quality indicators could be reliably estimated, while overcoming the computational limitations of typical Bayesian inference with MCMC methods.

The Bayesian framework used herein enabled direct means for uncertainty quantification (e.g., with credible intervals), in contrast to “crude” standard error approximations (Hastie, 1992) for comparable frequentist models in the GAMLSS software suite. This feature is primarily useful for interpreting model trust directly from the data. Potential applications in sensitivity analysis and detailed modeling of GNSSs, which rely on such uncertain stochastic models, were discussed (and reserved for future investigation) to improve the robustness and confidence of computed position estimates and bounds.

The model we presented for the case study characterizing the distribution of pseudorange residuals (conditioned on the signal carrier-to-noise-density ratio and window consistency) was able to fit the data well, including estimating difficult-to-observe distribution tail weights, even in signal environments with limited data. It was observed that 98% of measurements in the data set were able to be well modeled, with the remaining 2% of data unmodeled because of skewness or sparsity being identified and removed. Detailed goodness-of-fit analysis and a synthetic study showed that the model framework was effective in adapting to the changing conditional distribution across ranging values of C/N₀ and window consistency.

For a system with a similar data-generating process (i.e., GNSS hardware, firmware, and operating environment), the fitted distributional regression model developed in the case study can be readily implemented using a small set of regression parameters to make predictions for the distribution of pseudorange residuals.

There is scope to extend this work in GNSS error modeling for many different sources of error (which can be conditioned on environmental and related factors), including local multipath on carrier-phase measurements, atmospheric corrections, or orbit/clock corrections. The conditional distribution assumption for error may be improved by using a distribution that can model skewness as well as tail weight (as briefly discussed in Section 3) or by applying wrapped distributions (Pewsey et al., 2007) for data with a periodic domain such as carrier-phase measurement errors.

HOW TO CITE THIS ARTICLE:

Dorahy, H., Hide, C., Julien, O., Martini, I., & Sheret, I. (2026). Scalable Bayesian distributional regression for GNSS multipath error modeling with uncertainty quantification. NAVIGATION, 73. https://doi.org/10.33012/navi.755

ACKNOWLEDGMENTS

Many thanks to the team at Delft University of Technology, in collaboration for the “I-GNSS positioning for assisted and automated driving” project (18305, funded by the Dutch Research Council), for insightful discussions that contributed to this work.

APPENDICES

A ADDITIONAL DATA SET DETAILS

In this section, additional information about the case study data set described in Section 2 is given.

A.1 Hardware Configuration

Vehicle: Mercedes-Benz Vito (2017)
Antenna: u-blox ANN-MB multi-band antenna
- Antenna mounted on the vehicle on a flat metal roof rack toward the driver’s side edge
Receiver: u-blox ZED-F9P high-precision multi-band GNSS receiver
- Using modified firmware with a fixed discriminator and modified bandwidth

A.2 Hardware Configuration

Pseudorange residuals (the part of the pseudorange measurement that is not accounted for by the functional model) are estimated offline in large batches using a factor graph estimation framework that solves for the entire system of measurements made by the receiver over a few hours, including those from the following bands: GPS L1, L2C; Galileo E1C, E5b; BeiDou B1ID1, B2ID1; GLONASS L1OF, L2OF. To aid in the accurate estimation of residuals, the following steps are applied:

Carrier-phase ambiguities are fixed to integers, enabling the use of precise range measurements with millimeter-level measurement errors.
Virtual reference station corrections, using a dense network of reference receivers, are applied to minimize ionosphere and troposphere delays.
- Reference station pseudorange multipath is mitigated through the use of divergence-free carrier-phase smoothing.
Precise satellite models, from the GFZ German Research Centre for Geosciences, are applied to minimize satellite errors.
Change-in-position estimates from an Applanix reference system using a high-performance inertial measurement unit are used to improve estimation robustness. Absolute position estimates were not used.

B WINDOW CONSISTENCY

The window consistency signal quality indicator utilized in this work is based on a nonlinear (robustified) weighted least-squares solution over a symmetric five-epoch (at 1 Hz) window around the epoch of interest. Delta-pseudorange and delta-phase measurements across the window, $\underline{y} \in ℝ^{m}$ , are used to estimate differences in receiver position, clock, and biases, $\underline{x} \in ℝ^{n}$ , over the window using a suitable linearized observation equation $\underline{y} = A \underline{x}$ , where A is an m × n matrix. The following solution is estimated using a robustified Gauss-Newton algorithm (Triggs et al., 1999):

$\underline{\hat{x}} = \arg \min \sum_{i = 1}^{m} L_{δ} (w_{i} z_{i})$ 20

where z_i are observation model residuals (elements of $\underline{z} = \underline{y} - A \underline{x}$ ), w_i are scalar measurement weights, and L_δ is the Huber loss function (Huber, 1964), for which we set the value of the quadratic/linear change point, δ, to 2. The window consistency, for each signal from each satellite, is then the root mean square of the measurement residuals over the window.

Quantifying the variability in measurements over a time window gives us a direct indicator of the signal’s quality, which is valuable for characterizing the stochastic properties of the measurement error. Note that the time window we have utilized here is symmetric, as it provides a more informative indicator than an asymmetric window, meaning that a GNSS that uses window consistency must handle measurement error models that come with a 2-s latency.

Footnotes

↵1 A response variable is the response we wish to measure/model dependent on related factors, also known as a dependent variable.
↵2 A covariate is a variable that affects the outcome of the response, also known as an independent variable.
↵3 Combined use of carrier phase with pseudorange observations to reduce noise in pseudorange measurements.
↵4 A method for fitting a smooth, continuous function to data.
↵5 Independent, uniformly randomly chosen variations for each parameter, over 60 tests.
↵6 A space of probability distributions, e.g., all multivariate normal distributions (see Section 3.4).
↵7 A log10 p-value greater than –1.3.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

REFERENCES

↵
Al Hage, J., Salvatico, N., & Bonnifait, P. (2023). High integrity localization of intelligent vehicles with Student’s t filtering and fault exclusion. In 2023 IEEE 26th International Conference on Intelligent Transportation Systems (ITSC), 1341–1348. https://doi.org/10.1109/ITSC57777.2023.10422598
Google Scholar
↵
Al Hage, J., Xu, P., Bonnifait, P., & Ibanez-Guzman, J. (2020). Localization integrity for intelligent vehicles through fault detection and position error characterization. IEEE Transactions on Intelligent Transportation Systems, 23(4), 2978–2990. https://doi.org/10.1109/TITS.2020.3027433
Google Scholar
↵
Bingham, E., Chen, J. P., Jankowiak, M., Obermeyer, F., Pradhan, N., Karaletsos, T., Singh, R., Szerlip, P., Horsfall, P. & Goodman, N. D. (2018). Pyro: Deep universal probabilistic programming. Journal of Machine Learning Research. http://jmlr.org/papers/v20/18-403.html
Google Scholar
↵
Blei, D. M., Kucukelbir, A., & McAuliffe, J. D. (2017). Variational inference: A review for statisticians. Journal of the American Statistical Association, 112(518), 859–877. https://doi.org/10.1080/01621459.2017.1285773
Google Scholar CrossRef
↵
Bottou, L. (1991). Stochastic gradient learning in neural networks. Proceedings of Neuro-Nimes, 91(8), 12
Google Scholar
↵
Bryant, R., Julien, O., Hide, C., Moridi, S., & Sheret, I. (2020). Novel snapshot integrity algorithm for automotive applications: Test results based on real data. In Proceedings of the 2020 IEEE/ION Position, Location and Navigation Symposium (PLANS), 670–681. https://doi.org/10.1109/PLANS46316.2020.9109830
Google Scholar
↵
Bryant, R., Julien, O., Hide, C., Skorupa, M., Dorahy, H., & Sheret, I. (2021). Safety-critical automotive positioning based on SEPB without atmospheric corrections. In Proc. of the 34th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2021), 1843–1858. https://doi.org/10.33012/2021.17985
Google Scholar
↵
Buuren, S. v., & Fredriks, M. (2001). Worm plot: A simple diagnostic device for modelling growth reference curves. Statistics in Medicine, 20(8), 1259–1277. https://doi.org/10.1002/sim.746
Google Scholar CrossRef PubMed Web of Science
↵
Caizzone, S., Circiu, M.-S., Elmarissi, W., Enneking, C., Rippl, M., & Sgammini, M. (2022). The role of antennas on GNSS pseudorange and multipath errors and their impact on DFMC multipath models for avionics. NAVIGATION, 69(3) https://doi.org/10.33012/navi.532
Google Scholar
↵
David, H. A., & Nagaraja, H. N. (2004) Order statistics( 3rd ed.). John Wiley & Sons https://doi.org/10.1002/9781118445112.stat00830
Google Scholar
↵
DeCleene, B. (2000). Defining pseudorange integrity-overbounding. Proceedings of the 13th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GPS 2000), 1916–1924. https://www.ion.org/publications/abstract.cfm?articleID=1603
Google Scholar
↵
Eilers, P. H., & Marx, B. D. (2010). Splines, knots, and penalties. Wiley Interdisciplinary Reviews: Computational Statistics, 2(6), 637–653. https://doi.org/10.1002/wics.125
Google Scholar
↵
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014) Bayesian data analysis( 3rd ed.). Chapman;Hall/CRC
Google Scholar
↵
Groves, P. D., & Jiang, Z. (2013). Height aiding, C/N0 weighting and consistency checking for GNSS NLOS and multipath mitigation in urban areas. The Journal of Navigation, 66(5), 653–669. https://doi.org/10.1017/S0373463313000350
Google Scholar
↵
Hastie, T. J. (1992)Generalized additive models.In M. Chambers, J., & J. Hastie, T.(Eds.), Statistical models in S(pp. 303–304). Chapman & Hall/CRC https://doi.org/10.1201/9780203738535
Google Scholar
↵
Hoffman, M. D., Blei, D. M., Wang, C., & Paisley, J. (2013). Stochastic variational inference. Journal of Machine Learning Research. https://doi.org/10.48550/arXiv.1206.7051
Google Scholar
↵
Huang, Y., Zhang, Y., Li, N., & Chambers, J. (2016). Robust Student’s t based nonlinear filter and smoother. IEEE Transactions on Aerospace and Electronic Systems, 52(5), 2586–2596. https://doi.org/10.1109/TAES.2016.150722
Google Scholar
↵
Huber, P. J. (1964). Robust estimation of a location parameter. The Annals of Mathematical Statistics, 35(1), 73–101. https://doi.org/10.1214/aoms/1177703732
Google Scholar
↵
Jones, M. C., & Pewsey, A. (2009). Sinh-arcsinh distributions. Biometrika, 96(4), 761–780. https://doi.org/10.1093/biomet/asp053
Google Scholar CrossRef
↵
Kaufmann, T., & Julien, O. (2021). Automotive GNSS antenna platform-induced group delay biases on GPS L1 signals. 2021 15th European Conference on Antennas and Propagation (EuCAP), 1–5. https://doi.org/10.23919/EuCAP51087.2021.9411115
Google Scholar
↵
Khanafseh, S., Kujur, B., Joerger, M., Walter, T., Pullen, S., Blanch, J., Doherty, K., Norman, L., de Groot, L., & Pervan, B. (2018). GNSS multipath error modeling for automotive applications. Proceedings of the 31st International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2018), 1573–1589. https://doi.org/10.33012/2018.16107
Google Scholar
↵
Klein, N. (2024). Distributional regression for data analysis. Annual Review of Statistics and Its Application, 11 https://doi.org/10.1146/annurev-statistics-040722-053607
Google Scholar
↵
Kneib, T., Silbersdorff, A., & Säfken, B. (2023). Rage against the mean—a review of distributional regression approaches. Econometrics and Statistics, 26,99–123. https://doi.org/10.1016/j.ecosta.2021.07.006
Google Scholar
↵
Koenker, R., & Bassett Jr, G. (1978). Regression quantiles. Econometrica: Journal of the Econometric Society, 33–50. https://doi.org/10.2307/1913643
Google Scholar
↵
Kucukelbir, A., Ranganath, R., Gelman, A., & Blei, D. (2015). Automatic variational inference in Stan. Advances in Neural Information Processing Systems, 28 https://doi.org/10.48550/arXiv.1506.03431
Google Scholar
↵
Lang, S., & Brezger, A. (2004). Bayesian P-splines. Journal of Computational and Graphical Statistics, 13(1), 183–212. https://doi.org/10.1198/1061860043010
Google Scholar
↵
Lindgren, F., Rue, H., & Lindström, J. (2011). An explicit link between Gaussian fields and Gaussian Markov random fields: The stochastic partial differential equation approach. Journal of the Royal Statistical Society Series B: Statistical Methodology, 73(4), 423–498. https://doi.org/10.1111/j.1467-9868.2011.00777.x
Google Scholar CrossRef Web of Science
↵
Magee, L. (1998). Nonlocal behavior in polynomial regressions. The American Statistician, 52(1), 20–22. https://doi.org/10.2307/2685560
Google Scholar CrossRef Web of Science
↵
Mejia, A. F., Yue, Y. R., Bolin, D., Lindgren, F., & Lindquist, M. A.. (2019). A Bayesian general linear modeling approach to cortical surface fMRI data analysis. Journal of the American Statistical Association. https://doi.org/10.1080/01621459.2019.1611582
Google Scholar
↵
Navarro, P. F. (2016). Method for computing an error bound of a Kalman filter based GNSS position solution (European Patent No. EP3009860A1). European Patent Office. https://patents.google.com/patent/EP3009860A1/en
Google Scholar
↵
O’Neill, M., & Burke, K.(2023). Variable selection using a smooth information criterion for distributional regression models. Statistics and Computing, 33(3), 71 https://doi.org/10.1007/s11222-023-10204-8
Google Scholar
↵
Oriol, A.-P., Virgile, A., Colin, C., Larry, D., J,, C, F., Maxim, K., Ravin, K., Jupeng, L., C, C, L., A,, O, M., Michael, O., Ricardo, V., Thomas, W., & Robert, Z. (2023). Pymc: A modern and comprehensive probabilistic programming framework in Python. PeerJ Computer Science, 9,e1516. https://doi.org/10.7717/peerj-cs.1516
Google Scholar
↵
Pewsey, A., Lewis, T., & Jones, M. (2007). The wrapped t family of circular distributions. Australian & New Zealand Journal of Statistics, 49(1), 79–91. https://doi.org/10.1111/j.1467-842X.2006.00465.x
Google Scholar
↵
Pyro,. (2017) SVI part III.Retrieved October 7, 2023, from https://pyro.ai/examples/svi_part_iii.html
Google Scholar
Rife, J., Pullen, S., Enge, P., & Pervan, B. (2006). Paired overbounding for nonideal LAAS and WAAS error distributions. IEEE Transactions on Aerospace and Electronic Systems, 42(4), 1386–1395. https://doi.org/10.1109/taes.2006.314579
Google Scholar
↵
Rigby, R. A., & Stasinopoulos, M. D. (2005). Generalized additive models for location, scale and shape. Journal of the Royal Statistical Society Series C: Applied Statistics, 54(3), 507–554. http://www.jstor.org/stable/3592732
Google Scholar CrossRef
↵
Roth, M., özkan, E., & Gustafsson, F. (2013). A Student’s t filter for heavy tailed process and measurement noise. 2013 IEEE International Conference on Acoustics, Speech and Signal Processing, 5770–5774. https://doi.org/10.1109/ICASSP.2013.6638770
Google Scholar
↵
Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., & Tarantola, S. (2008). Global sensitivity analysis: The primer.John Wiley & Sons https://doi.org/10.1002/9780470725184
Google Scholar
↵
Shapiro, S. S., & Wilk, M. B. (1965). An analysis of variance test for normality (complete samples). Biometrika, 52(3/4), 591–611. https://doi.org/10.2307/2333709
Google Scholar CrossRef PubMed Web of Science
↵
Sidén, P. (2020). Scalable Bayesian spatial analysis with Gaussian Markov random fields (Vol. 790). Linköping University Electronic Press. https://doi.org/10.3384/diss.diva-165872
Google Scholar
↵
Stasinopoulos, M. D., Rigby, R. A., Heller, G. Z., Voudouris, V., & De Bastiani, F.. (2017). Flexible regression and smoothing: Using GAMLSS in R.CRC Press
Google Scholar
↵
Taylor, C., & Gross, J.(2024). Factor graphs for navigation applications: A tutorial. NAVIGATION, 71(3) https://doi.org/10.33012/navi.653
Google Scholar
↵
Teunissen, P. J., & Montenbruck, O. (2017) Springer handbook of global navigation satellite systems (Vol. 10).Springer https://doi.org/10.1007/978-3-319-42928-1
Google Scholar
↵
Tran, D., Kucukelbir, A., Dieng, A. B., Rudolph, M., Liang, D., & Blei, D. M. (2016). Edward: A library for probabilistic modeling, inference, and criticism. arXiv. https://doi.org/10.48550/arXiv.1610.09787
Google Scholar
↵
Triggs, B., McLauchlan, P. F., Hartley, R. I., & Fitzgibbon, A. W. (1999). Bundle adjustment—a modern synthesis. International Workshop on Vision Algorithms, 298–372. https://doi.org/10.1007/3-540-44480-7_21
Google Scholar
↵
Waltrup, L. S., Sobotka, F., Kneib, T., & Kauermann, G. (2015) Expectile and quantile regression – David and Goliath? Statistical Modelling, 15(5), 433–456. https://doi.org/10.1177/1471082x14561155
Google Scholar
↵
Welte, A., Xu, S., & Bonnifait, P. (2017) Protection levels for high integrity localization for autonomous driving [Master’s thesis, ENSTA Bretagne, Univ. Angers]. https://www.hds.utc.fr/~bonnif/supervision/Report_2017_Anthony_Welte_Final.pdf
Google Scholar
↵
Wood, S. N. (2006). Low-rank scale-invariant tensor product smooths for generalized additive mixed models. Biometrics, 62(4), 1025–1036. https://doi.org/10.1111/j.1541-0420.2006.00574.x
Google Scholar CrossRef PubMed Web of Science
↵
Wörner, M., Schuster, F., Dölitzscher, F., Keller, C. G., Haueis, M., & Dietmayer, K. (2016). Integrity for autonomous driving: A survey. Proceedings of the 2016 IEEE/ION Position, Location and Navigation Symposium (PLANS), 666–671. https://doi.org/10.1109/plans.2016.7479759
Google Scholar
↵
Zieba, A. (2010). Effective number of observations and unbiased estimators of variance for autocorrelated data-an overview. Metrology and Measurement Systems, (1), pp. 3–16. http://www.metrology.pg.gda.pl/full/2010/M&MS_2010_003.pdf
Google Scholar

PDF

Previous Next

Scalable Bayesian Distributional Regression for GNSS Multipath Error Modeling with Uncertainty Quantification

Abstract