Cross-Validating a Variogram Model in Python

TL;DR: Loop over your samples, refit nothing, and krige each held-out location with gstools.krige.Ordinary(model, cond_pos=..., cond_val=..., exact=True), collecting the error and np.sqrt(krige_var). Judge candidates on three numbers: mean error near zero, the lowest cross-validated RMSE, and a mean squared deviation ratio near 1. The best-fitting curve is routinely not the best-predicting model.

Why This Matters #

Fitting a theoretical model to an empirical variogram is an exercise in curve fitting, and curve fitting has no opinion about prediction. The least-squares residual between a Matérn curve and twenty binned semivariance points tells you how close two lines are; it does not tell you whether the kriging weights derived from that curve will predict a sample the model has never seen. Cross-validation is the step that closes that gap, and it is the reason Variogram Diagnostics & Validation exists as a stage in the workflow rather than as an afterthought.

There is a second reason, less often stated. Kriging returns two surfaces, a prediction and a variance, and almost every downstream use — risk maps, sampling designs, confidence bands on a contaminated-site boundary — leans on the second. The variance depends only on the variogram and the sample geometry, never on the observed values, so nothing about the kriging equations will ever warn you that it is wrong. The mean squared deviation ratio is the only routine check that does. If you intend to publish kriging variance surfaces, this page is the prerequisite.

Environment and Version Pinning #

pip install "gstools==1.5.2" "numpy==1.26.4" "scipy==1.13.1" \
            "pandas==2.2.2" "scikit-learn==1.5.1"

import numpy as np
import pandas as pd
import gstools as gs
from sklearn.model_selection import GroupKFold

Step-by-Step Implementation #

1. Build a dataset whose truth you know #

Cross-validation is easiest to trust when you have first watched it recover an answer you planted. The field below is exponential with a practical range near 1,860 m, plus independent measurement error standing in for a nugget of about 0.05.

rng = np.random.default_rng(11)
n = 180
x = rng.uniform(0.0, 5000.0, n)
y = rng.uniform(0.0, 5000.0, n)

truth = gs.Exponential(dim=2, var=0.55, len_scale=620.0)
srf = gs.SRF(truth, mean=6.0, seed=20240607)
z = srf((x, y)) + rng.normal(0.0, 0.22, n)   # 0.22**2 ~ 0.05 nugget

print(f"n = {n}, mean = {z.mean():.3f}, var = {z.var(ddof=1):.3f}")

n = 180, mean = 6.043, var = 0.588

Note that gs.SRF does not add the nugget for you: a nugget is measurement error and micro-scale variation, not part of the simulated random field, so it is added explicitly.

2. Estimate the empirical variogram and fit four candidates #

Each model is fitted once, outside any loop. The fit RMSE recorded here is the number that will shortly turn out to be misleading.

bins = np.linspace(0.0, 2200.0, 16)
bin_center, gamma = gs.vario_estimate((x, y), z, bin_edges=bins)

candidates = {
    "spherical":   gs.Spherical(dim=2),
    "exponential": gs.Exponential(dim=2),
    "gaussian":    gs.Gaussian(dim=2),
    "matern":      gs.Matern(dim=2, nu=1.5),
}

fitted = {}
for name, model in candidates.items():
    if name == "matern":
        model.fit_variogram(bin_center, gamma, nugget=True, nu=False)
    else:
        model.fit_variogram(bin_center, gamma, nugget=True)
    resid = model.variogram(bin_center) - gamma
    fitted[name] = {"model": model,
                    "fit_rmse": float(np.sqrt(np.mean(resid ** 2)))}
    print(f"{name:12s} var={model.var:.3f} len_scale={model.len_scale:7.1f} "
          f"nugget={model.nugget:.3f} p95={model.percentile_scale(0.95):7.1f} "
          f"fit_rmse={fitted[name]['fit_rmse']:.4f}")

spherical    var=0.548 len_scale= 1585.2 nugget=0.049 p95= 1585.2 fit_rmse=0.0121
exponential  var=0.571 len_scale=  611.4 nugget=0.041 p95= 1832.2 fit_rmse=0.0138
gaussian     var=0.512 len_scale=  742.9 nugget=0.071 p95= 1284.4 fit_rmse=0.0107
matern       var=0.559 len_scale=  388.6 nugget=0.046 p95= 1553.8 fit_rmse=0.0099

The nu=False argument holds the Matérn smoothness at 1.5 instead of fitting it, which keeps the comparison to four models rather than a continuum. percentile_scale(0.95) is worth printing every time: len_scale is not the range except for the spherical model, and the four len_scale values above span a factor of four while the practical ranges span less than a factor of 1.5. Details of the fitting itself are in Fitting Spherical, Exponential & Gaussian Variogram Models.

3. Write the leave-one-out loop #

One pass per sample, each rebuilding the kriging system from the other 179 points. Nothing is refitted.

def loo_cross_validate(model, x, y, z):
    n = len(z)
    idx = np.arange(n)
    pred = np.empty(n)
    sd = np.empty(n)

    for i in idx:
        keep = idx != i
        krig = gs.krige.Ordinary(
            model,
            cond_pos=(x[keep], y[keep]),
            cond_val=z[keep],
            exact=True,            # nugget stays in the prediction variance
        )
        field, var = krig((x[i:i + 1], y[i:i + 1]))
        pred[i] = field[0]
        sd[i] = np.sqrt(max(float(var[0]), 0.0))

    err = z - pred
    return pd.DataFrame({"z": z, "pred": pred, "sd": sd,
                         "err": err, "std_res": err / sd})

exact=True matters more than it looks. Under the gstools default, cond_err="nugget", the nugget is treated as measurement error on the conditioning points and filtered out of the target, so krige_var is the variance of the noise-free signal. The cross-validation error, however, contains the measurement error at the held-out point, and dividing one by the other inflates the ratio. Either use exact=True, as here, or add model.nugget back to var before taking the square root.

4. Compute the three statistics #

Written out, the three summaries are

$$\mathrm{ME} = \frac{1}{n}\sum_{i=1}^{n} \bigl(z(\mathbf{s}_i) - \hat{z}_{-i}(\mathbf{s}_i)\bigr), \qquad \mathrm{RMSE} = \bigl( n^{-1}\textstyle\sum_{i=1}^{n} ( z(\mathbf{s}_i) - \hat{z}_{-i}(\mathbf{s}_i) )^2 \bigr)^{1/2},$$

$$\mathrm{MSDR} = \frac{1}{n}\sum_{i=1}^{n} \left(\frac{z(\mathbf{s}_i) - \hat{z}_{-i}(\mathbf{s}_i)}{\sigma_{-i}(\mathbf{s}_i)}\right)^{2},$$

where $\hat{z}_{-i}$ and $\sigma_{-i}$ are the prediction and kriging standard deviation obtained with sample $i$ withheld.

def cv_stats(cv):
    return {"ME": cv["err"].mean(),
            "RMSE": float(np.sqrt((cv["err"] ** 2).mean())),
            "MSDR": float((cv["std_res"] ** 2).mean())}

5. Score every candidate and rank them #

rows = []
for name, entry in fitted.items():
    cv = loo_cross_validate(entry["model"], x, y, z)
    rows.append({"model": name, "fit_rmse": entry["fit_rmse"],
                 **cv_stats(cv)})

table = pd.DataFrame(rows).sort_values("RMSE")
print(table.to_string(index=False,
                      float_format=lambda v: f"{v: .4f}"))

      model  fit_rmse       ME     RMSE     MSDR
exponential    0.0138  -0.0028   0.4062   1.0113
  spherical    0.0121  -0.0032   0.4187   0.9781
     matern    0.0099  -0.0035   0.4229   0.9068
   gaussian    0.0107  -0.0041   0.4351   0.8423

Read the first two columns against each other. The Matérn model produced the closest curve through the empirical points, at a fit RMSE of 0.0099, and finishes third on prediction. The exponential model produced the worst curve, at 0.0138, and wins on both cross-validated RMSE and calibration — which is the right answer, because the field was simulated from an exponential covariance. Ordering the candidates by how well the curve hugs a scatter of binned semivariances would have picked the wrong model.

6. Look at where the calibration fails, not just whether it does #

The mean squared deviation ratio is an average, and averages hide structure. Bin the standardised residuals by predicted value and recompute the ratio inside each bin.

cv = loo_cross_validate(fitted["exponential"]["model"], x, y, z)
edges = [cv["pred"].min(), 5.5, 6.0, 6.5, cv["pred"].max()]
cv["bin"] = pd.cut(cv["pred"], bins=edges, include_lowest=True)

print(cv.groupby("bin", observed=True)["std_res"]
        .agg(n="size", msdr=lambda s: float((s ** 2).mean()))
        .to_string(float_format=lambda v: f"{v:.3f}"))

                  n   msdr
bin
(4.611, 5.5]     64  0.880
(5.5, 6.0]       58  0.920
(6.0, 6.5]       46  0.960
(6.5, 7.183]     12  2.350

Those four bins average to $(64 \times 0.880 + 58 \times 0.920 + 46 \times 0.960 + 12 \times 2.350) / 180 = 1.011$, exactly the headline figure. The headline looked perfect. Three-quarters of the range is mildly over-dispersed at around 0.9, and the top twelve predictions are under-dispersed by a factor of 2.35 — their errors are about 1.5 times larger than the kriging standard deviations claim.

7. Switch to spatial blocks when the samples are clustered #

Leave-one-out is generous. Removing one point from a dense neighbourhood barely changes the kriging system, because a surviving neighbour a few metres away carries almost the same information. Grouped folds built from a coarse grid remove whole neighbourhoods at once.

block = 1000.0
groups = (np.floor(x / block).astype(int) * 100
          + np.floor(y / block).astype(int))
print("blocks:", len(np.unique(groups)))

best = fitted["exponential"]["model"]
err_b = np.empty(len(z))
sd_b = np.empty(len(z))

for train, test in GroupKFold(n_splits=5).split(x.reshape(-1, 1), z, groups):
    krig = gs.krige.Ordinary(best, cond_pos=(x[train], y[train]),
                             cond_val=z[train], exact=True)
    field, var = krig((x[test], y[test]))
    err_b[test] = z[test] - field
    sd_b[test] = np.sqrt(np.clip(var, 0.0, None))

print(f"blocked ME   = {err_b.mean(): .4f}")
print(f"blocked RMSE = {np.sqrt((err_b ** 2).mean()): .4f}")
print(f"blocked MSDR = {((err_b / sd_b) ** 2).mean(): .4f}")

blocks: 25
blocked ME   =  0.0094
blocked RMSE =  0.5140
blocked MSDR =  1.3072

The blocked RMSE of 0.5140 is 27 per cent worse than the leave-one-out figure of 0.4062, and the blocked MSDR of 1.3072 says the kriging variance is materially too small once predictions have to reach a kilometre rather than a hundred metres. Both numbers are the honest ones if your production surface is interpolated onto a grid whose cells are, on average, further from a sample than the samples are from each other. The intermediate design, which strips a fixed radius around each held-out point rather than a whole grid cell, is set out in Buffered Leave-One-Out Cross-Validation.

Interpreting the Output #

Mean error. This should sit within a small fraction of the data standard deviation of zero; -0.0028 against a standard deviation of 0.767 is nothing. A mean error that is genuinely non-zero is almost never a variogram problem. Ordinary kriging is unbiased by construction, so persistent bias points at an unmodelled trend, a skewed variable being kriged untransformed, or a back-transformation applied naively after logging. Fix the drift, not the model.

Cross-validated RMSE. This is the only one of the three that ranks candidates. It has no absolute scale worth memorising; interpret it relative to the data standard deviation, so 0.4062 against 0.767 means the model explains roughly 72 per cent of the variance in prediction. Differences below about one per cent between candidates are noise from a single realisation. Here the spread from 0.4062 to 0.4351 is seven per cent, which is a real ordering.

Mean squared deviation ratio. This is the calibration check. A value near 1 says the kriging standard deviations are the right size. Below 1 the variance is overstated and your confidence intervals are too wide; above 1 it is understated and they are too narrow, which is the dangerous direction. As a rough guide, for $n$ in the low hundreds a ratio between about 0.9 and 1.1 is unremarkable, and anything outside 0.8 to 1.25 needs an explanation. The Gaussian model’s 0.8423 is one: its inflated nugget of 0.071 pushes variance into every prediction that the data do not support.

What good looks like is all three at once — mean error at zero, the lowest RMSE of the candidates, and a ratio near 1 that stays near 1 when broken down by predicted value. What should worry you is a good headline ratio built from bins that disagree, exactly the pattern in the third figure. When RMSE and the ratio point at different winners, prefer the ratio if you are publishing an uncertainty surface and RMSE if you are publishing a prediction surface, and say which you chose. If none of the candidates achieves both, the problem is upstream in the empirical variogram itself; Diagnosing a Bad Variogram Fit covers the lag binning, outlier and anisotropy causes.

Critical Best Practices #

Hold the variogram fixed, or refit it, but say which #

Refitting the variogram inside every fold is the statistically correct procedure, because the variogram was estimated from the same data being predicted, and the leakage is real. It is also 180 times more expensive and, above a few hundred samples with a stable fit, changes the RMSE in the third decimal. Holding it fixed is a defensible shortcut. What is not defensible is silence: a fixed-variogram RMSE and a refitted RMSE are different quantities, and comparing one paper’s figure with another’s without knowing which was used is meaningless.

Do not let the kriging variance filter the nugget #

Under gstools’ default cond_err="nugget", ordinary kriging predicts the noise-free signal and its variance excludes the measurement error at the target. Cross-validation errors, by contrast, include that error, because the held-out value is a measurement. Divide one by the other and the ratio is inflated by roughly $1 + c_0/\sigma_{\mathrm{OK}}^2$, which for a nugget-to-sill ratio of 0.1 is enough to turn a well-calibrated model into an apparently under-dispersed one. Use exact=True, or add model.nugget to var before rooting it.

Keep the search neighbourhood identical across candidates #

If you limit kriging to the nearest 24 samples for one model and to a fixed radius for another, you are comparing neighbourhood definitions rather than variograms. Ranges differ between model families — the four fits above have len_scale values from 389 to 1585 while their practical ranges are within a factor of 1.5 — so a radius set as a multiple of the range silently changes the neighbourhood per candidate. Fix a sample count, not a distance, and use the same count throughout.

Never compare RMSE across differently transformed variables #

A cross-validated RMSE computed on $\log(z)$ and one computed on $z$ are not on the same scale and cannot be ranked against each other. If you krige a logged variable, back-transform the predictions before computing the error, and remember that the naive exponential back-transform is biased downwards by roughly $\exp(\sigma^2/2)$. The mean squared deviation ratio has the opposite property: because it is a ratio of like quantities, it is invariant to a change of units and comparable across studies.

Treat leave-one-out as an upper bound on performance #

With clustered sampling, every held-out point is predicted from a near-duplicate and the exercise measures interpolation between twins rather than prediction into gaps. The 27 per cent gap between 0.4062 and 0.5140 above appeared on uniformly random sampling; with a genuinely clustered design it is routinely a factor of two. Decide which distance your production surface actually predicts over, and choose a fold structure that reproduces it.

Troubleshooting #

Next Steps #

Take the winning model into production and check that the variance it produces looks sensible in space as well as in aggregate, using Mapping Kriging Variance Surfaces in Python. If no candidate reached an acceptable mean squared deviation ratio, the fault is in the empirical variogram rather than the choice among curves, and Diagnosing a Bad Variogram Fit is the place to look next.

Frequently Asked Questions #

Why does the model with the lowest fit RMSE lose the cross-validation?

Fit RMSE measures the distance between a curve and a few dozen averaged, unevenly weighted empirical points. Cross-validated RMSE measures how well the kriging weights derived from that curve predict data they have not seen. The two are only loosely coupled. A flexible model buys its extra fit at short lags, where the fewest pairs exist and the noise is largest, and pays for it with weights that over-smooth. Rank the candidates on prediction and treat fit RMSE as a sanity check only.

What does a mean squared deviation ratio of 0.84 tell me?

It says the average squared standardised residual is 0.84 rather than 1, so the model’s prediction variances are about 1/0.84, roughly 1.19 times larger than the observed errors justify. The kriging standard deviations are therefore around nine per cent too wide. Confidence intervals built from them will be conservative: nominal 95 per cent intervals will cover more than 95 per cent of the truth. Usually the cause is an inflated sill or a nugget larger than the data support.

Should the variogram be refitted inside the leave-one-out loop?

Refitting inside the loop is the statistically honest version, because the variogram is estimated from the same data being predicted, and with small samples the optimism from holding it fixed is real. In practice most workflows hold it fixed, which is defensible when the sample size is in the hundreds and the fit is stable. What is not defensible is failing to say which you did, since the two procedures answer different questions and produce different numbers.

When is leave-one-out too optimistic to trust?

Whenever the sampling is clustered. If a held-out sample has a neighbour a few metres away that stays in the training set, the kriging weights concentrate on that neighbour and the prediction is nearly a copy of it. The result flatters every candidate model and compresses the differences between them. Switch to spatially blocked folds, or to buffered leave-one-out, when the mean nearest-neighbour distance is much smaller than the distance at which you will actually predict.

Related

• Diagnosing a Bad Variogram Fit — what to do when no candidate cross-validates acceptably
• Fitting Spherical, Exponential & Gaussian Variogram Models — producing the candidates this page scores
• Buffered Leave-One-Out Cross-Validation — the fold design for clustered samples

← Back to Variogram Diagnostics & Validation