Most treatments of Gaussian process regression begin with a comfortable assumption: you know where your data came from. The spatial coordinates of your observations are fixed, known quantities, and the inferential challenge is purely about learning the latent field. A recent post by Christopher Krapu dismantles this assumption in a principled way, constructing a full Bayesian model in which the spatial coordinates themselves are latent variables subject to measurement noise. The application domain is mineral exploration, but the modelling ideas generalise broadly to any setting where georeferencing error is non-trivial.
The Problem: Coordinate Uncertainty in Spatial Models
The motivating scenario is geologic sampling during mineral exploration. Drill cores are extracted, ore concentrations are measured, and the resulting dataset is spatially referenced. The trouble is that in field conditions, the recorded coordinates may be substantially wrong. GPS drift, transcription errors, or the physical displacement between a drill collar and the actual subsurface sample location can all introduce meaningful positional uncertainty. The dataset used here consists of uranium and vanadium concentrations from Walker Lake, drawn from Isaaks and Srivastava's classic geostatistics text and distributed via the R package gstat.
The standard GP regression framework cannot handle this gracefully. If you simply treat noisy coordinates as ground truth, your covariance matrix is computed at the wrong locations, and the resulting predictions will be miscalibrated in ways that are difficult to diagnose. The geostatistical literature has acknowledged this problem, with Cressie and Kornak and Cervone and Pillai both emphasising that location error changes the covariance structure in a fundamental way, not merely adding noise to the response.
The Hierarchical Model Structure
The proposed model is a hierarchical extension of the standard marginal GP. Each observation has a recorded coordinate s̃i and a latent true coordinate si, related by an additive Gaussian displacement:
- Location displacement: Δi ~ Normal(0, σs² I2), so si = s̃i + Δi
- Global mean: μ ~ Normal(2, 2)
- GP marginal variance: σ ~ HalfNormal(1)
- Length scale: ℓ ~ HalfNormal(100)
- Observation noise: σ0 ~ HalfNormal(0.5)
- Latent field: f(·) | σ, ℓ ~ GP(0, σ² c(·, ·; ℓ)) with a Matern-5/2 kernel
- Likelihood: Yi | f, μ, σ0, si ~ Normal(μ + f(si), σ0)
The GP is evaluated at the latent coordinates, not the observed ones. The scale σs is treated as known and varied across experiments to simulate different levels of coordinate degradation. The GP values at observation locations are marginalised out using pm.gp.Marginal, which is the standard Cholesky-based approach that avoids explicitly sampling the latent field at training points.
This is a clean and well-motivated construction. The prior on Δi is isotropic Gaussian in Cartesian coordinates, which is a natural default when you have no directional information about the error source. The author notes that polar coordinates with priors on magnitude and angle would be an equally valid parameterisation, which is correct and worth considering when the error has a known directional bias, for instance from systematic GPS drift along a survey transect.
Computational Challenges and Sampling Diagnostics
The computational cost here is substantially higher than for a fixed-location GP. In the standard case, the covariance matrix is computed once and factored once per likelihood evaluation. Here, because the latent coordinates change with every NUTS proposal, the covariance matrix must be recomputed and re-factored at every step. For a dataset of size n, this is O(n³) per sample, which is the dominant cost. The sampling logs confirm this: each of the three experiments (at σs = 12, 25, and 40 metres) took roughly 6,000 to 7,700 seconds of wall-clock time with only two chains and 1,000 draws each.
The diagnostics are concerning in the expected ways:
- Divergences ranged from 24 at σs = 12 to 79 at σs = 25, indicating a posterior geometry that NUTS struggles to navigate
- R-hat statistics reached 1.03 at the intermediate noise level, suggesting incomplete chain mixing
- Minimum bulk ESS fell to 44 at σs = 40, which is well below the threshold needed for reliable inference
- Mean displacement of the inferred Δi scaled roughly linearly with σs, as expected
These are not flaws in the modelling approach; they are honest reflections of a genuinely hard posterior. The coordinate displacements Δi are weakly identified when σs is large relative to the GP length scale, because many different coordinate configurations can produce similar marginal likelihoods. The posterior over locations becomes diffuse, and the coupling between location and field value creates a highly correlated, funnel-shaped geometry that is notoriously difficult for HMC-based samplers. A non-centred parameterisation of the displacements, or a variational approximation as a warm-start, would likely help substantially.
Results and Comparison with Nadaraya-Watson Smoothing
The qualitative results are striking. Even at σs = 40 metres, which represents severe positional degradation, the posterior mean surface retains recognisable spatial structure: a lighter region in the lower left and darker regions in the upper left and lower right. This is compared against a Nadaraya-Watson kernel smoother with bandwidth set to the posterior mean length scale estimate, which produces little more than a spatially uniform average with negligible local variation.
The second row of the main figure is particularly informative. It shows posterior density contours for the inferred true locations of four selected points. As σs increases, these contours expand gracefully, correctly reflecting growing uncertainty about where each measurement actually occurred. This is exactly the behaviour you want from a well-calibrated uncertainty model: the posterior does not collapse to a point estimate at the recorded coordinate, but instead spreads across the plausible region consistent with the known error scale.
The comparison with Nadaraya-Watson is somewhat unfair to the kernel smoother, which was not designed for this problem, but it makes a useful pedagogical point. A naive interpolation method has no mechanism for representing or propagating coordinate uncertainty. It will smooth the data as if the recorded locations were exact, producing predictions that are confidently wrong in a way that the Bayesian model avoids.
Broader Implications and Open Directions
The methodological contribution here is less about the specific application and more about demonstrating that Bayesian hierarchical modelling can absorb measurement error in the covariates themselves, not just in the response. This connects to a broader literature on errors-in-variables models, Berkson error, and covariate measurement error in regression, but the spatial GP setting introduces the additional complication that covariate error directly modifies the covariance structure rather than simply attenuating a regression coefficient.
Several extensions would be worth pursuing. First, treating σs as an unknown parameter rather than a fixed hyperparameter would allow the model to learn the scale of coordinate error from the data itself, though this would require informative priors to avoid non-identifiability. Second, the isotropic Gaussian error model is a strong assumption; in practice, GPS errors often have anisotropic structure correlated with satellite geometry or terrain. Third, scalability is a genuine concern. Sparse GP approximations such as inducing-point methods or Hilbert space GPs could reduce the O(n³) bottleneck, though care would be needed to ensure the approximation remains valid when training locations are themselves uncertain.
The post is a good example of what principled probabilistic modelling buys you: the ability to encode domain knowledge about error sources directly into the model structure, then let the inference machinery propagate that uncertainty through to predictions. The sampling difficulties are real and should not be minimised, but they are a consequence of the problem's inherent difficulty rather than a modelling mistake. As variational inference methods and better HMC geometry adaptation mature, this class of models will become more tractable. For now, the approach is computationally expensive but epistemically honest, which is often the right trade-off when the stakes of miscalibrated predictions are high.