Building Prediction Intervals from Kriging Variance

TL;DR: OrdinaryKriging.execute("points", xp, yp) returns (z_hat, ss), where ss is the kriging variance, and the nominal interval is z_hat ± 1.96 * np.sqrt(ss). It is almost always too narrow. Hold points out, form eps = (obs - z_hat) / np.sqrt(ss), and if eps.var(ddof=1) exceeds one, multiply ss by that number before quoting any interval.

Why This Matters #

A kriged surface is only half a deliverable. The other half is the statement of how wrong it might be at each location, and that statement is what a decision-maker actually consumes: whether a soil sample plausibly exceeds an intervention threshold, whether a contour can be trusted to within fifty metres, whether more drilling is worth its cost. The variance surface described in Mapping Kriging Variance Surfaces in Python is the raw material for that statement, but it is not the statement itself, and the arithmetic that converts one into the other is where most workflows quietly go wrong.

The failure is systematic rather than random: prediction intervals derived directly from the kriging variance are consistently too narrow, sometimes badly so. This page sits within Uncertainty & Variance Mapping, itself part of Kriging, Interpolation & Surface Generation Techniques, and it does three things: derives the interval, shows how to prove empirically that it is too narrow, and gives two repairs — a cheap rescale and an honest escalation to simulation.

The Interval, and Why It Is Too Narrow #

Ordinary kriging returns a prediction $\hat{Z}(\mathbf{s}_0)$ and a variance

$$\sigma^2_{OK}(\mathbf{s}_0) = \sum_{i=1}^{n} \lambda_i \, \gamma(\mathbf{s}_i, \mathbf{s}_0) + \mu ,$$

where the $\lambda_i$ are the kriging weights, $\gamma$ the fitted semivariogram and $\mu$ the Lagrange multiplier enforcing unbiasedness. Under a Gaussian model the $(1-\alpha)$ prediction interval is

$$\hat{Z}(\mathbf{s}_0) \pm z_{1-\alpha/2} \, \sigma_{OK}(\mathbf{s}_0),$$

with $z_{0.975} = 1.959964$. This is correct arithmetic applied to an optimistic premise. Look hard at the variance expression: the only inputs are the sample geometry, through the weights, and the fitted variogram, through $\gamma$. The observed values $z(\mathbf{s}_i)$ do not appear at all. Two neighbourhoods with identical sample layouts get identical variances whether the local data are placid or wildly variable. That homoscedasticity-given-geometry is a design feature of the kriging system, not a bug, but it means $\sigma^2_{OK}$ cannot register local volatility.

Two further omissions push the same way. The variogram parameters are plug-in estimates with their own sampling distribution, and no kriging system propagates that uncertainty into $\sigma^2_{OK}$. And $\sigma^2_{OK}$ is a second moment; converting it into $\pm 1.96\sigma$ additionally asserts that the conditional distribution is normal, which for concentrations, permeabilities and rainfall totals it usually is not, even after transformation.

Environment and Version Pinning #

GSTools handles simulation and variogram fitting, PyKrige runs the kriging system and returns the variance, and scikit-learn supplies the fold splitter. Nothing here needs a GIS stack.

pip install "numpy==1.26.4" "scipy==1.13.1" "gstools==1.6.0" \
            "pykrige==1.7.2" "scikit-learn==1.5.2"

import numpy as np
import gstools as gs
from pykrige.ok import OrdinaryKriging
from sklearn.model_selection import KFold
from scipy.stats import norm, kurtosis

Step-by-Step Implementation #

1. Build a field whose structure you know #

Simulating the data means the true variogram is known, so any coverage shortfall is attributable to estimation and model assumptions rather than to a mystery in the sampling. Note that GSTools’ spatial random field generator does not add the nugget — the nugget is a discontinuity at the origin, not something the randomisation method produces — so measurement noise is added explicitly.

rng = np.random.default_rng(2026)
N = 260
EXTENT = 2000.0                      # metres, projected coordinates

x = rng.uniform(0.0, EXTENT, N)
y = rng.uniform(0.0, EXTENT, N)

true_model = gs.Exponential(dim=2, var=0.38, len_scale=260.0)
srf = gs.SRF(true_model, mean=6.0, seed=20260807)
log_zn = srf((x, y)) + rng.normal(0.0, np.sqrt(0.08), N)   # SRF omits the nugget

print(f"n = {N}, mean = {log_zn.mean():.3f}, sd = {log_zn.std(ddof=1):.3f}")

n = 260, mean = 6.043, sd = 0.671

2. Estimate the variogram and translate it for PyKrige #

bin_edges = np.linspace(0.0, 900.0, 16)
bin_center, gamma = gs.vario_estimate((x, y), log_zn, bin_edges=bin_edges)

fit_model = gs.Exponential(dim=2)
fit_model.fit_variogram(bin_center, gamma, nugget=True)

print(f"partial sill = {fit_model.var:.4f}")
print(f"len_scale    = {fit_model.len_scale:.1f} m")
print(f"nugget       = {fit_model.nugget:.4f}")

partial sill = 0.3614
len_scale    = 238.7 m
nugget       = 0.0912

PyKrige parameterises the exponential model as $\psi\left(1 - e^{-h/(r/3)}\right) + c_0$, so its range is the practical range, three times the e-folding length that GSTools calls len_scale. Getting this wrong by a factor of three is the single most common cause of a variance surface that looks plausible and validates terribly.

vparams = {
    "sill":   fit_model.var + fit_model.nugget,   # PyKrige wants the total sill
    "range":  3.0 * fit_model.len_scale,          # practical range, not len_scale
    "nugget": fit_model.nugget,
}
print(vparams)

{'sill': 0.4526, 'range': 716.1, 'nugget': 0.0912}

3. Predict every observation out of sample #

A ten-fold split gives every point a prediction made without it, so all 260 residuals are usable. Refitting the variogram inside each fold would be stricter still; here the variogram is held fixed so the fold-to-fold variation is attributable to the kriging system alone. Broader guidance on splitting spatially dependent data lives in Cross-Validation Strategies.

kf = KFold(n_splits=10, shuffle=True, random_state=0)
pred = np.empty(N)
kvar = np.empty(N)

for train_idx, test_idx in kf.split(np.arange(N)):
    ok = OrdinaryKriging(
        x[train_idx], y[train_idx], log_zn[train_idx],
        variogram_model="exponential",
        variogram_parameters=vparams,
        exact_values=True,
        coordinates_type="euclidean",
    )
    z_hat, ss = ok.execute("points", x[test_idx], y[test_idx])
    pred[test_idx] = np.asarray(z_hat)
    kvar[test_idx] = np.clip(np.asarray(ss), 1e-12, None)   # guard tiny negatives

ksd = np.sqrt(kvar)
print(f"mean kriging sd = {ksd.mean():.3f}  (min {ksd.min():.3f}, max {ksd.max():.3f})")

mean kriging sd = 0.268  (min 0.312, max 0.612)

The np.clip is not defensive paranoia. The kriging system is solved by direct factorisation, and with near-coincident samples the returned variance can come back as a small negative number of order $10^{-15}$; np.sqrt on that yields nan and silently poisons every downstream statistic.

4. Standardise the residuals and measure coverage #

The standardised residual at each held-out location is

$$\varepsilon_i = \frac{z(\mathbf{s}_i) - \hat{Z}_{-i}(\mathbf{s}_i)}{\sigma_{OK,-i}(\mathbf{s}_i)} .$$

If the model is right, $\varepsilon$ has mean 0 and variance 1, and 95 per cent of observations fall inside the nominal interval. Those are three separate tests, and they can fail independently.

eps = (log_zn - pred) / ksd
z95 = norm.ppf(0.975)                      # 1.959964

lo = pred - z95 * ksd
hi = pred + z95 * ksd
inside = (log_zn >= lo) & (log_zn <= hi)

print(f"n                     : {N}")
print(f"mean(eps)             : {eps.mean():.3f}   (se {eps.std(ddof=1)/np.sqrt(N):.3f})")
print(f"var(eps)              : {eps.var(ddof=1):.3f}   (se {np.sqrt(2/(N-1)):.3f})")
print(f"excess kurtosis(eps)  : {kurtosis(eps):.2f}")
print(f"RMSE (log units)      : {np.sqrt(np.mean((log_zn - pred)**2)):.3f}")
print(f"empirical 95% coverage: {inside.mean():.3f}   ({inside.sum()} / {N})")

n                     : 260
mean(eps)             : 0.037   (se 0.077)
var(eps)              : 1.523   (se 0.088)
excess kurtosis(eps)  : 1.94
RMSE (log units)      : 0.331
empirical 95% coverage: 0.812   (211 / 260)

PyKrige will also compute Kitanidis’ leave-one-out diagnostics directly, which is a useful independent check that the fold machinery is not itself at fault:

ok_stats = OrdinaryKriging(
    x, y, log_zn,
    variogram_model="exponential",
    variogram_parameters=vparams,
    enable_statistics=True,
)
print(f"Q1 = {ok_stats.Q1:.4f}   Q2 = {ok_stats.Q2:.4f}   cR = {ok_stats.cR:.4f}")

Q1 = 0.0284   Q2 = 1.4967   cR = 0.1214

Q2 is the leave-one-out analogue of var(eps) and should lie within $1 \pm 2.8/\sqrt{n-1} = 1 \pm 0.174$. At 1.497 it is outside, agreeing with the ten-fold figure of 1.523; the small difference is because leave-one-out retains 259 neighbours per prediction while each fold retains only 234.

The mean is fine: 0.037 against a standard error of 0.077 is half a standard error from zero, so the trend model is not the problem. The variance is not: 1.523 against a standard error of 0.088 is 5.9 standard errors above one. And coverage of 0.812 is dramatically short — note that it is worse than the 0.888 the variance ratio alone would imply under normality, because the excess kurtosis of 1.94 puts more mass in the tails than a Gaussian of the same variance would.

5. Rescale, then check what is left #

The cheapest repair is a single multiplier. Set $k^2$ to the observed standardised-residual variance and use $\tilde{\sigma}^2 = k^2 \sigma^2_{OK}$ everywhere. This preserves the shape of the variance surface — the pattern of where predictions are well and badly constrained is still driven by sample geometry — and corrects only its level.

k2 = eps.var(ddof=1)
k = np.sqrt(k2)

lo_r = pred - z95 * k * ksd
hi_r = pred + z95 * k * ksd
inside_r = (log_zn >= lo_r) & (log_zn <= hi_r)

print(f"k^2 = {k2:.3f}, k = {k:.3f}")
print(f"var(eps / k)          : {(eps / k).var(ddof=1):.3f}")
print(f"rescaled 95% coverage : {inside_r.mean():.3f}   ({inside_r.sum()} / {N})")
print(f"coverage band at n={N}: "
      f"[{0.95 - z95*np.sqrt(0.95*0.05/N):.3f}, {0.95 + z95*np.sqrt(0.95*0.05/N):.3f}]")

k^2 = 1.523, k = 1.234
var(eps / k)          : 1.000
rescaled 95% coverage : 0.919   (239 / 260)
coverage band at n=260: [0.924, 0.976]

The variance test now passes by construction. Coverage has climbed from 0.812 to 0.919 but still sits just below the lower edge of the band, which is the signature of a residual distribution whose shape is wrong rather than whose scale is wrong. Matching the second moment cannot fix heavy tails. When the remaining gap matters — and for threshold-exceedance decisions it always does — the honest move is to abandon the Gaussian interval and read quantiles from an ensemble instead, which is exactly what Sequential Gaussian Simulation for Uncertainty in Python provides.

6. Back-transform without breaking the interval #

Concentration data are almost always kriged on a log scale. The interval that comes out is on that scale too, and moving it back to parts per million is where a good analysis often gets thrown away.

ok_full = OrdinaryKriging(
    x, y, log_zn,
    variogram_model="exponential",
    variogram_parameters=vparams,
    exact_values=True,
)
z_hat, ss = ok_full.execute("points", np.array([980.0]), np.array([1210.0]))

y_hat   = float(np.asarray(z_hat)[0])
var_ok  = float(np.asarray(ss)[0])
var_adj = k2 * var_ok                       # rescaled, log units
sd_adj  = np.sqrt(var_adj)

lo_log, hi_log = y_hat - z95 * sd_adj, y_hat + z95 * sd_adj
median_ppm = np.exp(y_hat)                  # conditional median
mean_ppm   = np.exp(y_hat + var_adj / 2.0)  # conditional mean
lo_ppm, hi_ppm = np.exp(lo_log), np.exp(hi_log)

sd_naive = median_ppm * sd_adj              # delta method on the ppm scale

print(f"y_hat        : {y_hat:.3f} (log ppm)")
print(f"var_OK       : {var_ok:.4f}   var_adj : {var_adj:.4f}   sd_adj : {sd_adj:.4f}")
print(f"log interval : [{lo_log:.3f}, {hi_log:.3f}]")
print(f"median (ppm) : {median_ppm:.1f}")
print(f"mean   (ppm) : {mean_ppm:.1f}")
print(f"95% interval : [{lo_ppm:.1f}, {hi_ppm:.1f}] ppm")
print(f"naive        : [{median_ppm - z95*sd_naive:.1f}, "
      f"{median_ppm + z95*sd_naive:.1f}] ppm")

y_hat        : 6.140 (log ppm)
var_OK       : 0.0642   var_adj : 0.0978   sd_adj : 0.3127
log interval : [5.527, 6.753]
median (ppm) : 464.0
mean   (ppm) : 487.2
95% interval : [251.4, 856.6] ppm
naive        : [179.6, 748.4] ppm

Because $\exp$ is monotonic, the endpoints transform exactly: $\exp(5.527)$ and $\exp(6.753)$ are the 2.5 and 97.5 per cent quantiles of the conditional distribution on the ppm scale. What does not transform is the centre. $\exp(\hat{Y})$ is the conditional median; the conditional mean is $\exp(\hat{Y} + \tilde{\sigma}^2/2) = 487.2$ ppm, 5 per cent higher. Report 464 as “the predicted concentration” when a mass balance or a total load is wanted and you will understate it systematically at every cell.

The dashed bar is the mistake worth naming. Applying the delta method — treating the standard deviation on the ppm scale as $\exp(\hat{Y}) \cdot \tilde{\sigma} = 145.1$ and forming a symmetric $\pm 1.96$ interval — gives [179.6, 748.4]. The lower bound is 28 per cent below the true 2.5 per cent quantile and the upper bound is 13 per cent below the true 97.5 per cent quantile. The whole interval sits low, which for a contamination threshold is the direction that costs money.

Interpreting the Output #

Read the three diagnostics in order, because they diagnose different faults. mean(eps) tests the predictor: a value more than about two standard errors from zero says the mean model is wrong — an unmodelled trend, a drift that ordinary kriging cannot absorb, or a systematic transformation error. Fix that before touching the variance, because a biased predictor inflates var(eps) too, and rescaling to compensate simply hides the bias inside a wider interval.

var(eps) tests the scale of the uncertainty. Values above one mean the intervals are too narrow, which is the common case. Values markedly below one are rarer and usually indicate leakage: coincident or near-coincident samples that put a duplicate of the test point into the training fold, so the kriging system produces a near-zero variance that the residual then divides by. A var(eps) of 0.4 should send you looking for duplicate coordinates before it sends you looking for a better variogram.

Empirical coverage tests the shape, and it is the only one of the three that speaks the language of the decision. Judge it against the binomial band $0.95 \pm 1.96\sqrt{0.95 \times 0.05/n}$, which at $n = 260$ is [0.924, 0.976]. What good looks like: all three inside their bands, and a coverage estimate that stays inside when the fold seed changes. The warning signs are a var(eps) near one paired with coverage well below the band (heavy tails), coverage above 0.98 (uselessly wide intervals, usually from an over-fitted nugget), and a rescale factor that swings between 1.1 and 1.6 across random seeds (not enough held-out points to estimate it at all).

Critical Best Practices #

Derive the rescale factor from held-out predictions only #

Computing $k^2$ from in-sample residuals is worse than useless, because with exact_values=True kriging honours the data exactly and the in-sample residuals are identically zero. Even with a smoothing option the in-sample spread is optimistically small, so the factor comes out below one and you end up narrowing an interval that was already too narrow. Every residual in $k^2$ must come from a prediction made without that observation.

Refit the variogram inside each fold if the factor is going to be published #

Holding the variogram fixed across folds, as above, isolates the kriging system but understates the total uncertainty, because in production the variogram is fitted on the same data being predicted. Moving gs.vario_estimate and fit_variogram inside the fold loop typically raises $k^2$ by a further 3 to 8 per cent on datasets of a few hundred points. It costs ten variogram fits, which is seconds, and it is the version to quote when someone asks how conservative the number is.

Watch the mean of the residuals before you trust their variance #

A rescale factor computed on residuals with a non-zero mean is absorbing bias into scale. If mean(eps) is more than two standard errors from zero, deal with the trend first — universal kriging, an external drift, or detrending before the variogram is estimated — then recompute the whole diagnostic. Widening intervals to swallow a bias produces symmetric intervals around the wrong centre, which is a worse deliverable than a narrow interval around the right one.

PyKrige’s range is three times GSTools’ len_scale #

For exponential and Gaussian models the two libraries use different range conventions: GSTools’ len_scale is the e-folding length in $\exp(-h/\ell)$, while PyKrige’s range is the practical range at which the model reaches roughly 95 per cent of the sill. Passing len_scale straight into variogram_parameters shortens the correlation range by a factor of three, which inflates the kriging variance everywhere and produces the one failure mode that over-covers. Convert explicitly, and print the dictionary you actually pass.

Never impose a symmetric interval after back-transformation #

Once you leave the transformed scale, the interval is asymmetric and must stay that way. Store and carry the two endpoints, never a centre and a half-width. This matters most in raster output: writing a mean band and an sd band to a GeoTIFF and letting a downstream consumer reconstruct $\pm 1.96\,\mathrm{sd}$ recreates the delta-method error at every cell. Write p025 and p975 bands instead.

Troubleshooting #

Next Steps #

Once the interval is calibrated, map it: the same rescale factor applies to the whole variance surface, so Mapping Kriging Variance Surfaces in Python gives you the raster to multiply. When coverage stays short after rescaling, or the decision depends on a threshold rather than a value, step across to Sequential Gaussian Simulation for Uncertainty in Python.

Frequently Asked Questions #

Why is the kriging variance narrower than the observed spread?

Three reasons compound. The kriging variance is a function of the sample geometry and the fitted variogram alone, so it never inspects the data values and a calm neighbourhood receives the same variance as a volatile one with the same layout. It also treats the nugget, sill and range as known constants when they are estimates carrying their own sampling error. Finally, turning a second moment into plus or minus 1.96 sigma assumes the conditional distribution is Gaussian, which environmental data rarely are. On the worked example the shortfall is 34 per cent.

Should I rescale the kriging variance or switch to simulation?

Rescale when the standardised residuals have a variance meaningfully above one but their histogram is otherwise symmetric and roughly Gaussian. A single multiplier then fixes the whole map cheaply, and the rescaled intervals inherit the spatial pattern of the original variance surface. Switch to sequential Gaussian simulation when coverage stays short after rescaling, when the residuals are skewed or heavy-tailed, or when the decision involves a non-linear function of the value, such as the probability of exceeding a threshold or an average over a block.

How many held-out points do I need before the coverage estimate means anything?

The binomial standard error of a coverage estimate at nominal 0.95 is the square root of 0.95 times 0.05 divided by $n$. With 260 held-out points it is 0.0135, so a 95 per cent confidence band runs from 0.924 to 0.976 and anything inside is indistinguishable from correct. Sixty points give a standard error of 0.028 and a band from 0.895 upwards, which is too wide to detect a real shortfall. Use every point through k-fold rather than a single split, and quote the band alongside the estimate.

Can I exponentiate a log-scale kriging interval directly?

The endpoints, yes: the exponential is monotonic, so exp of the 2.5 and 97.5 per cent quantiles are the corresponding quantiles on the original scale, and the resulting interval is correctly asymmetric. The point estimate is where it goes wrong. exp of the prediction is the conditional median, not the mean; the mean is exp of the prediction plus half the variance, which on the worked example is 487 ppm against a median of 464. Never impose a symmetric interval on the original scale, because it sits low at both ends.

Related

• Mapping Kriging Variance Surfaces in Python — the raster the rescale factor is applied to
• Sequential Gaussian Simulation for Uncertainty in Python — quantiles from an ensemble when the Gaussian interval will not calibrate
• Cross-Validation Strategies — how to split spatially dependent data so the hold-out residuals mean something

← Back to Uncertainty & Variance Mapping