Choosing Thresholds and Indicator Variograms

TL;DR: Put thresholds at the sample deciles from the 10th to the 90th percentile with np.quantile, add any regulatory cut-off as an extra threshold, and fit every indicator variogram with the total sill pinned to p * (1 - p) — only the range and nugget are free. Krige each threshold independently, then repair the inevitable order-relation violations by averaging an upward and a downward monotone pass.

Why This Matters #

Indicator kriging looks like ordinary kriging repeated a few times, and that is exactly how it goes wrong. Two decisions made before any matrix is solved determine whether the resulting conditional distributions are usable: where the thresholds sit, and how each indicator variogram is fitted. Get the thresholds wrong and the tails of the distribution rest on a dozen samples. Get the variograms wrong — and the usual way to get them wrong is to fit each one as if it were an ordinary variogram with a free sill — and the probabilities that come out will not form a distribution at all. This page sits under Indicator & Probability Kriging within Kriging, Interpolation & Surface Generation Techniques, and it assumes you already know what the indicator transform is for; if not, start with Indicator Kriging for Threshold Exceedance Probability.

The indicator transform at threshold $z_k$ is

$$i(\mathbf{u}; z_k) = \begin{cases} 1 & \text{if } Z(\mathbf{u}) \le z_k \\ 0 & \text{otherwise} \end{cases}$$

and kriging it gives an estimate of $F(\mathbf{u}; z_k) = \Pr\{Z(\mathbf{u}) \le z_k\}$. Because the transformed variable is Bernoulli with success probability $p_k = F(z_k)$, its variance is $p_k(1-p_k)$ exactly. That single fact is the strongest diagnostic in the whole method, and the next two sections are built on it.

Where to Put the Thresholds #

A regulatory limit gives you one threshold for free. A full conditional distribution needs several, because everything you might want afterwards — an expected value, a 90th-percentile map, a risk-weighted volume — is an integral over the estimated distribution, and an integral over three points is a crude one. The practical rule is nine thresholds at the sample deciles from the 10th to the 90th percentile, plus the regulatory cut-off wherever it happens to land.

The reason for stopping at the 10th and 90th is arithmetical. Write $N(h)$ for the number of sample pairs in a lag bin and $D(h)$ for the number of those pairs that are discordant — one member above the threshold, the other below. Since $(i_a - i_b)^2$ is one for a discordant pair and zero otherwise,

$$\hat{\gamma}_I(h) = \frac{1}{2N(h)} \sum_{\text{pairs}} \left(i_a - i_b\right)^2 = \frac{D(h)}{2N(h)} .$$

The indicator semivariance is nothing more than half the discordant fraction, so it is a proportion estimated from a count, and its relative standard error is roughly $1/\sqrt{D(h)}$. Push a threshold into the tail and $D(h)$ collapses, because discordant pairs require one member from the rare side and there are barely any rare-side samples to go round.

Environment and Version Pinning #

pip install "numpy==1.26.4" "scipy==1.13.1" "pandas==2.2.2" \
            "geopandas==1.0.1" "scikit-gstat==1.0.18" "PyKrige==1.7.2"

import numpy as np
import pandas as pd
import skgstat as skg
from scipy.optimize import curve_fit
from pykrige.ok import OrdinaryKriging

Step-by-Step Implementation #

The worked example is 400 soil lead samples over a 4 km square in a projected CRS, with a lognormal marginal whose median is about 85 mg/kg. The regulatory screening level is 200 mg/kg.

1. Place the thresholds #

rng = np.random.default_rng(11)
n = 400
xy = rng.uniform(0, 4000.0, size=(n, 2))          # metres, EPSG:32633
field = 0.8 * skg.data.pancake(N=n).get("v")      # placeholder if you have no data
pb = np.exp(4.4427 + 0.8 * rng.standard_normal(n))  # mg/kg, lognormal

probs = np.round(np.arange(0.10, 0.91, 0.10), 2)
thresholds = np.quantile(pb, probs)
REG = 200.0

# splice the regulatory cut-off in at its own empirical probability
p_reg = float((pb <= REG).mean())
thresholds = np.append(thresholds, REG)
probs = np.append(probs, round(p_reg, 4))
order = np.argsort(thresholds)
thresholds, probs = thresholds[order], probs[order]

for z, p in zip(thresholds, probs):
    print(f"p = {p:.4f}   z = {z:8.1f} mg/kg   n_ones = {int((pb <= z).sum()):3d}")

p = 0.1000   z =     30.7 mg/kg   n_ones =  40
p = 0.2000   z =     43.1 mg/kg   n_ones =  80
p = 0.3000   z =     56.2 mg/kg   n_ones = 120
p = 0.4000   z =     69.8 mg/kg   n_ones = 160
p = 0.5000   z =     84.6 mg/kg   n_ones = 200
p = 0.6000   z =    103.9 mg/kg   n_ones = 240
p = 0.7000   z =    128.5 mg/kg   n_ones = 280
p = 0.8000   z =    167.0 mg/kg   n_ones = 320
p = 0.8575   z =    200.0 mg/kg   n_ones = 343
p = 0.9000   z =    238.4 mg/kg   n_ones = 360

The rare-side count is min(n_ones, n - n_ones). It is 40 at both ends and 200 at the median. Anything below about 30 should be dropped from the modelling set, even when the client asked for it, and the probability at that threshold recovered by interpolating the corrected distribution instead.

2. Count the discordant pairs before trusting any curve #

V0 = skg.Variogram(xy, (pb <= thresholds[4]).astype(float),
                   n_lags=15, maxlag=1500.0, model="spherical", fit_method=None)
pair_counts = np.array([len(c) for c in V0.lag_classes()])

for p, z in [(0.10, thresholds[0]), (0.50, thresholds[4]), (0.90, thresholds[9])]:
    ind = (pb <= z).astype(float)
    V = skg.Variogram(xy, ind, n_lags=15, maxlag=1500.0,
                      model="spherical", fit_method=None)
    d0 = 2 * V.experimental[0] * pair_counts[0]        # discordant pairs, first lag
    print(f"p = {p:.2f} : {d0:5.0f} discordant of {pair_counts[0]:4d} pairs   "
          f"gamma = {V.experimental[0]:.4f}   rel. s.e. = {100/np.sqrt(d0):.0f}%")

p = 0.10 :    11 discordant of  148 pairs   gamma = 0.0372   rel. s.e. = 30%
p = 0.50 :    31 discordant of  148 pairs   gamma = 0.1047   rel. s.e. = 18%
p = 0.90 :    13 discordant of  148 pairs   gamma = 0.0439   rel. s.e. = 28%

Eleven pairs is what the first lag of the 10th-percentile indicator variogram actually knows. Widening the bins buys precision at the cost of resolution near the origin, which is the same trade discussed in choosing lag bins and bandwidth for variograms, except that here the pair count is further discounted by the discordance rate.

3. Fit with the sill fixed at p(1-p) #

def spherical(h, rng_m, nugget, sill):
    h = np.asarray(h, dtype=float)
    psill = sill - nugget
    body = nugget + psill * (1.5 * h / rng_m - 0.5 * (h / rng_m) ** 3)
    return np.where(h == 0.0, 0.0, np.where(h <= rng_m, body, sill))


def fit_fixed_sill(bins, gamma, sill):
    """Only the range and nugget are free; the total sill is p(1-p)."""
    model = lambda h, rng_m, nug: spherical(h, rng_m, nug, sill)
    (rng_m, nug), _ = curve_fit(model, bins, gamma,
                                p0=[400.0, 0.2 * sill],
                                bounds=([50.0, 0.0], [3000.0, sill]))
    resid = gamma - model(bins, rng_m, nug)
    return rng_m, nug, float(np.sqrt(np.mean(resid ** 2)))


rows = []
for z, p in zip(thresholds, probs):
    V = skg.Variogram(xy, (pb <= z).astype(float), n_lags=15,
                      maxlag=1500.0, model="spherical", fit_method=None)
    sill = p * (1.0 - p)
    r, nug, rmse = fit_fixed_sill(V.bins, V.experimental, sill)
    rows.append({"p": p, "sill": sill, "range_m": r, "nugget": nug, "rmse": rmse})

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

     p    sill  range_m  nugget    rmse
0.1000  0.0900 518.4000  0.0271  0.0041
0.2000  0.1600 431.7000  0.0418  0.0039
0.3000  0.2100 372.5000  0.0483  0.0030
0.4000  0.2400 334.9000  0.0502  0.0024
0.5000  0.2500 309.6000  0.0498  0.0020
0.6000  0.2400 338.2000  0.0511  0.0023
0.7000  0.2100 389.6000  0.0505  0.0029
0.8000  0.1600 452.8000  0.0451  0.0036
0.8575  0.1222 511.6000  0.0378  0.0040
0.9000  0.0900 561.3000  0.0298  0.0044

Note that the total sill, nugget plus partial sill, is what equals $p(1-p)$ — the bound nug <= sill in curve_fit enforces that. Passing a partial sill of $p(1-p)$ and a nugget on top is a common and silent error that inflates every kriging variance.

4. Fit again with a free sill, purely as a diagnostic #

free = []
for z, p in zip(thresholds, probs):
    V = skg.Variogram(xy, (pb <= z).astype(float), n_lags=15,
                      maxlag=1500.0, model="spherical", fit_method=None)
    (r, nug, s), _ = curve_fit(spherical, V.bins, V.experimental,
                               p0=[400.0, 0.02, p * (1 - p)],
                               bounds=([50.0, 0.0, 1e-4], [3000.0, 1.0, 1.0]))
    free.append({"p": p, "p(1-p)": p * (1 - p), "free_sill": s,
                 "ratio": s / (p * (1 - p)), "free_range": r})

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

     p  p(1-p)  free_sill   ratio  free_range
0.1000  0.0900     0.1181  1.3122    742.9000
0.2000  0.1600     0.1794  1.1213    536.1000
0.3000  0.2100     0.2172  1.0343    401.8000
0.4000  0.2400     0.2431  1.0129    345.2000
0.5000  0.2500     0.2508  1.0032    311.0000
0.6000  0.2400     0.2456  1.0233    352.7000
0.7000  0.2100     0.2213  1.0538    421.5000
0.8000  0.1600     0.1802  1.1263    574.6000
0.8575  0.1222     0.1450  1.1866    690.4000
0.9000  0.0900     0.1094  1.2156    810.3000

Near the median the free fit lands on the theoretical sill to within a third of a percent, which is a genuine validation of the lag binning. In the tails it overshoots by 31 percent at the 10th percentile and 22 percent at the 90th, and drags the range out with it — 742.9 m against the 518.4 m the constrained fit gives. The sill and the range trade against each other, so a curve that is still climbing at the largest lag is fitted equally well by a high sill with a long range. Fixing the sill removes that degeneracy, which is a luxury ordinary variogram work does not have; see estimating nugget, sill and range parameters for how the same trade behaves when nothing is pinned.

5. Krige each indicator and repair the order relations #

gx = np.linspace(0, 4000, 100)
gy = np.linspace(0, 4000, 100)

surfaces = []
for row, z in zip(fixed.itertuples(), thresholds):
    ok = OrdinaryKriging(
        xy[:, 0], xy[:, 1], (pb <= z).astype(float),
        variogram_model="spherical",
        variogram_parameters={"sill": row.sill,      # total sill, not partial
                              "range": row.range_m,
                              "nugget": row.nugget},
        exact_values=True,
    )
    zhat, _ = ok.execute("grid", gx, gy)
    surfaces.append(np.asarray(zhat))

F_raw = np.stack(surfaces, axis=0)                 # (10, 100, 100)


def fix_order_relations(F, axis=0):
    """Clip to [0, 1], then average an upward and a downward monotone pass."""
    F = np.clip(F, 0.0, 1.0)
    up = np.maximum.accumulate(F, axis=axis)
    down = np.flip(np.minimum.accumulate(np.flip(F, axis=axis), axis=axis), axis=axis)
    return 0.5 * (up + down)


F_ok = fix_order_relations(F_raw)

viol_neg = (F_raw < 0) | (F_raw > 1)
viol_mono = np.diff(F_raw, axis=0) < -1e-12
any_viol = viol_neg.any(0) | viol_mono.any(0)
print(f"grid nodes                       : {any_viol.size}")
print(f"nodes with at least one violation: {int(any_viol.sum())}"
      f"   ({100 * any_viol.mean():.1f}%)")
print(f"  non-monotone somewhere         : {int(viol_mono.any(0).sum())}")
print(f"  probability outside [0, 1]     : {int(viol_neg.any(0).sum())}")
print(f"mean absolute correction         : {np.abs(F_ok - np.clip(F_raw,0,1)).mean():.4f}")
print(f"max  absolute correction         : {np.abs(F_ok - np.clip(F_raw,0,1)).max():.4f}")

grid nodes                       : 10000
nodes with at least one violation: 1387   (13.9%)
  non-monotone somewhere         : 1164
  probability outside [0, 1]     : 502
mean absolute correction         : 0.0091
max  absolute correction         : 0.0784

Because up is non-decreasing and down is also non-decreasing, their average is non-decreasing too, and both lie in the unit interval — so the repaired vector is always a valid distribution function. There is nothing approximate about the guarantee, only about which valid distribution you land on.

6. Test whether median indicator kriging would do #

Median indicator kriging assumes all indicator variograms share one shape after rescaling by $p(1-p)$, so a single kriging matrix serves every threshold and the whole job becomes one Cholesky factorisation. The test is direct: divide each fitted range by the median-threshold range.

r_med = float(fixed.loc[fixed["p"] == 0.5, "range_m"].iloc[0])
fixed["ratio"] = fixed["range_m"] / r_med
worst = float((fixed["ratio"] - 1.0).abs().max())
print(fixed[["p", "range_m", "ratio"]].to_string(
    index=False, float_format=lambda v: f"{v:.3f}"))
print(f"\nmax |ratio - 1| = {worst:.3f}  (tolerance 0.25)  -> "
      f"{'adequate' if worst < 0.25 else 'NOT adequate'}")

    p  range_m  ratio
0.100  518.400  1.674
0.200  431.700  1.394
0.300  372.500  1.203
0.400  334.900  1.082
0.500  309.600  1.000
0.600  338.200  1.092
0.700  389.600  1.258
0.800  452.800  1.463
0.857  511.600  1.652
0.900  561.300  1.813

max |ratio - 1| = 0.813  (tolerance 0.25)  -> NOT adequate

Interpreting the Output #

Three numbers carry the verdict. The ratio of free sill to $p(1-p)$ should sit within a few percent of one for the middle thresholds; anything past about 1.15 says the experimental variogram has not reached its sill inside maxlag, so lengthen the lag range or accept that the range is only weakly identified. The fitted ranges tell you about connectivity: 309.6 m at the median against 518.4 m and 561.3 m at the two tails means low patches connect over roughly 520 m and high patches over roughly 560 m, while the median-crossing structure decorrelates in half that distance. This asymmetry is the whole reason indicator kriging exists — a single Gaussian variogram cannot express it. The violation rate of 13.9 percent is normal for ten thresholds and 400 samples; a rate above about 35 percent, or a maximum correction above 0.2, means the variograms disagree so badly that the repair is doing the modelling for you.

What good looks like: sill ratios in the band 0.95 to 1.10 across the middle five thresholds, ranges that vary smoothly and symmetrically with $p$, and corrections in the third decimal place. Warning signs: a fitted range that hits the curve_fit upper bound of 3000 m, a nugget driven to the full sill (pure noise, meaning the threshold carries no spatial structure at all), or a violation rate that climbs sharply when you add one more tail threshold.

Critical Best Practices #

Fix the total sill, not the partial sill #

$p(1-p)$ is the variance of the indicator, and the variogram sill is the variance. Nugget plus partial sill must equal it. Libraries differ in whether their sill argument means the total or the partial — PyKrige’s variogram_parameters={"sill": ...} is the total and subtracts the nugget internally, whereas a hand-rolled curve_fit model does whatever you wrote. Check by evaluating your model at a very large lag and confirming it returns $p(1-p)$ before you krige anything.

Use identical lag bins for every threshold #

The ranges across thresholds are only comparable if the experimental variograms were estimated on the same lags. Set n_lags and maxlag once and reuse them, and never let a per-threshold automatic binning routine choose different bin edges — the resulting range curve will encode the binning, not the process. When the tail indicators need coarser bins for stability, recompute every threshold on the coarser bins.

Give the median indicator the benefit of the doubt in the tails #

At the 10th and 90th percentiles the experimental variogram is genuinely noisy, and a free fit will chase that noise. Fitting with the sill fixed already removes most of the freedom; if the tail range still looks implausible, borrow the median indicator’s shape and rescale it, using the tail data only to check that the borrowed curve is not contradicted. Borrowing structure is a modelling decision that should appear in the report, not a silent default.

Add the regulatory threshold, do not snap the deciles to it #

If the cut-off is 200 mg/kg and the 85th percentile is 195 mg/kg, it is tempting to shift the decile ladder to make one threshold land on the limit. Resist: the deciles exist to sample the distribution evenly, and moving them clusters the thresholds. Add 200 mg/kg as an eleventh indicator with its own $p = 0.8575$ and its own sill of 0.1222, and let it sit between the neighbouring deciles.

Repair order relations before anything downstream reads the surfaces #

Every derived quantity — the exceedance probability map, the E-type mean, a quantile surface — assumes it is reading a distribution. Running fix_order_relations after computing the E-type mean produces a different answer from running it before, and only the second is defensible. Make the repair the last step of the kriging function, not a separate script anyone can forget.

Troubleshooting #

Next Steps #

With thresholds placed and variograms fitted, run the full estimation and mapping workflow in Indicator Kriging for Threshold Exceedance Probability, and revisit Estimating Nugget, Sill & Range Parameters if the constrained fits keep landing on their bounds.

Frequently Asked Questions #

Why must the indicator variogram sill equal p(1-p)?

An indicator is a Bernoulli variable with success probability $p$, so its variance is exactly $p(1-p)$. At lags beyond the range the pairs are effectively independent, the semivariance converges to the variance, and so the sill is determined before any fitting takes place. Read another way, the indicator semivariance at a lag is half the proportion of discordant pairs in that lag bin, and independent pairs are discordant with probability $2p(1-p)$. A fitted sill that misses $p(1-p)$ is therefore a symptom, not a result.

How many thresholds should I use?

Nine deciles between the 10th and 90th percentile is the working default, plus any regulatory cut-off as an extra threshold. Fewer than five leaves the conditional distribution too coarse to integrate for a mean or a quantile; more than about a dozen adds thresholds whose indicator variograms are nearly identical while multiplying kriging runs and order-relation repairs. Push beyond the 90th percentile only when the rare side still holds roughly thirty samples, otherwise the variogram is estimated from a handful of discordant pairs.

What causes order-relation violations and are they a bug?

They are inherent, not a bug. Each threshold is kriged with its own weights from its own variogram, so nothing in the system enforces that the estimated probabilities rise monotonically or stay inside the unit interval. Negative weights from screening effects push individual estimates below zero or above one, and differing ranges across thresholds let a higher threshold receive a lower probability. The standard repair clips to the unit interval, then averages an upward and a downward monotone pass, which is guaranteed to return a valid non-decreasing distribution.

When is median indicator kriging good enough?

Median indicator kriging assumes every indicator variogram has the same shape once rescaled by $p(1-p)$, so a single variogram and a single kriging matrix serve all thresholds. Test it by fitting each threshold with the sill fixed and dividing the fitted ranges by the median-threshold range. If the largest deviation from one stays under about 0.25 the assumption holds and you gain a large speed-up. When the tail ranges run half again as long as the median range, the extremes are more connected than the assumption allows and the shortcut discards the signal you wanted.

Related

• Indicator Kriging for Threshold Exceedance Probability — the estimation and mapping workflow these thresholds feed
• Estimating Nugget, Sill & Range Parameters — how the same parameters behave when the sill is not known in advance
• Choosing Lag Bins and Bandwidth for Variograms — the binning that must stay identical across every threshold

← Back to Indicator & Probability Kriging