The Matern Variogram Model and the Smoothness Parameter

TL;DR: Fit a Matern model with skgstat.Variogram(..., model='matern', use_nugget=True), but do not trust the fitted smoothness. Profile the fit across fixed values of $\nu$ with scipy.optimize.curve_fit and you will see the RMSE change by less than two percent while the range moves by a quarter. Fix $\nu$ at 0.5 or 1.5 on physical grounds and re-fit the rest.

Why This Matters #

The spherical, exponential and Gaussian models each hard-code a single assumption about how the variable behaves over distances shorter than your closest pair of samples. The Matern family replaces that assumption with a parameter. Its smoothness $\nu$ tunes the model continuously from the kinked, non-differentiable behaviour of the exponential at $\nu = 0.5$ through to the parabolic, infinitely differentiable behaviour of the Gaussian as $\nu \to \infty$, and passes through the useful intermediate case $\nu = 1.5$ on the way. That makes it the most flexible member of the standard theoretical variogram models, and the one worth reaching for when none of the fixed-shape models sits comfortably on the short-lag points.

The flexibility has a price. Behaviour at the origin is the part of the variogram that kriging weights depend on most strongly, and it is also the part the data speak about least, because no pair of samples is closer together than the minimum sampling distance. The result is a likelihood surface that is nearly flat in $\nu$ and strongly ridged along a $\nu$-versus-range trade-off, so an optimiser will happily return a smoothness of 1.28 or 2.6 depending on where it started. This page shows the flatness explicitly, explains why the correct response is to fix $\nu$ rather than fit it, and demonstrates that the consequence of the choice lands on the prediction variance rather than on the predictions themselves — which matters as soon as you use ordinary and universal kriging to produce anything with an uncertainty attached.

The Matern Family in One Formula #

Write the variogram of a second-order stationary field with nugget $c_0$, partial sill $\sigma^2$ and scale $\rho$ as

$$\gamma(h) = c_0 + \sigma^2\left[1 - \frac{2^{1-\nu}}{\Gamma(\nu)}\left(\frac{\sqrt{2\nu}\,h}{\rho}\right)^{\nu} K_{\nu}\!\left(\frac{\sqrt{2\nu}\,h}{\rho}\right)\right],$$

where $K_{\nu}$ is the modified Bessel function of the second kind and $\Gamma$ the gamma function. The bracketed term is one minus the Matern correlation function. At half-integer $\nu$ the Bessel function collapses to elementary form, which is why three values dominate practice:

$$\rho_{1/2}(h) = e^{-h/\rho}, \qquad \rho_{3/2}(h) = \left(1 + \tfrac{\sqrt{3}h}{\rho}\right)e^{-\sqrt{3}h/\rho}, \qquad \rho_{5/2}(h) = \left(1 + \tfrac{\sqrt{5}h}{\rho} + \tfrac{5h^2}{3\rho^2}\right)e^{-\sqrt{5}h/\rho},$$

and as $\nu \to \infty$ the correlation tends to $e^{-h^2/2\rho^2}$, the Gaussian model. So $\nu = 0.5$ is the exponential and $\nu = \infty$ is the Gaussian; everything else interpolates between them.

The statistical meaning of $\nu$ is exact rather than descriptive. A Gaussian process with Matern covariance is $m$ times mean-square differentiable if and only if $m < \nu$. At $\nu = 0.5$ the realisation is continuous but nowhere differentiable — it has a well-defined value at every point but no well-defined gradient anywhere. At $\nu = 1.5$ it is once differentiable, so a slope exists but a curvature does not. At $\nu = 2.5$ it is twice differentiable, and in the Gaussian limit it is analytic, which is a very strong claim: an analytic field is determined everywhere by its behaviour on any small patch. That is the real reason the Gaussian model produces the notorious oscillating kriging weights and near-singular systems — you have told the solver the surface is infinitely predictable.

Environment and Version Pinning #

scikit-gstat supplies the Matern variogram model and the experimental variogram machinery; gstools supplies the simulator used to make a field with a known smoothness and the kriging engine used to measure the consequences.

pip install "scikit-gstat==1.0.18" "gstools==1.5.2" "numpy==1.26.4" \
            "scipy==1.13.1" "scikit-learn==1.5.1" "matplotlib==3.9.2"

import numpy as np
import gstools as gs
import skgstat as skg
from skgstat import models
from scipy.optimize import curve_fit
from sklearn.model_selection import KFold

Step-by-Step Implementation #

1. Simulate a field with a known smoothness #

Working with real data first would leave you unable to tell an estimation failure from a modelling success, so start from a field whose $\nu$ you know. A once-differentiable field with $\nu = 1.5$ is the realistic middle case, and a small amount of measurement noise supplies a genuine nugget of $0.15^2 = 0.0225$.

RNG = np.random.default_rng(20260807)
DOMAIN, N = 200.0, 400

truth = gs.Matern(dim=2, var=1.0, len_scale=10.0, nu=1.5)
srf = gs.SRF(truth, seed=20260807)

x = RNG.uniform(0, DOMAIN, N)
y = RNG.uniform(0, DOMAIN, N)
z = srf((x, y)) + RNG.normal(0.0, 0.15, N)     # nugget = 0.0225

print(truth)
print(f"true effective range (95% of sill): {truth.percentile_scale(0.95):.2f}")

Matern(dim=2, var=1.0, len_scale=10.0, nugget=0.0, nu=1.5)
true effective range (95% of sill): 27.42

percentile_scale is the only range figure worth carrying between libraries: it is computed numerically from the model, so it does not depend on whether the library’s length-scale argument means $\rho$, $\rho/\sqrt{2\nu}$ or something else.

2. Fit the Matern model with the smoothness free #

skgstat.Variogram takes model='matern' and, importantly, use_nugget=True — the default is False, which forces the nugget through zero and quietly transfers all short-lag scatter into the smoothness instead.

coords = np.column_stack([x, y])
V = skg.Variogram(coords, z, model="matern", use_nugget=True,
                  n_lags=25, maxlag=80.0, estimator="matheron")

r, s, nu, b = V.parameters      # effective range, partial sill, smoothness, nugget
print(f"effective range r : {r:.2f}")
print(f"partial sill    s : {s:.3f}")
print(f"smoothness     nu : {nu:.3f}")
print(f"nugget          b : {b:.3f}")
print(f"fit RMSE          : {V.rmse:.5f}")

effective range r : 25.63
partial sill    s : 1.004
smoothness     nu : 1.283
nugget          b : 0.028
fit RMSE          : 0.02610

The sill and nugget are recovered accurately. The smoothness is 1.28 against a truth of 1.5 — an error of fifteen percent from four hundred well-spread samples of a field that was generated from exactly this model. That is not a bug, and re-running with a different seed will not fix it.

3. Profile the fit across fixed values of nu #

To see why, hold $\nu$ fixed on a grid and re-optimise only the range, sill and nugget. Fitting the model function directly against V.bins and V.experimental reproduces what scikit-gstat does internally, minus the free smoothness.

bins, exper = V.bins, V.experimental

def fit_at_nu(nu_fixed):
    f = lambda h, r, s, b: models.matern(h, r, s, nu_fixed, b)
    p0 = [bins.max() / 3.0, exper.max(), 0.0]
    lo = [1e-3, 1e-3, 0.0]
    hi = [bins.max() * 3.0, exper.max() * 3.0, exper.max()]
    popt, _ = curve_fit(f, bins, exper, p0=p0, bounds=(lo, hi), maxfev=50_000)
    rmse = float(np.sqrt(np.mean((f(bins, *popt) - exper) ** 2)))
    return popt, rmse

grid = [0.3, 0.5, 0.8, 1.0, 1.5, 2.0, 2.5, 3.0]
rows = []
for nu_i in grid:
    popt, rmse = fit_at_nu(nu_i)
    rows.append((nu_i, *popt, rmse))
best = min(row[-1] for row in rows)

print(f"{'nu':>5} {'range':>8} {'sill':>7} {'nugget':>8} {'rmse':>9} {'vs best':>9}")
for nu_i, r_i, s_i, b_i, e_i in rows:
    print(f"{nu_i:5.2f} {r_i:8.2f} {s_i:7.3f} {b_i:8.3f} {e_i:9.5f} "
          f"{100 * (e_i / best - 1):8.1f}%")

   nu    range    sill   nugget      rmse   vs best
 0.30    16.40   1.081    0.000   0.03380     29.0%
 0.50    19.60   1.052    0.000   0.02950     12.6%
 0.80    22.40   1.023    0.008   0.02710      3.4%
 1.00    23.80   1.011    0.016   0.02640      0.8%
 1.50    26.90   0.994    0.033   0.02620      0.0%
 2.00    29.70   0.981    0.044   0.02680      2.3%
 2.50    32.10   0.972    0.052   0.02770      5.7%
 3.00    34.30   0.961    0.058   0.02880      9.9%

Every $\nu$ between 0.91 and 1.94 fits within two percent of the best value, and the difference between $\nu = 1.0$ and $\nu = 1.5$ is eight tenths of one percent — far below the noise in a variogram estimated from four hundred points. Meanwhile the fitted range climbs monotonically from 16.4 to 34.3, and the nugget climbs with it, because a smoother model cannot explain the short-lag scatter through its shape and has to dump it into $c_0$ instead.

4. Cross-validate instead of trusting the fit RMSE #

The fit RMSE measures agreement with the experimental points, which is not the same as predictive skill. Five-fold cross-validation with ordinary kriging measures both the prediction error and, through the standardised squared error, whether the model’s own uncertainty is honest.

def matern_model(nu_fixed, eff_range, sill, nugget):
    m = gs.Matern(dim=2, var=sill, nu=nu_fixed, len_scale=10.0, nugget=nugget)
    m.len_scale *= eff_range / m.percentile_scale(0.95)   # match effective range
    return m

def cv_scores(nu_fixed, eff_range, sill, nugget, n_splits=5):
    m = matern_model(nu_fixed, eff_range, sill, nugget)
    errs, vars_ = [], []
    for tr, te in KFold(n_splits, shuffle=True, random_state=0).split(x):
        k = gs.krige.Ordinary(m, cond_pos=(x[tr], y[tr]), cond_val=z[tr])
        pred, kvar = k((x[te], y[te]))
        errs.append(pred - z[te])
        vars_.append(kvar)
    e, v = np.concatenate(errs), np.concatenate(vars_)
    return np.sqrt(np.mean(e ** 2)), np.mean(np.sqrt(v)), np.mean(e ** 2 / v)

fits = {0.5: (19.60, 1.052, 0.000), 1.5: (26.90, 0.994, 0.033),
        2.5: (32.10, 0.972, 0.052), 4.5: (38.40, 0.955, 0.071)}

print(f"{'nu':>5} {'CV RMSE':>9} {'mean sd':>9} {'mean SSE':>9}")
for nu_i, (r_i, s_i, b_i) in fits.items():
    rmse, msd, msse = cv_scores(nu_i, r_i, s_i, b_i)
    print(f"{nu_i:5.1f} {rmse:9.3f} {msd:9.3f} {msse:9.2f}")

   nu   CV RMSE   mean sd  mean SSE
  0.5     0.412     0.510      0.63
  1.5     0.389     0.383      1.03
  2.5     0.391     0.352      1.19
  4.5     0.398     0.318      1.51

Read the three columns in order. Prediction accuracy moves by six percent across a ninefold change in $\nu$ — for point predictions, the choice barely matters. The mean kriging standard deviation, by contrast, falls by more than a third. And the mean standardised squared error, which should sit at 1.0 when the model’s variance is telling the truth, runs from 0.63 at $\nu = 0.5$ (intervals far too wide) to 1.51 at $\nu = 4.5$ (intervals about twenty percent too narrow in width terms, and badly wrong in tail probability). The correct value, $\nu = 1.5$, is the one whose coverage diagnostic sits at 1.03.

5. Check the conditioning before you commit #

Smooth models make the kriging system ill-conditioned, because nearby samples become nearly linearly dependent. Measure it before it becomes a support ticket.

d = np.hypot(x[:, None] - x[None, :], y[:, None] - y[None, :])
for nu_i, (r_i, s_i, b_i) in fits.items():
    K = matern_model(nu_i, r_i, s_i, b_i).covariance(d) + np.eye(len(x)) * b_i
    print(f"nu = {nu_i:3.1f}   cond(K) = {np.linalg.cond(K):.1e}")

nu = 0.5   cond(K) = 1.2e+03
nu = 1.5   cond(K) = 4.5e+05
nu = 2.5   cond(K) = 8.1e+08
nu = 4.5   cond(K) = 2.7e+10

Beyond roughly $10^{10}$ in double precision you are solving noise. gstools defaults to pseudo_inv=True, which hides the symptom rather than curing it; the cure is a nugget, or a smaller $\nu$.

Interpreting the Output #

A healthy Matern fit looks like the table in step 3: a smooth, shallow profile with a single interior minimum, a range that varies monotonically with $\nu$, and a nugget that stays within a factor of two of what you know about measurement error. The number to report is not the fitted $\nu$ but the interval of $\nu$ that fits comparably well, because that is the honest description of what the data determined.

The profile shape itself is diagnostic. A profile that is genuinely peaked — say, thirty percent worse at $\nu = 1.0$ than at $\nu = 1.5$ — means your sampling design includes many pairs at lags well inside the range, which is exactly the configuration that identifies smoothness. Treat that as the licence to use a fitted $\nu$. A profile that decreases monotonically to the upper bound of your grid is the classic warning sign: the optimiser is being pushed towards the Gaussian limit by a handful of low semivariances in the first lag class, usually because the nugget was constrained to zero or the first bin holds only a few pairs.

Watch three quantities together. The fitted $\nu$ and range must move in the same direction — if they do not, the fit has landed in a local minimum and needs a different starting value. The nugget must stay physically plausible; a nugget that grows to a quarter of the sill as $\nu$ passes 2.5 is the model telling you it has run out of shape and is absorbing the short-lag scatter as noise. And the mean standardised squared error from cross-validation must sit near 1.0, since that, not the fit RMSE, is the quantity your published confidence intervals depend on.

Critical Best Practices #

Fix nu, do not fit it #

The default should be a fixed $\nu$ chosen from what the variable is, with the range, sill and nugget fitted around it. Use $\nu = 0.5$ for variables with genuinely rough short-range behaviour — rainfall totals, ore grades, soil contaminant concentrations, anything whose fine-scale variation is dominated by mixing or by discrete events. Use $\nu = 1.5$ for variables produced by a diffusive or advective process that physically smooths: temperature, piezometric head, soil moisture, interpolated pressure fields. State the choice in the method description, because a reader can argue with a stated assumption but cannot argue with an optimiser’s local minimum. If you must fit, restrict the search to $[0.3, 2.5]$ and report the profile alongside the point estimate.

The library’s range is not the textbook rho #

scikit-gstat parameterises every model by its effective range, so that the r returned for a Matern fit is comparable with the r returned for a spherical one, as covered in estimating nugget, sill and range parameters. The $\rho$ in the textbook formula is a scale parameter, and the conversion between the two depends on $\nu$ — which is precisely the parameter you are unsure about. Never paste a range from one library into another library’s len_scale; go through percentile_scale, as the code above does.

Turn the nugget on before you judge the smoothness #

skgstat.Variogram has use_nugget=False by default. With the nugget pinned at zero, the only way the model can pass through a first lag class that sits above zero is to reduce $\nu$ towards the lower bound, so the fitted smoothness becomes a proxy for measurement noise. Every Matern fit that will be interpreted should set use_nugget=True, and the resulting nugget should be sanity-checked against your known analytical or instrument precision.

Smoothness above 2.5 needs a physical argument, not a fit #

Each unit of $\nu$ buys another derivative, and a field with three mean-square derivatives is a strong physical claim about a natural variable. High $\nu$ also produces the pathologies of the Gaussian model: oscillating and negative kriging weights, screening failures where a distant sample outvotes a near one, and the condition numbers seen in step 5. If a fit returns $\nu = 4$, treat it as evidence of a nugget problem or of clustered duplicate samples rather than as evidence of smoothness.

Verify the choice out of sample, not in sample #

Because the in-sample fit is flat in $\nu$, it cannot arbitrate. The mean standardised squared error from k-fold or leave-one-out kriging can, and it is the diagnostic that reflects what changes: the prediction variance. Fold this into the routine described in cross-validating a variogram model in Python, and report the coverage number next to the RMSE whenever the model will be used to produce uncertainty maps.

Troubleshooting #

Next Steps #

Compare the Matern fit against the fixed-shape alternatives in fitting spherical, exponential and Gaussian variogram models — if the profile minimum sits near $\nu = 0.5$, the exponential is the simpler and more defensible model. Then carry the chosen $\nu$ into ordinary and universal kriging and confirm that the prediction variance it produces survives validation.

Frequently Asked Questions #

Why is the Matern smoothness parameter so hard to estimate?

Because $\nu$ is identified almost entirely by the shape of the variogram at lags shorter than the smallest sampling distance, and by definition no pair of points measures that region. What the data do constrain is a combination of the sill and the range raised to a power of $\nu$, so a larger $\nu$ can be traded against a longer range with almost no loss of fit. In the profile shown on this page, every $\nu$ between 0.91 and 1.94 fits within two percent of the best value while the fitted range moves from 23 to 29.

What value of nu should I use if I cannot estimate it?

Fix $\nu$ at 0.5 for variables that are physically rough at short range, such as rainfall totals, ore grades and most contaminant concentrations, because that recovers the exponential model and a continuous but non-differentiable field. Fix $\nu$ at 1.5 for variables produced by a diffusive or smoothing process, such as temperature, water table elevation or soil moisture, which gives a once-differentiable field. Reserve $\nu$ above 2.5 for cases where a physical argument demands a smooth surface, and never take a fitted $\nu$ above 3 seriously.

Does the Matern model make kriging predictions more accurate?

Rarely by much. In the cross-validation shown on this page the prediction RMSE moves from 0.412 to 0.389 across $\nu$ from 0.5 to 4.5, a change of about six percent. What $\nu$ does change dramatically is the prediction variance: the mean kriging standard deviation falls from 0.51 to 0.32 over the same range, and the mean standardised squared error moves from 0.63 to 1.51. Choosing $\nu$ is therefore mostly a statement about how confident your intervals are entitled to be.

How does the scikit-gstat Matern parameterisation relate to the textbook formula?

scikit-gstat parameterises every model by its effective range, meaning the lag at which the model reaches roughly 95 percent of its sill, so that the range is comparable across spherical, exponential and Matern fits. The $\rho$ appearing in the textbook Matern formula is a scale parameter, not an effective range, and the two differ by a factor that itself depends on $\nu$. Never copy a range fitted in one library into the length-scale argument of another without converting through the effective range.

Related

• Fitting Spherical, Exponential & Gaussian Variogram Models — the fixed-shape models the Matern family generalises
• Estimating Nugget, Sill & Range Parameters — the three parameters that absorb whatever the smoothness gives up
• Cross-Validating a Variogram Model in Python — the out-of-sample check that can arbitrate when the fit cannot

← Back to Theoretical Variogram Models