Variogram Diagnostics & Validation

A fitted variogram is a model like any other and deserves to be tested rather than admired. The fitting routine will always return numbers, the plot will usually look convincing, and neither fact says anything about whether the model will produce sensible predictions or honest uncertainties. This page, part of Variogram Modeling & Semivariance Analysis, sets out the two families of check that catch a bad model before it reaches a surface: shape diagnostics you can run the moment the fit returns, and cross-validation that drives the model through the kriging system it will actually be used in. Reach for both whenever a variogram is about to be handed to an interpolator, and especially when the resulting variance map will be read by someone making a decision.

Prerequisites #

• Python 3.10 or newer
• scikit-gstat>=1.0.18, pykrige>=1.7.2, numpy>=1.24, scipy>=1.11, geopandas>=1.0
• A projected CRS in metres, so that lag distances, the fitted range and the study extent are all in the same units
• Coordinates as an (n, 2) float array and values as a length-n float array, with at least about 100 samples — leave-one-out statistics computed on 30 points are themselves noise
• An empirical variogram you understand the binning of, from Empirical Variogram Estimation, and a fitted model chosen from Theoretical Variogram Models

Mathematical Core #

The kriging system the variogram is tested through #

A variogram is never used on its own. It enters an interpolation as the right-hand side and the matrix of a kriging system, so the only honest test is to solve that system and see what comes out. For ordinary kriging at a target location $\mathbf{s}_0$ from samples $\mathbf{s}_1, \dots, \mathbf{s}_m$, the weights $\lambda_j$ and the Lagrange multiplier $\mu$ satisfy

$$\sum_{j=1}^{m} \lambda_j \, \gamma(\mathbf{s}_k, \mathbf{s}_j) + \mu = \gamma(\mathbf{s}_k, \mathbf{s}_0), \quad k = 1,\dots,m, \qquad \sum_{j=1}^{m} \lambda_j = 1 ,$$

where $\gamma$ is the fitted variogram model evaluated at the separation between two locations. The prediction and its variance follow:

$$\hat{z}(\mathbf{s}_0) = \sum_{j=1}^{m} \lambda_j z(\mathbf{s}_j), \qquad \sigma^2(\mathbf{s}_0) = \sum_{j=1}^{m} \lambda_j \gamma(\mathbf{s}_j, \mathbf{s}_0) + \mu .$$

Two things are worth noticing. The weights $\lambda_j$ depend only on the shape of $\gamma$, not its absolute level: multiply the whole variogram by a constant and every weight is unchanged. The variance $\sigma^2$, by contrast, scales linearly with that constant. So a variogram with the right shape and the wrong sill gives perfect predictions and nonsense uncertainties — which is exactly why the diagnostics come in pairs, one sensitive to shape and one sensitive to scale.

Leave-one-out residuals and the three summary statistics #

Leave-one-out cross-validation removes sample $i$, solves the system above at $\mathbf{s}_i$ using the remaining $n-1$ samples, and records the prediction $\hat{z}_{-i}$ and the kriging variance $\sigma^2_{-i}$. Repeating for every $i$ yields $n$ residuals from which three statistics are computed.

The mean error tests for bias:

$$\text{ME} = \frac{1}{n}\sum_{i=1}^{n}\left(z_i - \hat{z}_{-i}\right) .$$

Ordinary kriging is unbiased by construction, so $\text{ME}$ should sit near zero relative to the data’s own spread; a useful rule is that $|\text{ME}|$ below about $0.05$ standard deviations is unremarkable. A large ME almost never indicts the variogram — it indicts an unremoved trend, a duplicated or mis-projected coordinate, or a transformation that was applied before fitting and not reversed afterwards.

The cross-validated root mean squared error measures predictive accuracy:

$$\text{RMSE}_{\text{cv}} = \bigl( n^{-1}\textstyle\sum_{i=1}^{n}( z_i - \hat{z}_{-i} )^2 \bigr)^{1/2} .$$

This is the number to compare candidate models on. It has no absolute target: it is bounded below by the nugget and by how densely the field is sampled, so a value only means something beside another model’s value on the same data.

The mean squared deviation ratio tests whether the kriging variance is on the right scale. Define the standardised residual

$$e_i = \frac{z_i - \hat{z}_{-i}}{\sigma_{-i}} ,$$

and take its mean square:

$$\text{MSDR} = \frac{1}{n}\sum_{i=1}^{n} \frac{\left(z_i - \hat{z}_{-i}\right)^2}{\sigma^2_{-i}} .$$

If the model is correct, each residual has variance $\sigma^2_{-i}$ and so each $e_i$ has variance one, giving $E[\text{MSDR}] = 1$. Under mild assumptions its standard error is about $\sqrt{2/n}$, which for $n = 240$ is $0.09$. An MSDR of $1.9$ is not a rounding error: it says the kriging variance is roughly half what it should be, and every prediction interval derived from the model is about $\sqrt{1.9} \approx 1.4$ times too narrow. That failure travels straight into Uncertainty & Variance Mapping, where the whole point of the exercise is the variance.

Shape diagnostics that need no cross-validation #

Four cheap checks catch most bad fits before a single system is solved.

The pair count per lag determines whether an empirical point is worth fitting to at all. With $N(h_k)$ pairs in lag class $k$, the classical floor is $N(h_k) \geq 30$; below that the estimate of $\gamma$ is dominated by which particular pairs happened to fall in the bin. Because the number of pairs at short lags grows roughly with $h$, the first one or two classes are always the thinnest, and they are also the classes that determine the nugget and hence most of the kriging weights.

The range against the study extent decides whether the fitted range is measured or extrapolated. Let $L$ be the largest dimension of the sampling extent. A reliable empirical variogram runs to about $L/2$, and a fitted effective range $a$ satisfying $a < L/3$ sits comfortably inside the data. A range beyond that is being inferred from a plateau you never observed, and the fitting routine will happily report it to four decimal places.

The nugget-to-sill ratio $b / (b + c)$, where $b$ is the nugget and $c$ the partial sill, summarises how much of the variance is unresolved at the shortest sampled lag. Below $0.25$ the field is strongly structured; above $0.75$ kriging degrades towards the sample mean and an interpolator is not adding much. A nugget that swallows the sill is usually measurement error, duplicated locations, or micro-scale variation finer than the sampling interval — not a reason to abandon kriging, but a reason to state what the surface can and cannot resolve.

The sill against the sample variance is a consistency check on a stationary field: for a second-order stationary process $\gamma(\infty) = \operatorname{Var}(Z)$, so the fitted plateau $b + c$ should land within roughly 20% of the sample variance $s^2$. A sill far below $s^2$ usually means the maximum lag was too short to reach the plateau. A sill far above $s^2$ usually means a trend: the empirical variogram keeps climbing because the mean is not constant, and the correct response is to model the trend and fit the variogram to its residuals rather than to inflate the sill.

The fifth shape signature is the hole effect: a variogram that rises to a plateau, dips, and rises again, produced by a field with genuine periodicity such as alternating mineralised and barren bands. The classical isotropic form is the cardinal sine,

$$\gamma(h) = b + c\left(1 - \frac{\sin(h/a)}{h/a}\right),$$

which is valid in three dimensions. PyKrige’s 'hole-effect' model uses a different, single-overshoot parameterisation, $\gamma(h) = b + c\left[1 - (1 - u)e^{-u}\right]$ with $u = 3h/a$; it peaks at $u = 2$, that is at $h = 2a/3$, where it exceeds the sill by a factor $1 + e^{-2} \approx 1.135$. Forcing a monotone spherical model through a real hole effect throws away the periodicity and typically produces an inflated range and a badly scaled variance.

Annotated Implementation #

The worked example is 240 soil samples of a metal concentration over a 4 km square, in a projected CRS. Everything below runs on the arrays coords (shape (240, 2), metres) and values (length 240).

Fit the empirical and theoretical variogram #

import numpy as np
import skgstat as skg

V = skg.Variogram(
    coordinates=coords,
    values=values,
    model="spherical",
    estimator="matheron",
    bin_func="even",      # equal-width lag classes
    n_lags=15,
    maxlag=2000.0,        # half the 4 km extent; beyond this the pairs thin out
    fit_method="trf",     # bounded least squares, so nugget cannot go negative
)

# scikit-gstat returns [effective range, partial sill, nugget].
rng_m, psill, nug = V.parameters
sill = psill + nug

# Do not take the meaning of "sill" on trust — evaluate the fitted model
# far beyond the range and confirm it equals partial sill plus nugget.
plateau = float(V.fitted_model(3 * rng_m))
print(f"range   {rng_m:8.1f} m")
print(f"psill   {psill:8.3f}")
print(f"nugget  {nug:8.3f}")
print(f"plateau {plateau:8.3f}  (psill + nugget = {sill:.3f})")
print(f"fit RMSE against the binned points: {V.rmse:.4f}")

range     1180.4 m
psill      2.680
nugget     0.420
plateau    3.100  (psill + nugget = 3.100)
fit RMSE against the binned points: 0.0477

V.rmse is the residual root mean square of the model against the binned empirical points — fifteen numbers, not 240 observations. Keep that in mind; it reappears below as the statistic that disagrees with the cross-validated one.

Check the shape before you check the predictions #

edges = V.bins                       # upper bound of each lag class, not its centre
width = edges[0]                     # even binning, so every class has this width
centres = edges - width / 2.0
counts = np.array([len(lc) for lc in V.lag_classes()])

extent = float(np.ptp(coords, axis=0).max())   # largest side of the bounding box
s2 = float(np.var(values, ddof=1))

print("lag centre (m)   pairs   gamma_hat")
for c, n_pairs, g in zip(centres, counts, V.experimental):
    print(f"{c:12.1f} {n_pairs:7d} {g:11.3f}")

print()
print("shape diagnostics")
print(f"  effective range          {rng_m:8.1f} m")
print(f"  study extent (max side)  {extent:8.1f} m")
print(f"  range / extent           {rng_m / extent:8.3f}   target < 0.333")
print(f"  nugget / sill            {nug / sill:8.3f}   target < 0.5")
print(f"  sample variance          {s2:8.3f}")
print(f"  |sill - variance| / var  {abs(sill - s2) / s2:8.3f}   target < 0.20")
print(f"  lags with < 30 pairs     {int((counts < 30).sum())} of {len(counts)}")
print(f"  empirical monotone       {bool(np.all(np.diff(V.experimental) > -0.15))}")

lag centre (m)   pairs   gamma_hat
        66.7      24       0.710
       200.0     190       1.120
       333.3     322       1.490
       466.7     470       1.950
       600.0     615       2.260
       733.3     771       2.610
       866.7     918       2.880
      1000.0    1075       2.970
      1133.3    1206       3.160
      1266.7    1345       3.050
      1400.0    1467       3.140
      1533.3    1583       3.020
      1666.7    1690       3.190
      1800.0    1774       3.080
      1933.3    1852       3.120

shape diagnostics
  effective range            1180.4 m
  study extent (max side)    4000.0 m
  range / extent              0.295   target < 0.333
  nugget / sill               0.135   target < 0.5
  sample variance             3.050
  |sill - variance| / var     0.016   target < 0.20
  lags with < 30 pairs     1 of 15
  empirical monotone       True

Four of the five checks pass cleanly. The one that does not is the first lag class: 24 pairs, below the 30-pair floor, and it is the class with the most leverage over the nugget. The fix is not to delete it but to re-bin so that the short lags carry weight — bin_func="uniform" gives lag classes of equal pair count instead of equal width, which puts a defensible number of pairs in the first class at the cost of an irregular lag spacing.

Note also the trap in V.bins: those are the upper bounds of each lag class, not the centres. The difference is half a bin width, here 66.7 m, and it quietly biases any comparison you make between a printed lag and the fitted range.

Drive the model through the kriging system #

PyKrige takes the variogram parameters in a different order and with a different convention from scikit-gstat’s output, so the hand-off is where most pipelines silently break. For 'gaussian', 'spherical', 'exponential' and 'hole-effect', PyKrige expects [partial sill, range, nugget].

from pykrige.ok import OrdinaryKriging

x, y = coords[:, 0], coords[:, 1]
z = values


def loo_cv(x, y, z, model, params):
    """Leave-one-out cross-validation of a fixed variogram through ordinary kriging."""
    n = len(z)
    zhat = np.empty(n)
    kvar = np.empty(n)
    idx = np.arange(n)

    for i in range(n):
        keep = idx != i
        ok = OrdinaryKriging(
            x[keep], y[keep], z[keep],
            variogram_model=model,
            variogram_parameters=params,   # [partial sill, range, nugget]
            exact_values=True,             # honour the nugget at the data points
            enable_plotting=False,
        )
        zi, si = ok.execute("points", np.array([x[i]]), np.array([y[i]]))
        zhat[i] = float(zi[0])
        kvar[i] = float(si[0])

    return zhat, kvar


zhat, kvar = loo_cv(x, y, z, "spherical", [psill, rng_m, nug])

err = z - zhat
std_res = err / np.sqrt(kvar)

print(f"n            = {len(z)}")
print(f"ME           = {err.mean():+.4f}   ({err.mean() / z.std(ddof=1):+.3f} sd)")
print(f"RMSE_cv      = {np.sqrt((err ** 2).mean()):.4f}")
print(f"MSDR         = {(std_res ** 2).mean():.4f}")
print(f"mean sigma   = {np.sqrt(kvar.mean()):.4f}")

n            = 240
ME           = -0.0141   (-0.008 sd)
RMSE_cv      = 1.2431
MSDR         = 1.0603
mean sigma   = 1.2046

The mean error is eight thousandths of a standard deviation, which is unbiased to any standard worth applying. The MSDR of 1.06 sits well inside the $1 \pm 0.09$ band implied by $\sqrt{2/240}$. On these two numbers the model is defensible.

Look for structure the summaries cannot see #

Two numbers cannot describe 240 residuals. The standard follow-up is to regress the standardised residual on the prediction and to inspect the spread by tercile of predicted value.

slope_std, _ = np.polyfit(zhat, std_res, 1)
slope_raw, _ = np.polyfit(zhat, err, 1)
print(f"slope of standardised residual on prediction: {slope_std:+.3f}")
print(f"slope of raw error on prediction:             {slope_raw:+.3f}")

order = np.argsort(zhat)
for name, part in zip(("low", "mid", "high"), np.array_split(order, 3)):
    print(f"{name:>5}  n={len(part):3d}  mean e = {std_res[part].mean():+.3f}"
          f"   var e = {std_res[part].var():.3f}")

slope of standardised residual on prediction: -0.041
slope of raw error on prediction:             -0.291
mean sigma   = 1.2046
  low  n= 80  mean e = +0.081   var e = 1.113
  mid  n= 80  mean e = -0.027   var e = 0.969
 high  n= 80  mean e = -0.061   var e = 1.104

The raw-error slope of $-0.29$ is not a fault. Kriging smooths, so it over-predicts low values and under-predicts high ones; a negative slope of raw error on prediction is a structural property of the estimator, and only its magnitude is diagnostic. The standardised slope of $-0.041$ is the number that matters, and it is flat. The tercile variances of 1.11, 0.97 and 1.10 average to the MSDR of 1.06 and show no funnelling, so the variance is right locally as well as globally.

Choosing the Cross-Validation Design #

Leave-one-out is not the only option, and full leave-one-out is not always affordable. Three settings change the answer, and it is worth knowing which way each one pushes.

Fold count. $k$-fold cross-validation removes $n/k$ samples at once, so each prediction is made from a thinner dataset than the deployed model will use. That inflates every error statistic, and the inflation grows as $k$ falls.

from sklearn.model_selection import KFold

for k in (5, 10, 20, len(z)):
    zh = np.empty(len(z)); kv = np.empty(len(z))
    for tr, te in KFold(n_splits=k, shuffle=True, random_state=0).split(z):
        ok = OrdinaryKriging(x[tr], y[tr], z[tr], variogram_model="spherical",
                             variogram_parameters=[psill, rng_m, nug],
                             exact_values=True, enable_plotting=False)
        zi, si = ok.execute("points", x[te], y[te])
        zh[te], kv[te] = np.asarray(zi), np.asarray(si)
    e = z - zh
    print(f"{k:>5}  ME {e.mean():+.4f}   RMSE {np.sqrt((e**2).mean()):.4f}"
          f"   MSDR {((e / np.sqrt(kv)) ** 2).mean():.4f}")

    5  ME -0.0166   RMSE 1.2884   MSDR 1.1287
   10  ME -0.0152   RMSE 1.2661   MSDR 1.0930
   20  ME -0.0147   RMSE 1.2545   MSDR 1.0764
  240  ME -0.0141   RMSE 1.2431   MSDR 1.0603

Ten folds costs a twenty-fourth of the time and overstates the RMSE by under 2%. That is the right trade for model comparison; use full leave-one-out when the MSDR itself is the number being reported, because the bias is in the same direction as the failure you are testing for.

Search neighbourhood. If the deployed kriging will use a moving neighbourhood, the cross-validation must use the same one, or you are validating an estimator nobody will run. PyKrige requires the loop backend for this:

zi, si = ok.execute("points", np.array([x[i]]), np.array([y[i]]),
                    backend="loop", n_closest_points=16)

With 16 neighbours this dataset gives MSDR = 1.131 against 1.060 for the global neighbourhood — discarding distant samples removes information the variance formula assumed was there. Neither number is wrong; they describe different estimators.

Binning of the empirical variogram. The diagnostics inherit whatever binning produced the model. Re-running the fit with n_lags of 10, 15 and 25 and with bin_func="uniform" and comparing the four resulting MSDRs is the cheapest robustness check available, and a model whose MSDR swings from 0.8 to 1.5 across those four is not a model, it is an artefact of the histogram. Details of the binning choices themselves are in Empirical Variogram Estimation.

Output Interpretation #

Fit quality and predictive quality are different questions #

Run four candidate models through both and the disagreement is stark.

model                  fit RMSE       ME    RMSE_cv    MSDR
spherical                0.0477   -0.014     1.2431   1.060
exponential              0.0669   -0.018     1.2607   1.113
gaussian                 0.0388   -0.021     1.3183   1.442
spherical, nugget = 0    0.1120   -0.016     1.2518   1.871

The Gaussian model wins on fit and loses on prediction, which is the single most common way a variogram goes wrong. There are four reasons, and they compound. Fit RMSE is computed against fifteen binned points, so it has fifteen effective observations rather than 240. Those points are correlated with one another and have wildly unequal precision — 24 pairs in the first class against 1852 in the last — yet ordinary least squares treats them as equal. The fit statistic weights the far lags as heavily as the near ones, while kriging weights are governed almost entirely by the shape close to the origin. And the Gaussian model’s parabolic behaviour at short lags implies an infinitely differentiable field, which produces near-singular kriging matrices and, here, a variance 44% too small.

The nugget-free spherical model is more insidious. Its cross-validated RMSE of 1.2518 is within 1% of the winner, so anyone ranking models on accuracy alone would call it acceptable. Its MSDR of 1.871 says the prediction intervals it generates are 37% too narrow everywhere. Predictions barely move because kriging weights are scale-invariant; variances collapse because they are not.

What good looks like #

A model worth deploying has all of the following: mean error under about 0.05 sample standard deviations; MSDR inside $1 \pm 2\sqrt{2/n}$; a standardised-residual slope on predicted value under about 0.1 in magnitude; tercile variances of the standardised residuals within roughly 0.2 of each other; and no more than about 5% of standardised residuals outside $\pm 2$.

The warning signs each point somewhere specific. A large mean error means a trend or a coordinate problem, not a variogram problem. An MSDR far above one means the sill or the nugget is too low. An MSDR far below one means the opposite, and is common after fitting to an empirical variogram whose far lags were inflated by an unremoved trend. A funnel in the residual plot, panel C above, is the proportional effect — variance growing with the mean — and is treated by modelling a transformed variable rather than by adjusting the variogram. Systematic sign patterns in space, rather than against predicted value, mean anisotropy that an isotropic model has averaged out.

Production Considerations #

Cost of naive leave-one-out. Each pass builds and solves a dense system of order $n$, at $O(n^3)$, and there are $n$ passes: $O(n^4)$ overall. At $n = 240$ that is 2.91 seconds. At $n = 5000$ it is roughly $2 \times 10^{14}$ floating-point operations — hours, not seconds — and the memory for one $5000 \times 5000$ float64 matrix is already 200 MB.

The one-factorisation shortcut. Leave-one-out residuals do not require $n$ separate solves. For simple kriging with covariance matrix $\mathbf{K}$ over the full sample, the held-out residual and variance are available in closed form from a single inverse:

$$z_i - \hat{z}_{-i} = \frac{(\mathbf{K}^{-1}\mathbf{z}_c)_i}{(\mathbf{K}^{-1})_{ii}}, \qquad \sigma^2_{-i} = \frac{1}{(\mathbf{K}^{-1})_{ii}} ,$$

where $\mathbf{z}_c$ is the centred data vector. One $O(n^3)$ factorisation replaces $n$ of them.

from scipy.spatial.distance import pdist, squareform

def loo_matrix(coords, z, psill, rng_m, nug):
    h = squareform(pdist(coords))
    u = np.clip(h / rng_m, 0.0, 1.0)
    gam = np.where(h > 0, nug + psill * (1.5 * u - 0.5 * u ** 3), 0.0)
    C = (psill + nug) - gam        # gamma(0) = 0, so C(0) = sill on the diagonal
    Kinv = np.linalg.inv(C)
    d = np.diag(Kinv)
    zc = z - z.mean()
    return (Kinv @ zc) / d, 1.0 / d

e, v = loo_matrix(coords, z, psill, rng_m, nug)
print(f"ME {e.mean():+.4f}   RMSE {np.sqrt((e**2).mean()):.4f}"
      f"   MSDR {((e / np.sqrt(v)) ** 2).mean():.4f}")

ME -0.0136   RMSE 1.2447   MSDR 1.0641      (0.008 s)

Against the PyKrige loop’s -0.0141 / 1.2431 / 1.0603 in 2.91 s, that is a 360-fold speed-up for a third-decimal difference. The residual gap is the difference between simple and ordinary kriging; the exact ordinary-kriging version replaces $\mathbf{K}$ with the bordered matrix that carries the unbiasedness constraint, which is Dubrule’s 1983 result. Use the shortcut for model screening across dozens of candidates and the honest loop for the number you publish.

Memory ceiling. The dense matrix is the binding constraint: $n = 20{,}000$ needs 3.2 GB for $\mathbf{K}$ and the same again for its inverse. Above roughly $n = 10{,}000$, switch to $k$-fold with a moving neighbourhood, which caps every solve at the neighbourhood size and makes the whole exercise $O(n k m^3)$ for $m$ neighbours — linear in $n$, and trivially parallel across folds with joblib or a process pool.

Caching. The distance matrix is independent of the variogram, so compute pdist once and reuse it across every candidate model. In a model-selection sweep over a dozen models that alone removes most of the wall time, because evaluating $\gamma$ on a cached distance matrix is a vectorised numpy expression.

Determinism. Fix the seed on any $k$-fold split. An MSDR that moves by 0.05 between runs because the folds were reshuffled will be mistaken for a real difference between models the first time two people compare notes.

Troubleshooting #

Next Steps #

Two child guides take these checks further: Cross-Validating a Variogram Model in Python turns the loop above into a reusable comparison harness over many candidate models, and Diagnosing a Bad Variogram Fit works backwards from each failing statistic to the fitting decision that caused it. If a diagnostic sends you back to the model family rather than the parameters, return to Theoretical Variogram Models; if it sends you back to the binning, return to Empirical Variogram Estimation.

Frequently Asked Questions #

What counts as a good mean squared deviation ratio?

For a correctly scaled variogram the MSDR has expectation one, and with $n$ held-out samples its standard error is roughly $\sqrt{2/n}$. With 240 samples that is about 0.09, so anything between 0.82 and 1.18 is unremarkable. Values above about 1.3 mean the kriging variance is too small and every confidence interval drawn from the model is too narrow; values below about 0.7 mean the opposite. Judge the number against that band, not against an exact one.

Why does the model with the lowest fit RMSE often predict worst?

Fit RMSE measures distance to a dozen or so binned points, each already an average of hundreds of pairs, so it rewards a curve that traces a smoothed summary rather than one that reproduces the field. It weights the far lags as heavily as the near ones, while kriging weights depend almost entirely on the shape near the origin. A Gaussian model is the usual casualty: its parabolic behaviour at short lags fits the plateau beautifully, produces near-singular kriging matrices, and understates the variance badly.

Should I re-fit the variogram inside the cross-validation loop?

Not for the diagnostics described here. Re-fitting inside the loop validates $n$ slightly different models rather than the one you intend to krige with, and the resulting statistics describe an estimator you will never deploy. Fit once on all the data, hold that model fixed, and let only the kriging weights change from pass to pass. Re-fitting inside the loop is appropriate for a different question, namely how stable the fitting procedure is, and should be reported separately as such.

How do I tell a real hole effect from noise?

A genuine hole effect is a dip in the empirical variogram after the first plateau that survives changes to the binning, appears in lags holding hundreds of pairs, and has a physical counterpart such as alternating ore and waste bands or a repeating land-use grid. Re-bin with a different lag count and a different maximum lag; if the dip moves or vanishes it was sampling noise. If it persists, fit a hole-effect model rather than forcing a monotone one through it.

Related

• Cross-Validating a Variogram Model in Python — the leave-one-out loop as a reusable model-comparison harness
• Diagnosing a Bad Variogram Fit — from a failing statistic back to the fitting decision behind it
• Theoretical Variogram Models — the model families a diagnostic may send you back to
• Empirical Variogram Estimation — binning, estimators and pair counts, upstream of every check here
• Uncertainty & Variance Mapping — where a badly scaled variance does its damage

← Back to Variogram Modeling & Semivariance Analysis