Robust Variogram Estimators for Outlier-Heavy Data

TL;DR: Swap the reduction rule, not the data. Compute matheron, cressie and dowd on identical lag bins; if the Matheron plateau is several times the Dowd plateau you have contamination, not a large sill. In SciKit-GStat that is Variogram(coords, vals, estimator='dowd', use_nugget=True). Then find the offending samples in the variogram cloud before you decide anything.

Why This Matters #

The classical Matheron estimator averages squared differences. Squaring is what makes it an unbiased estimator of semivariance under second-order stationarity, and it is also what makes it fragile: a difference twice as large contributes four times as much, and a difference ten times as large contributes a hundred times as much. Worse, a single contaminated sample is not confined to one lag. With $n$ samples it takes part in $n-1$ pairs, and those pairs are spread across every lag class its neighbours happen to fall into. One transcription error does not add a bump to the variogram — it lifts the whole curve, flattens its shape, and leaves you fitting a model to an artefact.

That failure propagates. The variogram is the input to every kriging weight, so an inflated sill inflates the kriging variance everywhere, and a range destroyed by contamination produces a map that is smoother or rougher than the process warrants. This page sits inside Empirical Variogram Estimation and assumes you already have a working empirical variogram; if the curve looks wrong for reasons other than outliers, the wider Variogram Modeling & Semivariance Analysis section covers the other causes.

The Three Estimators #

All three take the same input — the set $N(h)$ of sample pairs whose separation falls in a lag class — and differ only in how they reduce that set to one number. Write $d_{ij} = z(\mathbf{s}_i) - z(\mathbf{s}_j)$ for the difference of a pair.

The Matheron estimator is the mean of the squares, halved:

$$\hat{\gamma}_{M}(h) = \frac{1}{2|N(h)|}\sum_{(i,j) \in N(h)} d_{ij}^{2}$$

The Cressie–Hawkins estimator averages $|d_{ij}|^{1/2}$ — equivalently the fourth root of the squared difference — raises the average to the fourth power, and divides by a bias correction derived for Gaussian differences:

$$\hat{\gamma}_{CH}(h) = \frac{1}{2}\cdot\frac{\left[\dfrac{1}{|N(h)|}\sum_{(i,j) \in N(h)} |d_{ij}|^{1/2}\right]^{4}}{0.457 + \dfrac{0.494}{|N(h)|} + \dfrac{0.045}{|N(h)|^{2}}}$$

Taking the square root before averaging compresses the scale, so a difference ten times the typical one contributes about three times the typical amount rather than a hundred times. The fourth power restores the units. The correction term removes the bias that this round trip introduces, and it depends on $|N(h)|$ because the bias is worse in sparse bins.

The Dowd estimator abandons the mean altogether and uses the median of the squared differences:

$$\hat{\gamma}_{D}(h) = \tfrac{1}{2} \cdot 2.198 \cdot \operatorname{med}\left\{d_{ij}^{2}\right\} = 1.099 \cdot \operatorname{med}\left\{d_{ij}^{2}\right\}$$

The constant follows directly: if $d_{ij} \sim N(0, 2\gamma)$ then $d_{ij}^{2} = 2\gamma \chi^{2}_{1}$, and the median of a chi-square with one degree of freedom is $0.4549$, so $2\gamma = \operatorname{med}\{d^{2}\}/0.4549 = 2.198 \operatorname{med}\{d^{2}\}$. Because the median ignores the size of the values above it, half the pairs in a bin would have to be contaminated before the estimate moved much. That resistance is bought with efficiency: on clean Gaussian data Dowd is the noisiest of the three, which is why it is a diagnostic first and a production estimator second.

Environment and Version Pinning #

Only NumPy is needed for the estimators themselves. SciKit-GStat provides the same three as named options once you move to production, and SciPy is used by its fitter.

pip install "numpy>=1.26,<3" "scipy>=1.11" "scikit-gstat>=1.0.18" "matplotlib>=3.8"

import numpy as np
import skgstat as skg

Step-by-Step Implementation #

1. Build a contaminated dataset #

Simulate a spherical Gaussian random field with a 300 m range on 250 irregularly placed samples, then corrupt five of them with a fixed +25 mg/kg offset — the kind of error a mis-set calibration or a units mix-up produces.

rng = np.random.default_rng(11)
N, SIDE = 250, 1000.0
xy = rng.uniform(0, SIDE, size=(N, 2))

# Pairwise distances, then a spherical covariance: range 300 m,
# partial sill 4.0, nugget 0.2.
D = np.sqrt(((xy[:, None, :] - xy[None, :, :]) ** 2).sum(-1))
a, psill, nugget = 300.0, 4.0, 0.2
u = np.clip(D / a, 0.0, 1.0)
C = np.where(D <= a, psill * (1 - (1.5 * u - 0.5 * u ** 3)), 0.0)
C[np.diag_indices(N)] = psill + nugget

L = np.linalg.cholesky(C + 1e-8 * np.eye(N))
z_clean = 20.0 + L @ rng.standard_normal(N)

z_dirty = z_clean.copy()
bad = np.sort(rng.choice(N, size=5, replace=False))
z_dirty[bad] = z_clean[bad] + 25.0

print(f"n = {N}   clean: mean {z_clean.mean():.2f}  sd {z_clean.std(ddof=1):.2f}  max {z_clean.max():.2f}")
print(f"contaminated: mean {z_dirty.mean():.2f}  sd {z_dirty.std(ddof=1):.2f}  max {z_dirty.max():.2f}")
print("planted indices:", bad.tolist())

n = 250   clean: mean 19.70  sd 1.70  max 24.88
contaminated: mean 20.20  sd 3.87  max 46.53
planted indices: [70, 117, 184, 188, 215]

Two per cent of the samples have moved the mean by half a unit and more than doubled the standard deviation. That alone is a hint, but a univariate summary cannot tell you whether the spread is contamination or genuine heterogeneity.

2. Write the three estimators #

Each takes a one-dimensional array of differences and returns a semivariance, so they are interchangeable in the binning loop that follows.

def matheron(d):
    """Classical estimator: mean of squared differences, halved."""
    return np.sum(d ** 2) / (2 * d.size)


def cressie(d):
    """Cressie-Hawkins: fourth power of the mean of |d|**0.5, bias-corrected."""
    n = d.size
    num = 0.5 * np.mean(np.abs(d) ** 0.5) ** 4
    return num / (0.457 + 0.494 / n + 0.045 / n ** 2)


def dowd(d):
    """Dowd: 2.198 x median of the squared differences, halved -> 1.099."""
    return 1.099 * np.median(d ** 2)


# Twelve differences: eleven ordinary, one contaminated.
demo = np.array([0.3, -1.9, 0.7, -0.2, 2.6, -1.1, 0.9, -3.1, 1.4, -0.5, 1.8, 25.0])
for name, fn in (("matheron", matheron), ("cressie", cressie), ("dowd", dowd)):
    print(f"{name:9s} first 11: {fn(demo[:-1]):7.3f}   all 12: {fn(demo):8.3f}")

matheron  first 11:   1.276   all 12:   27.211
cressie   first 11:   1.324   all 12:    3.865
dowd      first 11:   1.330   all 12:    1.742

3. Compare the three curves lag by lag #

Hold the bins, the pair set and the maximum lag identical so that any difference between the columns is attributable to the estimator alone. Choosing those bins is a separate decision with its own trade-offs, covered in choosing lag bins and bandwidth for variograms.

iu = np.triu_indices(N, k=1)
dist = D[iu]
MAXLAG, N_LAGS = 500.0, 10
edges = np.linspace(0.0, MAXLAG, N_LAGS + 1)
lags = 0.5 * (edges[:-1] + edges[1:])
masks = [(dist > edges[i]) & (dist <= edges[i + 1]) for i in range(N_LAGS)]
counts = np.array([m.sum() for m in masks])


def curve(z, fn):
    d = z[iu[0]] - z[iu[1]]
    return np.array([fn(d[m]) for m in masks])


mat_clean = curve(z_clean, matheron)
mat_dirty = curve(z_dirty, matheron)
cre_dirty = curve(z_dirty, cressie)
dow_dirty = curve(z_dirty, dowd)

print("  lag_m   n_pairs   matheron_clean   matheron   cressie     dowd")
for i in range(N_LAGS):
    print("  %5.0f   %7d   %14.3f   %8.3f   %7.3f   %6.3f"
          % (lags[i], counts[i], mat_clean[i], mat_dirty[i], cre_dirty[i], dow_dirty[i]))

for name, g in (("matheron_clean", mat_clean), ("matheron", mat_dirty),
                ("cressie", cre_dirty), ("dowd", dow_dirty)):
    print("  %-15s plateau %6.3f   inflation x%.2f"
          % (name, g[-4:].mean(), g[-4:].mean() / mat_clean[-4:].mean()))

  lag_m   n_pairs   matheron_clean   matheron   cressie     dowd
     25       237            0.984     12.865     1.790    1.079
     75       670            1.508     15.780     2.846    1.807
    125       922            2.033     14.330     3.180    2.205
    175      1346            2.443     12.372     3.632    2.653
    225      1646            2.429     13.533     3.895    2.849
    275      1802            2.667     14.030     4.150    2.979
    325      1926            2.371     13.494     3.853    2.739
    375      2135            2.669     14.798     4.279    3.029
    425      2163            2.632     16.326     4.327    2.868
    475      2146            2.830     16.438     4.474    2.828
  matheron_clean  plateau  2.626   inflation x1.00
  matheron        plateau 15.264   inflation x5.81
  cressie         plateau  4.233   inflation x1.61
  dowd            plateau  2.866   inflation x1.09

The Matheron column has lost its shape entirely: instead of rising from 0.98 to a plateau near 2.6 it oscillates between 12 and 16 from the very first bin. Dowd tracks the clean curve to within nine per cent at the plateau and reproduces the rise. Cressie–Hawkins sits between them, sixty per cent high — enough to matter, and a good illustration that a compressed average is not the same thing as an order statistic.

4. Find the offending samples in the variogram cloud #

The variogram cloud is every pair plotted individually: separation distance against $\tfrac{1}{2}d_{ij}^{2}$, the quantity that gets averaged. A contaminated sample appears as a horizontal streak, because its difference against every other sample is roughly the size of the contamination regardless of how far apart they are. Ranking samples by the median contribution of the pairs they belong to turns that visual signature into a list.

diff = z_dirty[iu[0]] - z_dirty[iu[1]]
contrib = 0.5 * diff ** 2          # one point of the variogram cloud per pair
inlag = dist <= MAXLAG
global_med = np.median(contrib[inlag])

ratio = np.array([
    np.median(contrib[((iu[0] == i) | (iu[1] == i)) & inlag]) / global_med
    for i in range(N)
])
order = np.argsort(ratio)[::-1]

print(f"global median cloud contribution = {global_med:.3f}")
print("  idx      z    ratio")
for i in order[:8]:
    flag = "   <-- planted" if i in bad else ""
    print("  %3d  %6.2f  %7.1f%s" % (i, z_dirty[i], ratio[i], flag))
print("median ratio %.2f   95th percentile %.2f"
      % (np.median(ratio), np.percentile(ratio, 95)))

global median cloud contribution = 1.229
  idx      z    ratio
  215   46.53    300.8   <-- planted
   70   44.04    243.8   <-- planted
  184   44.94    242.7   <-- planted
  117   44.63    227.6   <-- planted
  188   42.31    207.0   <-- planted
  107   24.88      8.8
  212   22.99      6.7
  102   15.03      6.5

median ratio 0.82   95th percentile 5.20

The five planted samples occupy the top five places, separated from the sixth by a factor of twenty-three. There is no threshold to tune: the gap identifies itself. Sample 107 at ratio 8.8 is the largest legitimate value in the field, and it is exactly the kind of observation you must not discard reflexively.

5. Confirm the diagnosis with an h-scatterplot #

The h-scatterplot draws $z(\mathbf{s}_i)$ against $z(\mathbf{s}_j)$ for the pairs in a single lag class. Well-behaved short-lag data cling to the diagonal; contaminated samples throw off horizontal and vertical arms, and they collapse the correlation. Symmetrising the pairs keeps the plot and the correlation independent of which member of a pair you call $i$.

k = 0                                   # first lag class, 0-50 m
m = masks[k]


def lag_corr(z):
    a_, b_ = z[iu[0]][m], z[iu[1]][m]
    return np.corrcoef(np.r_[a_, b_], np.r_[b_, a_])[0, 1]


involved = (iu[0][m] == bad[:, None]).any(0) | (iu[1][m] == bad[:, None]).any(0)
print(f"lag {lags[k]:.0f} m, {m.sum()} pairs")
print(f"  r on clean values        : {lag_corr(z_clean):.3f}")
print(f"  r on contaminated values : {lag_corr(z_dirty):.3f}")
print(f"  pairs touching a planted sample: {involved.sum()} "
      f"({100 * involved.sum() / m.sum():.1f}%)")

lag 25 m, 237 pairs
  r on clean values        : 0.623
  r on contaminated values : 0.062
  pairs touching a planted sample: 9 (3.8%)

Nine pairs out of 237 take the short-lag correlation from 0.62 to 0.06. That single number is the clearest statement of the problem on this page: the evidence for short-range continuity — the thing a variogram exists to measure — was destroyed by under four per cent of the pairs.

6. Fit a model to each curve #

Fit the same spherical model to each empirical curve with pair-count weights, so the comparison is between estimators rather than between fitters. A coarse grid search keeps this block dependency-free; it takes under a minute.

def spherical(h, nug, psill_, rng_m):
    u_ = np.clip(h / rng_m, 0.0, 1.0)
    return nug + psill_ * (1.5 * u_ - 0.5 * u_ ** 3)


def fit_spherical(g):
    best = None
    for n0 in np.linspace(0.0, 16.0, 161):
        for ps in np.linspace(0.05, 26.0, 260):
            for rm in np.linspace(20.0, 600.0, 59):
                sse = float(np.sum(counts * (g - spherical(lags, n0, ps, rm)) ** 2))
                if best is None or sse < best[0]:
                    best = (sse, n0, ps, rm)
    _, n0, ps, rm = best
    return n0, n0 + ps, rm


print("  estimator        nugget    sill   range_m")
for name, g in (("matheron_clean", mat_clean), ("matheron", mat_dirty),
                ("cressie", cre_dirty), ("dowd", dow_dirty)):
    n0, sill_, rm = fit_spherical(g)
    print("  %-15s  %6.2f  %6.2f  %8.0f" % (name, n0, sill_, rm))

  estimator        nugget    sill   range_m
  matheron_clean     0.80    2.65       270
  matheron          12.30   15.76       600
  cressie            2.30    4.35       430
  dowd               0.80    2.85       240

7. Use the pairwise-relative estimator when the spread tracks the mean #

A different pathology needs a different tool. When the local variance grows with the local mean — the proportional effect, endemic in assay, rainfall and pollutant data — high-grade neighbourhoods dominate the variogram simply because they are high grade. The pairwise-relative estimator normalises each squared difference by the squared mean of its own pair:

$$\hat{\gamma}_{PR}(h) = \frac{1}{2|N(h)|}\sum_{(i,j) \in N(h)} \frac{d_{ij}^{2}}{\left[\tfrac{1}{2}\left(z(\mathbf{s}_i) + z(\mathbf{s}_j)\right)\right]^{2}}$$

def pairwise_relative(z):
    d = z[iu[0]] - z[iu[1]]
    mid = 0.5 * (z[iu[0]] + z[iu[1]])
    return np.array([np.sum((d[m] ** 2) / (mid[m] ** 2)) / (2 * m.sum())
                     for m in masks])


prv_clean = pairwise_relative(z_clean)
prv_dirty = pairwise_relative(z_dirty)
print("  lag_m   prv_clean   prv_contaminated")
for i in range(N_LAGS):
    print("  %5.0f     %7.5f            %7.5f" % (lags[i], prv_clean[i], prv_dirty[i]))
print("plateau %.5f x mean^2 (%.1f) = %.3f  vs Matheron on clean data %.3f"
      % (prv_clean[-4:].mean(), z_clean.mean() ** 2,
         prv_clean[-4:].mean() * z_clean.mean() ** 2, mat_clean[-4:].mean()))

  lag_m   prv_clean   prv_contaminated
     25     0.00258            0.01447
     75     0.00397            0.01780
    125     0.00532            0.01710
    175     0.00628            0.01579
    225     0.00623            0.01683
    275     0.00677            0.01761
    325     0.00605            0.01657
    375     0.00684            0.01822
    425     0.00671            0.01973
    475     0.00720            0.02009
plateau 0.00670 x mean^2 (388.1) = 2.600  vs Matheron on clean data 2.626

On this field, whose mean is constant, the pairwise-relative curve is just the Matheron curve divided by the squared mean — 0.00670 times 388.1 gives 2.600 against 2.626, as it should. That identity is the check: if dividing the pairwise-relative sill by the squared mean does not reproduce the raw sill, you have a genuine proportional effect. Note also that the estimator is not robust; contamination still inflates it two and a half fold, because a value of 46 raises the numerator far faster than the denominator.

Interpreting the Output #

The number to read first is the ratio of the Matheron plateau to the Dowd plateau. Here it is $15.264 / 2.866 = 5.3$. On clean, roughly Gaussian data that ratio sits close to one, typically inside 0.9 to 1.15; anything above about 1.5 says a small number of extreme differences are carrying the classical estimate. The ratio is a diagnostic, not a correction — it tells you to go looking, not what to conclude.

The fitted parameters make the damage concrete. The clean-data fit gives nugget 0.80, sill 2.65 and range 270 m. Dowd applied to the contaminated data returns 0.80, 2.85 and 240 m: the same model, recovered from data with five gross errors still in it. Cressie–Hawkins returns 2.30, 4.35 and 430 m — a plausible-looking model that is wrong in every parameter, and wrong in the most dangerous way, because nothing about it advertises a problem. Matheron returns nugget 12.30, sill 15.76 and a range that walks to the 600 m ceiling of the search: the curve has so little shape left that the fitter cannot locate a range at all.

What good looks like: three curves within roughly ten per cent of each other across all lags, all rising from a small nugget to a stable plateau. Warning signs, in order of severity: a Matheron curve that is high at the very first bin and stays flat, a Matheron-to-Dowd ratio above 1.5, a robust curve that rises while the classical one does not, and a fitted range that lands on a search boundary. Any of these means the diagnosis belongs in the cloud, not in the fitting routine. If the shape is wrong for other reasons — anisotropy, a trend, too few pairs per bin — diagnosing a bad variogram fit covers the alternatives.

Critical Best Practices #

Compute all three; report the one you can defend #

Running the three estimators over the same bins costs one extra pass over the pair array and turns a single opaque number into a comparison. Make it a standing part of the empirical variogram step, not something you reach for once a map already looks wrong. On clean Gaussian data Matheron is the most efficient of the three, so the comparison also protects you from reaching for a robust estimator when you do not need one.

A robust estimator hides the problem; it does not fix it #

This is the discipline point and it is worth stating flatly. Dowd gave a variogram indistinguishable from the clean-data fit while five samples reading 42 to 47 mg/kg were still sitting in the array. Those values will go into the kriging system, appear in cross-validation residuals and be drawn on the map. The robust curve buys you a trustworthy model while you investigate; it is not the end of the investigation. Name the flagged samples, check them against the field sheets or the laboratory records, classify each as an error or a real extreme, and write the decision down alongside the model parameters.

Never delete a high value because it is high #

A 300-fold cloud ratio is evidence of a data problem. A ratio of 8.8, as sample 107 shows here, is what the tail of a skewed distribution looks like — and in ore-grade, contaminant and rainfall work the tail is usually the part that matters. Deleting real extremes biases the sill downward and produces confidently wrong kriging variances over exactly the locations you care about. The rule is: robust estimation is reversible, deletion is not.

Set use_nugget=True when you move to SciKit-GStat #

The estimators here exist as named options in SciKit-GStat, which is where this belongs in production. The trap is that Variogram fits without a nugget unless you ask for one, so a curve with a real nugget will be fitted with a distorted range and sill.

V = skg.Variogram(xy, z_dirty, estimator="dowd", model="spherical",
                  bin_func="even", n_lags=10, maxlag=500.0, use_nugget=True)
print(V.parameters)          # [range, sill, nugget] for the fitted model

The full workflow, including bin functions and fit methods, is in fitting empirical variograms with SciKit-GStat.

Watch the pair count when using Dowd #

The median needs pairs. A bin holding thirty pairs gives a Dowd estimate with visibly more scatter than a Matheron estimate on the same bin, and below about ten pairs the median is nearly meaningless. Since robust estimation and thin bins pull in opposite directions, prefer wider bins when you are working robustly — check counts before you trust the first lag, where pair numbers are always smallest.

Troubleshooting #

Next Steps #

Once the estimator is chosen and the flagged samples resolved, take the empirical curve through to a fitted model with fitting empirical variograms with SciKit-GStat, and if the fit still refuses to settle, work through diagnosing a bad variogram fit.

Frequently Asked Questions #

Should I always use a robust estimator?

No. On clean, roughly Gaussian data the Matheron estimator is the most efficient of the three and the robust alternatives buy nothing. Compute all three routinely as a diagnostic instead: when they agree, report the classical curve, and when they diverge you have learned something about the data that a single curve would have hidden. The robust estimator earns its place only once you have confirmed the divergence is caused by a small number of extreme differences rather than by genuine short-range structure.

What is the difference between the Cressie–Hawkins and Dowd estimators?

Cressie–Hawkins averages the square root of the absolute difference and raises the mean to the fourth power, with a bias correction that assumes Gaussian differences. It is still an average, so a large enough outlier still moves it. Dowd takes the median of the squared differences and rescales by 2.198, which is a true order statistic with a breakdown point near fifty per cent of the pairs. Dowd is the more resistant of the two and the more variable when the data are clean.

Does a robust estimator fix contaminated data?

It does not. A robust estimator changes how the differences are summarised, so the contaminated values still sit in the dataset and will be carried into kriging, cross-validation and any map you produce. Its role is to give you a variogram you can trust while the data problem is being investigated. Use the divergence between estimators as a trigger to find the offending samples, then repair, remove or retain them explicitly and record which choice you made.

When should I use the pairwise-relative estimator instead?

Use it when the local spread of the variable grows with the local mean, the proportional effect common in assay, rainfall and pollutant data. The pairwise-relative estimator divides each squared difference by the square of the mean of the pair, producing a dimensionless variogram whose sill is roughly the sill of the raw variogram divided by the squared mean. It is not a substitute for a robust estimator, because a single very large value inflates numerator and denominator at different rates.

Related

• Empirical Variogram Estimation — the binning and pair-counting these estimators sit on top of
• Choosing Lag Bins and Bandwidth for Variograms — why a robust estimator needs wider bins than a classical one
• Fitting Empirical Variograms with SciKit-GStat — taking the chosen curve through to a fitted model

← Back to Empirical Variogram Estimation