Fitting Spherical, Exponential & Gaussian Variogram Models

TL;DR: Build one Variogram with a fixed binning, then reassign V.model to "spherical", "exponential", and "gaussian" in a loop, calling V.fit() and reading V.parameters and V.rmse each time. Rank the models by V.rmse, then confirm the winner’s shape near the origin matches how smooth your variable really is before committing.

Why This Matters #

The experimental variogram is a scatter of points; a theoretical model turns that scatter into a smooth, positive-definite function that kriging can invert. Choosing among the spherical, exponential, and Gaussian families is the core task of working with theoretical variogram models, and it sits at the heart of the wider variogram modeling and semivariance analysis workflow. The three families differ mostly in how they approach the sill and how they behave near the origin, and that near-origin behaviour controls how smooth your kriged surface will be.

The three model functions, each with sill $c$, nugget $c_0$, and effective range $a$, are:

$$\gamma_{sph}(h) = c_0 + c\left(\frac{3h}{2a} - \frac{h^3}{2a^3}\right), \quad h \le a$$

$$\gamma_{exp}(h) = c_0 + c\left(1 - e^{-3h/a}\right)$$

$$\gamma_{gau}(h) = c_0 + c\left(1 - e^{-3h^2/a^2}\right)$$

The spherical model reaches the sill exactly at the range; the exponential and Gaussian models approach it asymptotically, and scikit-gstat rescales each so that $a$ is the effective range where the curve reaches about 95 percent of the sill.

Environment #

pip install "scikit-gstat>=1.0" "numpy>=1.23" "matplotlib>=3.7"

import numpy as np
import matplotlib.pyplot as plt
from skgstat import Variogram

Step-by-Step Implementation #

Step 1 — Build the experimental variogram once #

rng = np.random.default_rng(19)
coords = rng.uniform(0, 500, size=(200, 2))

# Smooth medium-range structure plus small nugget-scale noise
vals = (
    20.0
    + 6.0 * np.exp(-((coords[:, 0] - 250) ** 2 + (coords[:, 1] - 250) ** 2) / 2.0e4)
    + rng.normal(0, 0.6, size=len(coords))
)

V = Variogram(coords, vals, n_lags=14, maxlag="median", use_nugget=True)

Fix the binning before comparing models so every candidate is fitted to the identical experimental cloud. use_nugget=True lets each model fit its own nugget rather than forcing it through zero.

Step 2 — Fit each model and record the fit #

results = {}
for name in ("spherical", "exponential", "gaussian"):
    V.model = name          # reassigning the model invalidates the cached fit
    V.fit()                 # explicit re-fit (also happens lazily on access)
    rng_, sill, nugget = V.parameters
    results[name] = {
        "range": rng_,
        "sill": sill,
        "nugget": nugget,
        "rmse": V.rmse,
    }

Reassigning V.model clears the fit cache, and V.fit() re-runs the least-squares optimisation against the same experimental points. V.parameters and V.rmse then describe that specific model.

Step 3 — Rank the candidates by RMSE #

for name, r in sorted(results.items(), key=lambda kv: kv[1]["rmse"]):
    print(f"{name:12s} range={r['range']:7.1f}  sill={r['sill']:6.2f}  "
          f"nugget={r['nugget']:5.2f}  rmse={r['rmse']:.4f}")

best = min(results, key=lambda k: results[k]["rmse"])
print("lowest RMSE:", best)

The lowest rmse marks the curve that tracks the experimental points most closely. Keep the full table, since a near-tie on RMSE is common and the decision then rests on process knowledge.

Step 4 — Inspect the shapes side by side #

fig, axes = plt.subplots(1, 3, figsize=(15, 4), sharey=True)
for ax, name in zip(axes, ("spherical", "exponential", "gaussian")):
    V.model = name
    V.plot(axes=ax, show=False)
    ax.set_title(f"{name}  (rmse={results[name]['rmse']:.3f})")
fig.savefig("model_comparison.png", dpi=150, bbox_inches="tight")

Look at the origin: the Gaussian curve leaves the nugget flat then rises steeply (parabolic), while the exponential rises fastest immediately. That origin shape, not just the RMSE, tells you which model respects your variable’s smoothness.

Step 5 — Commit the chosen model #

V.model = best
final_range, final_sill, final_nugget = V.parameters
print(f"selected model  : {best}")
print(f"effective range : {final_range:.1f}")
print(f"sill            : {final_sill:.2f}")
print(f"nugget          : {final_nugget:.2f}")

Interpreting the Output #

Each model tells a different story about spatial continuity. The spherical model, with its linear rise near the origin and a hard sill at the range, is the workhorse for many geological and environmental variables. The exponential model rises fastest at short lags and never quite flattens, suiting rough, irregular phenomena. The Gaussian model’s flat-then-steep origin implies extreme short-range continuity, appropriate only for genuinely smooth fields such as elevation or a modelled potential surface.

A lower V.rmse is better, but the Gaussian model frequently wins on RMSE while implying a physically implausible smoothness that produces oscillating, numerically fragile kriging weights. Treat RMSE as the first filter and the origin shape as the deciding factor.

Critical Best Practices #

Fit all models to identical bins #

Comparing RMSE only makes sense if every model sees the same experimental points. Build the Variogram once and reassign V.model; do not rebuild with different n_lags or maxlag between candidates, or you are comparing binnings rather than models.

Distrust a Gaussian that wins by a hair #

The Gaussian model’s parabolic origin often shaves RMSE while encoding unrealistic continuity. Unless your variable is truly smooth, prefer the spherical or exponential fit even at a marginally higher RMSE, and always add a small nugget to stabilize the kriging system.

Let each model fit its own nugget #

Set use_nugget=True so the nugget is a free parameter. Forcing a zero nugget distorts the fitted range and sill and hides measurement error that belongs in the model. The nugget-to-sill ratio is itself a diagnostic, covered in estimating nugget, sill and range parameters.

Cross-check RMSE with leave-one-out kriging #

RMSE against the experimental variogram measures curve fit, not prediction quality. Before locking a model, validate it with leave-one-out or spatial cross-validation so the choice is judged on how well it predicts held-out samples, following the principles in cross-validation strategies.

Troubleshooting #

Next Steps #

After picking a model, interpret and constrain its coefficients in estimating nugget, sill and range parameters, and revisit the estimation of the cloud these models fit in fitting empirical variograms with SciKit-GStat.

Related

• Theoretical Variogram Models — the model families and their properties
• Estimating Nugget, Sill & Range Parameters — reading and constraining the fitted coefficients
• Fitting Empirical Variograms with SciKit-GStat — building the experimental cloud these models fit

← Back to Theoretical Variogram Models

Frequently Asked Questions #

How do I compare variogram models in scikit-gstat?

Fit each candidate to the same experimental variogram by reassigning V.model, then compare V.rmse across models. The lowest RMSE indicates the closest fit to the experimental points, but you should also confirm the origin behaviour matches the physical smoothness of the process.

When should I prefer an exponential over a Gaussian variogram model?

Use the exponential model when the phenomenon varies roughly and reaches its sill gradually, such as rainfall or ore grades. Reserve the Gaussian model for very smooth, continuous fields, because its parabolic behaviour near the origin implies extreme short-range continuity that few natural variables exhibit and that can make kriging numerically unstable.

Is the lowest RMSE always the best variogram model?

Not automatically. RMSE ranks how closely each curve tracks the experimental points, but the model must also respect the process: the Gaussian model can win on RMSE yet imply implausible smoothness. Balance the numeric fit against the near-origin shape and domain knowledge.