Accurately characterising datasets is essential for effective statistical modelling, particularly when dealing with asymmetric and heavy-tailed data samples. Global Navigation Satellite System (GNSS) latitude data are examples of such samples [1]. GNSS is a satellite constellation that provides signals from space, transmitting positioning and timing data to receivers on Earth [2]. This system includes GPS (USA), GLONASS (Russia), Galileo (Europe), and BeiDou (China) [3]. It is essential for applications such as navigation, surveying, mapping, and scientific research, where precise location data are required. However, GNSS positioning errors arise from multiple sources, including ionospheric and tropospheric delays, satellite clock drift, multipath interference, and receiver noise [4]. To account for these errors, classical statistical models often assume a Gaussian distribution of positioning errors. However, empirical studies have revealed deviations from normality, suggesting that errors exhibit heavy-tailed and skewed characteristics [5]. These deviations challenge conventional error modelling approaches and highlight the need for alternative probability distributions that better capture the statistical properties of positioning errors.
Residual pseudorange errors in single-frequency differential GNSS systems arise from sources such as signal multipath, increased receiver noise due to interference, residual differential troposphere error, and ionospheric gradients [6]. To address these errors, the study proposed modelling the total pseudo-range error distribution by combining the probability density functions of individual error sources, moving beyond classical Gaussian assumptions. As GNSS errors originate from various sources, modelling their total impact requires alternative statistical approaches. The application of non-classical error theory of measurements (NETM) was proposed to better handle these deviations, particularly for large datasets (n > 500) [7]. In urban environments, GNSS positioning performance is affected by outlier measurements, including multipath effects and non-line-of-sight (NLOS) receptions.
These outliers present a significant challenge in GNSS modelling. To address these effects, Wen et al. [8] proposed a graduated non-convexity factor graph optimisation technique, which effectively reduced the impact of outliers and improved positioning accuracy. In addition to outlier mitigation, enhancing state estimation under non-Gaussian noise conditions has also been a key research focus. Raitoharju et al. [9] introduced an algorithm that integrates empirical noise models into Kalman filtering, enabling improved estimation accuracy in scenarios where Gaussian assumptions fail. This approach improved estimation accuracy by incorporating non-Gaussian noise characteristics, making it particularly useful when measurement errors deviate from Gaussian assumptions.
Although methods such as Kalman filtering can improve state estimation, selecting an appropriate statistical distribution remains a critical challenge. Several studies have examined the suitability of alternative distributions, including Cauchy, Student's t, Laplace, lognormal, and skew-normal, for data fitting and modelling applications. For example, Sun-Yong Choi et al. [10] applied twelve different distributions, including Cauchy, Laplace, normal, and Student's t, to model stock index returns, using information criteria and goodness-of-fit tests to identify the best-fit distribution. Alzaatreh et al. [11] proposes a family of generalised Cauchy distributions, denoted as T-Cauchy(Y), using the T-(RY) framework. This family is generated using the quantile functions of uniform, exponential, log-logistic, logistic, extreme value, and Fréchet distributions. It provides a flexible approach for modelling heavy-tailed data, demonstrating greater adaptability than classical Gaussian models. Aryal et al. [12] investigated the Laplace probability distribution and its derivatives, focusing on their applicability in modelling real-world problems. The study highlights the growing interest in constructing flexible parametric classes of probability distributions and examines how the Laplace distribution and its variants have been widely applied in statistical modelling.
Building on these findings, this study systematically evaluates the suitability of non-Gaussian distributions, including the Laplace, skew-normal, skew-t, and generalised hyperbolic (GH) distributions, for modelling GNSS latitude data from a stationary receiver. These models are applied using weighted maximum likelihood estimation (WMLE) to achieve robust parameter estimation and confidence interval (CI) construction.
The novelty of this work lies in developing a unified statistical framework that combines empirical GNSS data with non-Gaussian modelling and WMLE-based estimation to quantify measurement uncertainty more accurately under real-world, heavy-tailed conditions. This approach moves beyond the classical Gaussian assumption, providing a comparative evaluation of flexible asymmetric distributions using multiple goodness-of-fit criteria such as log-likelihood, Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), and root mean squared error (RMSE).
The added value of this study lies in its direct application of metrological principles from the Guide to the Expression of Uncertainty in Measurement (GUM) and the International Vocabulary of Metrology (VIM) to non-Gaussian error modelling. This enhances the statistical reliability and interpretability of GNSS uncertainty assessment, bridging the gap between theoretical distribution fitting and practical measurement evaluation.
To achieve this, we analyse empirical GNSS latitude datasets and employ rigorous statistical methodologies for model evaluation. The rest of this paper is structured as follows: Section 2 describes the methods and materials used in our research. Section 3 details the selected synthetic data distribution; Section 4 presents goodness-of-fit tests to compare the selected models; Section 5 presents the results and discussion; and Section 6 provides the conclusion.
The Laplace distribution is a continuous probability distribution characterised by a sharp peak at the mean and heavier tails than the normal distribution [13]. It can be derived as the difference between two independent exponential distributions:
The probability density function (PDF) of the Laplace distribution is [14]
The skew-normal distribution extends the normal distribution by introducing a shape parameter (α) to control skewness. It is defined by three parameters: location (μ0), which controls the centre of distributions; scale (σ), which determines the spread; and shape (α) [16]. This distribution is constructed by multiplying the normal PDF by a transformation of its cumulative distribution function (CDF). A random variable X follows a skew-normal distribution if X = μ0 + σZ, where Z follows a standard skew-normal distribution. Thus, the skew-normal PDF is given by
Substituting (5) into (2), the general skew-normal PDF is
It is a flexible class of continuous probability distributions capable of modelling heavy tails and skewness. It includes many other distributions as special cases, such as the normal, Laplace, hyperbolic, variance gamma, and normal-inverse Gaussian (NIG) distributions [17]. The key distinction between the GH distribution and the standard hyperbolic distribution is the inclusion of the index parameters α, resulting in a five-parameter family (α, λ, β, σ, and μ). The PDF of the GH distribution arises from the expression of X as a normal mean-variance mixture, where the mixing distribution is a generalised inverse Gaussian (GIG)
Overall, the GH has the following PDF after substituting (11) into (10)
It is an extension of the Student's t-distribution that allows for skewness as well as heavy tails. The Student's t-distribution models the sample mean of a normally distributed population when the sample size is small and the standard deviation is unknown [19]. The skew-t distribution is useful for modelling asymmetric data with outliers or heavy-tailed behaviour [20]. Its distribution arises from a normal mean-variance mixture. The random variables Y and X are then defined as follows
The Laplace, skew-normal, skew-t, and GH distributions were selected for their complementary ability to represent key non-Gaussian characteristics commonly observed in GNSS data, namely asymmetry, heavy tails, and excess kurtosis.
The Laplace model serves as a simple reference with moderate tail behaviour, while the skew-normal introduces a shape parameter to accommodate data asymmetry. The skew-t extends this flexibility by modelling both skewness and pronounced tail heaviness. The GH distribution provides the most general framework, encompassing several of these models as special cases and offering superior adaptability for complex GNSS error structures.
Together, these distributions form a progressive hierarchy from simple to highly flexible models, enabling a systematic evaluation of non-Gaussian behaviour in GNSS measurements. Broader generalisations, such as the generalised alpha-skew-t, variance gamma, or NIG distributions, could further enhance modelling accuracy by offering additional control over tail behaviour and asymmetry, and are recommended for future investigation.
A goodness-of-fit test is a statistical method used to determine how well a set of observed data matches a specific theoretical distribution or model [21].
To evaluate the suitability of various statistical models for GNSS data, it is essential to consider distributions capable of capturing asymmetry and heavy tails. Among the candidate models are the Laplace, skew-normal, GH, and skew-t distributions, each offering different levels of flexibility in representing real-world data characteristics. Parameter estimation for these distributions is typically performed using the WMLE method. This is a modified form of maximum likelihood estimation (MLE) that applies weights to reduce bias, particularly for individuals with extreme response patterns [22]. The estimated parameter vector θ* is defined as
θ is the vector of parameters to be estimated,
xi is the i-th data point,
wi is the weight assigned to xi, satisfying
\Sigma_{i = 1}^n{w_i} = 1 ℒ(θ | xi) is the likelihood of observing xi given parameters of θ.
This approach provides robust parameter estimates and enables the construction of 95 % CIs for the fitted models. CIs are determined using the percentile method, which relies on the quantiles of the fitted distribution. The general formula for the CI is given by [23]
In this study, the constructed CIs correspond to the expected range of the measured quantity (GNSS latitude) as represented by each fitted distribution. Thus, the intervals indicate the range within which the true value of the measurand is expected to lie with a 95 % confidence level.
To evaluate the adequacy of statistical models in representing GNSS error distributions, various goodness-of-fit metrics are employed. These include the log-likelihood, AIC, BIC, and RMSE. Together, these measures provide a rigorous framework for identifying the most appropriate distributional model.
The log-likelihood function is a tool for evaluating how well a probability distribution fits the observed data. It quantifies the data compatibility with the assumed model by summing the logarithms of the PDF values at each observation. Higher values indicate a better fit.
For n independent observations X = {x1, x2, . . . ,xn} and a PDF f (x | θ), the log-likelihood is defined as [24]
The log-likelihood function of the skew-normal distribution is
After substituting (6) into (18), the log-likelihood function is
By substituting (12) into (18), the log-likelihood function of the GH distribution is as follows
Substituting the skew-t PDF (15) into (18), the for log-likelihood is
The AIC is a widely used metric for model selection, balancing goodness-of-fit with model complexity and fit quality. It is defined as [25]:
The BIC is a model selection tool that imposes a stronger penalty on model complexity than the AIC. The BIC is defined as [26]
The penalty term k ln(n) ensures that the BIC strongly discourages overfitting, especially for large datasets.
The Laplace distribution is characterised by two parameters (k = 2). Based on (25), the corresponding BIC values are calculated using the formula: BIC = −2ℒ(θ) + 2 · ln(n).
Similarly, the skew-normal and skew-t distributions each involve three parameters (k = 3). According to (25), their BIC values are computed as: BIC = −2ℒ(θ) + 3 · ln(n).
The GH distribution involves five parameters (k = 5) and BIC is expressed as: BIC = −2ℒ(θ) + 5 · ln(n). A lower BIC indicates a better model fit while imposing a stronger penalty on model complexity than the AIC, thereby helping to prevent overfitting.
RMSE is a widely used metric for quantifying the accuracy of a model or estimator by measuring the average deviation between observed and predicted values. By squaring the differences, RMSE penalises larger errors more heavily, providing a sensitive measure of overall fit [27]. It is particularly useful for comparing multiple models or candidate distributions in terms of their ability to represent the observed data. In our GPS latitude data analysis, we use RMSE to evaluate how closely each fitted probability distribution aligns with the observed measurements, thereby establishing a quantitative basis for comparing model performance. Formally, for a set of n observed values yi and corresponding fitted values ŷi derived from a candidate distribution, RMSE is defined as
The RMSE for the Laplace distribution with mean μ0 and variance 2σ2 can be directly computed as
For the skew-normal distribution with mean
For the GH distribution with parameters α, λ, β, σ, and μ, the mean and variance are
Then, the RMSE is
For a skew-t distribution with location parameter μ, scale σ , skewness λ , and degrees of freedom ν, the mean and variance (for ν > 2 and ν > 4, respectively) are given by
The presented formulas [(28), (29), (32), and (35)] explicitly define how the RMSE values were computed for each fitted distribution. Each expression corresponds to the RMSE calculation applied to a specific model, quantifying the average squared deviation between the observed GPS latitude data and the values predicted by the fitted distribution. These computed RMSE values are reported in Table 1, allowing direct comparison of the fitting accuracy among the candidate distributions.
Estimated parameters with 95 % CIs and performance metrics for the fitted models applied to GNSS latitude data.
| Model | Parameters | CI* | ℒ | AIC | BIC | RMSE | |||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| μ* | σ | ν | α | λ | β | ||||||
| ×10−5 | ×10−5 | ×103 | ×106 | ×106 | ×106 | ||||||
| Laplace | 864 | 20 | – | – | – | – | [805; 924] | 1.16 | −2.33 | −2.33 | 3.5 · 10−4 |
| Skew-normal | 857 | 21 | – | 0.50 | – | – | [825; 904] | 1.13 | −2.27 | −2.27 | 2.5 · 10−4 |
| GH | 845 | 0.06 | – | 0.05 | 1.39 | 0.43 | [841; 917] | 1.20 | −2.40 | −2.40 | 2.8 · 10−4 |
| Skew-t | 861 | 14 | 5.0 | – | – | – | [825; 897] | 1.16 | −2.32 | −2.32 | 2.6 · 10−4 |
Note: The values marked with an asterisk (μ and CI) are expressed as least significant digits relative to a base latitude of 50.28000° N and scaled by 10−5.
To collect GNSS sample data, we deployed an L76k GPS module at a fixed outdoor location. This module, known for its compact design and high positioning accuracy, was connected to a Raspberry Pi, which served as the control system. The Raspberry Pi handled module configuration, data logging, and real-time monitoring, ensuring continuous and reliable data collection throughout the experiment.
During 72 hours of continuous operation, the GPS module recorded 1 303 376 rows of National Marine Electronics Association (NMEA) sentences. These included various sentence types such as:
$GPRMC:Recommended Minimum Specific GNSS data,$GPGGA:Global Positioning System fix data,$GPGLL:Geographic Position Latitude / Longitude,$GPVTG:Track made good and ground speed,$GPGSA:GNSS DOP and Active Satellites,$GPGSV:GNSS Satellites in View.
From this dataset, 162 922 $GPRMC sentences were extracted for further analysis. These sentences were selected because they include essential positional and signal quality parameters such as fix quality, satellite count, horizontal dilution of precision (HDOP), and latitude / longitude coordinates. Python scripts were used to efficiently filter, parse, and structure the extracted data into a log file.
The GNSS latitude data distribution was first characterised before fitting multiple non-Gaussian distributions (Laplace, skew-normal, skew-t, and GH). Parameter estimation and CIs for each distribution were determined using WMLE to evaluate measurement reliability and quantify uncertainty. In addition, model performance was assessed using log-likelihood, AIC, BIC, and RMSE.
Fig. 1 presents the PDF and CDF for the latitude data, compared with four fitted distributions using WMLE: Laplace, skew-normal, GH, and skew-t. For numerical consistency, the latitude values were normalised by subtracting the reference latitude of 50.28000° N, so the values reported in Table 1 correspond to the least significant digits of the measured coordinates. For example, the Laplace model mean value 864 × 10−5 represents an actual latitude of approximately 50.28864° N, with a 95 % CI of [50.28805°, 50.28924°]. The same interpretation applies to the other models listed in Table 1. In Fig. 1(a), the empirical PDF (blue histogram) exhibits a right-skewed and heavy-tailed structure. The Laplace distribution fits the peak well but fails to capture the tail behaviour, while the skew-normal improves on the asymmetry but still under-represents the tail. In contrast, the skew-t distribution closely matches the empirical density across the entire range, accurately modelling the sharp peak and heavy right tail. The GH distribution exhibits the best alignment with the empirical data, particularly capturing the skewness and heavy tail behaviour, which aligns with the goodness-of-fit metrics in Table 1. Fig. 1(b) shows the empirical CDF along with the fitted CDFs. The Laplace distribution deviates significantly around the central region, while the skew-normal offers better central alignment but lacks precision in the tails. GH again provides the closest match to the empirical CDF, with skew-t exhibiting slightly superior tail conformity. Overall, both PDF and CDF analyses confirm that the GH distribution offers the best fit for the observed GNSS latitude data, effectively capturing its skewed and heavy-tailed nature.
Probability density function (a) and cumulative distribution function (b) for latitude data with fitted distributions using weighted MLE.
Table 1 presents the estimated parameters, CIs, and performance metrics for four fitted models applied to GPS latitude data. Among these, the GH distribution demonstrates the best overall performance, with the highest log-likelihood (ℒ = 1.20 × 106) and the lowest AIC and BIC values (−2.40 × 106) for both, indicating an optimal balance between goodness-of-fit and model complexity. It also achieves a low RMSE (2.8 × 10−4), confirming its accuracy in approximating the observed data. The skew-t model also performs competitively, with a slightly lower RMSE (2.6 × 10−4) and a similar log-likelihood, although its AIC and BIC values are marginally less favourable than those of the GH. The GH distribution incorporates additional parameters, α, λ , and β, to model heavy tails, kurtosis, and skewness, respectively. Notably, the small value of α = 0.05 suggests mild asymmetry, highlighting the distribution's suitability for capturing the non-Gaussian characteristics of the data. The skew-normal model performs moderately well, capturing skewness through its shape parameter (α = 0.50), but its RMSE (2.5 × 10−4) and lower likelihood suggest it is less flexible in modelling heavy tails. In contrast, the Laplace model performs worse, with the highest RMSE (3.5 × 10−4) and the lowest likelihood, failing to capture the data's asymmetry and heavy tails. Overall, the GH distribution provides the most reliable representations of the latitude data based on both fit quality and parameter efficiency. Fig. 1 and Table 1 together evaluate the statistical models' fit to GPS latitude data, supporting both qualitative and quantitative analysis.
This study uses probability–probability (P–P) plots to visually assess the agreement between empirical latitude data and theoretical distributions. These plots compare the CDF of the observed data with that of a specified theoretical model by plotting empirical probabilities against theoretical ones, making them a valuable tool for evaluating goodness-of-fit [28].
Fig. 2 shows the P–P plots for the residual analysis of four fitted distributions: Laplace, skew-normal, GH, and skew-t-against the empirical CDF of the latitude data. The closer the blue curve is to the red dashed line (representing a perfect fit), the better the distribution matches the empirical data. In Fig. 2(c), the GH distribution aligns most closely with the diagonal, indicating the best overall fit and demonstrating its effectiveness in capturing both skewness and kurtosis in the data. The skew-normal and skew-t distributions in Fig. 2(b) and (d) show a generally good fit, though slight deviations at the lower tails and centre suggest minor limitations in capturing the data's asymmetry or kurtosis. In contrast, the Laplace distribution in Fig. 2(a) deviates more noticeably from the diagonal, especially in the lower and upper tails, indicating a poor fit to the empirical data. This is consistent with Laplace's limited ability to capture asymmetry and heavy tails. Overall, the GH distribution demonstrates superior performance in residual behaviour, confirming its suitability for modelling the GPS latitude data compared to the other distributions analysed.
Residual analysis of the selected distributions using P–P plot.
In summary, the P–P plots confirm that the GH distribution provides the best empirical fit, followed by the skew-t and skew-normal models, with Laplace performing the worst. These visual results are consistent with the quantitative metrics from the likelihood and RMSE analyses.
This study critically assessed four probability distributions: Laplace, skew-normal, GH, and skew-t-for modelling skewed and heavy-tailed data using GNSS latitude data as an example. Using the WMLE, we estimated parameters and CIs to quantify uncertainty. We also evaluated model performance using log-likelihood, AIC, BIC, and RMSE.
The analysis showed that, although the Laplace model offers a simple symmetric representation, it performed the weakest overall. Its limited ability to capture heavy tails and asymmetry resulted in the highest RMSE, lowest likelihood, and the widest CIs. The skew-normal model improved fit by capturing asymmetry and producing tighter intervals, yet it struggled with heavy tails.
The skew-t distribution provided a more balanced fit, effectively modelling both skewness and moderate tail behaviour. It performed well across all metrics, with only minor deviations in the lower tails and centre. The GH distribution proved most robust, closely aligning with empirical data and excelling in both statistical and visual evaluations. Its additional shape, skewness, and kurtosis parameters enabled it to capture complex distributional features, outperforming all other models. Both GH and skew-t showed tighter CIs, indicating more reliable parameter estimates.
The novelty of this work lies in establishing a unified statistical framework that integrates WMLE-based parameter estimation with multiple non-Gaussian models to evaluate GNSS measurement uncertainty under real-world, non-Gaussian noise conditions. Unlike previous research, which focused mainly on Gaussian assumptions, this study demonstrates that heavy-tailed and asymmetric models, particularly GH and skew-t, significantly improve distributional fit and the accuracy of uncertainty quantification.
These results highlight the importance of flexible, asymmetric models for accurate GNSS error characterisation. While RMSE quantifies overall deviation, it overlooks distributional nuances evident in models such as skew-normal, which achieve low RMSE but misrepresent tail structure. In contrast, GH and skew-t models deliver superior performance across all criteria. AIC and BIC further support their suitability by balancing fit and complexity, ensuring robustness and generalisability to new datasets.
This work has practical implications for GNSS-based positioning, calibration, and metrological applications, particularly in environments where impulsive or non-Gaussian noise dominates. Future work will extend this framework to dynamic GNSS scenarios and explore its integration with advanced filtering and real-time uncertainty estimation techniques.