Chapter 08
Geostatistics is ML you already know
Mining geostatistics was developed in the 1950s–60s by Danie Krige and Georges Matheron, decades before anyone said "Gaussian process". It is the same mathematics, in different notation, with one genuinely important idea that the ML literature underweights. This chapter gives you the translation and then the idea.
The variogram is a covariance kernel
Geostatistics describes spatial continuity with the semivariogram: the expected squared difference between two values separated by a lag vector h.
γ(h) = ½ · E[ ( Z(x + h) − Z(x) )² ]
Under second-order stationarity — constant mean, covariance depending only on the lag — expand the square and the connection falls out immediately:
γ(h) = ½ E[Z(x+h)²] + ½ E[Z(x)²] − E[Z(x+h)Z(x)]
= C(0) − C(h)
So the variogram is just the covariance function turned upside down: C(h) = C(0) − γ(h). Everything a kernel does, a variogram does.
A fitted variogram model has three parameters, and each is a hyperparameter you already know:
| Variogram parameter | GP equivalent | Meaning |
|---|---|---|
| Nugget c₀ | σ²noise | Variance at zero lag — measurement error plus sub-sample-scale variability. Gold deposits have famously huge nuggets. |
| Sill c | signal variance | The plateau the variogram reaches; total spatial variance. |
| Range a | lengthscale | Lag beyond which samples are effectively uncorrelated. |
| Anisotropy ellipsoid | ARD lengthscales / Mahalanobis metric | Continuity is longer along bedding or a fault than across it. Almost always present. |
The spherical model is the field's workhorse and has a property Matérn/RBF kernels lack — compact support, reaching the sill exactly at range a:
γ(h) = c₀ + c · [ 1.5 (h/a) − 0.5 (h/a)³ ] for 0 < h ≤ a γ(h) = c₀ + c for h > a
Kriging is GP regression
Simple kriging — the case with known mean m — is, term for term, the Gaussian process posterior mean:
Ẑ(x₀) = m + k(x₀, X) · [ K(X, X) + σ²n I ]⁻¹ · ( z − m )
Var[Ẑ(x₀)] = k(x₀, x₀) − k(x₀, X) · [ K(X, X) + σ²n I ]⁻¹ · k(X, x₀)
Identical to the standard GP predictive equations. Geostatistics derives it as the best linear unbiased estimator (BLUE) minimising estimation variance; ML derives it as a posterior under a Gaussian prior. Same equations, same answer.
| Kriging flavour | GP equivalent |
|---|---|
| Simple kriging | GP with known, fixed mean |
| Ordinary kriging | GP with unknown constant mean (weights constrained to sum to 1 via a Lagrange multiplier — equivalently, a vague prior on the mean, marginalised) |
| Universal kriging / kriging with a trend | GP with a parametric mean function (linear/polynomial basis) |
| Co-kriging | Multi-output GP with a cross-covariance structure |
| Indicator kriging | Estimating P(Z ≤ z) by regressing thresholded indicators — a crude quantile regression |
- Fit hyperparameters by marginal likelihood instead of eyeballing a curve through a noisy experimental variogram — which is still, genuinely, how much of the industry does it.
- Non-stationary and deep kernels for fields whose continuity changes across the deposit, which is the norm, not the exception.
- Scalability: inducing points, structured kernel interpolation and conjugate-gradient solvers, instead of the moving-neighbourhood heuristic (only use the nearest ~30 samples) that classical kriging uses to dodge the O(n³) solve.
- Principled multi-output modelling for the many correlated elements you assay anyway.
Look at the posterior variance equation: the data values z do not appear in it. It depends only on sample geometry and the covariance model. So kriging reports identical uncertainty for a quiet, homogeneous zone and a wildly erratic one with the same drill spacing. This is an immediate consequence of the homoscedastic Gaussian likelihood, obvious in the GP framing and a well-known trap in the geostatistical one. If you need honest uncertainty, you need simulation — next section — or a heteroscedastic / warped model.
The idea ML underweights: smoothing and change of support
This is the important part of the chapter. Kriging minimises estimation variance, so it is a conditional expectation, so it is smooth: the distribution of kriged estimates has less variance than reality. Peaks are pulled down, troughs pulled up. For a "best estimate" map that is exactly right and exactly what you asked for.
Now recall from Chapter 4 that mining reports tonnes and grade above a cutoff. That is a nonlinear functional of the field, and for nonlinear g:
E[ g(Z) ] ≠ g( E[Z] )
Jensen's inequality, in the place it costs the most money. Computing "tonnes above 0.5% Cu" from a kriged block model systematically misstates the answer — and the direction and size of the error depend on the cutoff relative to the mean.
Take a smoothed model, threshold it at a cutoff, and sum the tonnes. Because the smoothed distribution is too narrow, you under-represent both the high-grade tail and the low-grade tail. Near typical cutoffs this generally overstates tonnage at an understated grade. Mines have been financed on this error. It is the single most consequential statistical mistake in the industry, and it follows directly from using a conditional mean where you needed a distribution.
Sequential Gaussian simulation = posterior sampling
The fix is to stop using the mean and start using samples. Sequential Gaussian simulation (SGS) draws realisations from the posterior, and the algorithm is exactly ancestral sampling from the conditional Gaussian:
normal-score transform the data # make the marginal Gaussian
visit every grid node in random order:
krige mean and variance from data + previously simulated nodes
draw a value ~ N(kriged_mean, kriged_variance)
ADD that value to the conditioning set # <- this is the key step
back-transform to grade unitsEach realisation has the right variance and the right spatial continuity — it looks like plausible geology rather than a blurred blob. Run a few hundred, apply your nonlinear functional to each, and you get a distribution of tonnage and grade above cutoff. Now you can say "there is an 80% chance this block of ground contains more than X tonnes", which is the sentence a mine planner and a financier both actually need.
SGS is Monte Carlo over the posterior, done sequentially so each draw conditions on the previous ones — a Gibbs-flavoured ancestral sampler exploiting the fact that Gaussian conditionals are closed-form. The normal-score transform is a marginal warping — a Gaussian copula, or a fixed-form warped GP. If you have ever sampled functions from a GP posterior instead of plotting the mean, you have done SGS.
Support, formally
A mining block is 10 × 10 × 5 m; a sample is a few kilograms. You are not predicting the field at a point but its average over a volume:
ZV(x) = (1/|V|) ∫V Z(x + u) du
Averaging reduces variance (volume-variance relation) and pulls the distribution toward its mean. Block kriging predicts this integral directly — a GP prediction under a linear functional of the field, which the GP formalism handles cleanly because linear functionals of GPs are Gaussian.
The practical rule: never compare or combine statistics computed on different supports without an explicit change-of-support step. Chapter 3's warning about validating rasters against assays is the same theorem in a different costume.
Two more translations worth having
Declustering = importance weighting
Drilling is dense where grade is high, because that is where people kept drilling. The naive sample mean of a drillhole database is therefore biased upward — it is a sample from p(x | someone drilled here), not from the deposit. Declustering assigns each sample a weight inversely proportional to its local sample density (cell or polygonal methods). It is importance weighting to correct covariate shift, and it is the same disease as the sampling bias in Chapter 9, at a different scale.
Domaining = stratification
Stationarity almost never holds across a whole deposit — an oxide cap and a fresh sulphide zone are different populations with different means, variances and continuity. Geologists partition the deposit into domains and model each separately, with hard boundaries preventing samples from one domain influencing estimates in another. This is stratification, or a mixture of experts with a hand-drawn gating function. The interesting open problem: those boundaries are themselves uncertain, and almost nobody propagates that uncertainty through to the resource estimate.
Where the ML framing stops being enough
A caution before you go and replace someone's kriging with a neural network. Resource estimates are regulated disclosures signed by a Qualified Person (Ch. 4). The method must be explicable, auditable, reproducible, and defensible to a competent peer years later. A gradient-boosted grade model with better hold-out error and no interpretable spatial-continuity model is, for that purpose, often worse than kriging even when it is more accurate.
The productive move is not replacement. It is to use the GP framing for what it is genuinely better at — learning hyperparameters honestly, handling non-stationarity, multi-output modelling, propagating uncertainty into the decision layer of Chapter 10 — while keeping the output in a form a QP can sign.
Sources for this chapter
- Matheron, "Principles of geostatistics", Economic Geology 58 (1963), building on Krige (1951). The origin of the field, in a mining journal, decades before the ML literature arrived at the same equations. Worth reading once for the history alone.
- Rasmussen & Williams, Gaussian Processes for Machine Learning (MIT Press, 2006) — free PDF. Read chapters 2 and 4 next to this chapter and the correspondence is immediate.
- Williams, "Prediction with Gaussian processes: from linear regression to linear prediction and beyond" (1998), and Stein, Interpolation of Spatial Data (Springer, 1999). Both state the kriging↔GP equivalence explicitly, from the two different directions.
- Chilès & Delfiner, Geostatistics: Modeling Spatial Uncertainty (Wiley, 2nd ed. 2012) — the rigorous reference; and Goovaerts, Geostatistics for Natural Resources Evaluation (OUP, 1997) — the readable one, excellent on simulation and indicator methods.
- Deutsch & Journel, GSLIB (OUP, 2nd ed. 1998). Still the clearest specification of the sequential Gaussian simulation algorithm.
- Rossi & Deutsch, Mineral Resource Estimation (Springer, 2014). Change of support, recoverable resources, and why the smoothing error in this chapter is the one that costs real money.