Example 06: fixed and selected Matérn regularity

Complete script: examples/example06_matern.py

Purpose

The script compares two REML models on the same one-dimensional observations. Model_ConstantMean_Maternp_REML fixes \(\nu=p+1/2\). Model_ConstantMean_Matern_REML selects \(\nu\) with the process variance and lengthscale. The role of Matérn regularity is discussed by Stein [9] and Chilès and Delfiner [2].

Parameter vectors

For an input dimension \(d\), the fixed-regularity model uses

[log(sigma2), -log(rho_0), ..., -log(rho_{d-1})]

The selected-regularity model adds log(nu):

[log(sigma2), log(nu), -log(rho_0), ..., -log(rho_{d-1})]

The logarithmic coordinate enforces \(\nu>0\). The Param object names this coordinate nu and stores it at path covparam/regularity.

What is computed

• sixteen observations of a sine function with a local absolute-value term.

• one fixed model with p=2, hence \(\nu=2.5\).

• one model with \(\nu\) selected by REML.

• posterior means and variances on a regular grid.

• the selected regularity read through model[0].get_param().

Outputs

The generated run prints the fixed and selected regularities followed by the free-regularity Param object:

Fixed Matérn regularity: nu = 2.5
Selected Matérn regularity: nu = 3.48996
     Name:                     Path       Norm         Bounds     Value    Denorm
z0_sigma2:       covparam->variance        log    (-inf, inf)     2.155     8.624
    z0_nu:     covparam->regularity        log    (-inf, inf)     1.250     3.490
 z0_rho_0:    covparam->lengthscale    log_inv    (-inf, inf)    0.0298     0.971

The numerical value selected for \(\nu\) depends on the observations, criterion, bounds, and numerical backend. The parameter name, path, and positive denormalized value are stable parts of the interface.

Both panels use the same observations. The upper model fixes \(\nu=2.5\); the lower model selects \(\nu\). Each shaded band is a pointwise 95 percent posterior interval. Differences between the panels come from selecting the regularity with the other covariance parameters.

Run python examples/example06_matern.py from the repository root to execute the script. Regenerate the documented output and figure with cd docs && python make_example_results.py.

Source excerpt

def main():
    xi = gnp.linspace(0.0, 1.0, 16).reshape(-1, 1)
    zi = gnp.sin(2.0 * gnp.pi * xi[:, 0]) + 0.3 * gnp.abs(xi[:, 0] - 0.45)
    xt = gnp.linspace(0.0, 1.0, 100).reshape(-1, 1)
    p = 2

    maternp_model = gpc.Model_ConstantMean_Maternp_REML(
        "fixed_regularity",
        1,
        mean_specification={"type": "constant"},
        covariance_specification={"p": p},
    )
    matern_model = gpc.Model_ConstantMean_Matern_REML(
        "selected_regularity",
        1,
        mean_specification={"type": "constant"},
    )

    maternp_model.select_params(xi, zi)
    matern_model.select_params(xi, zi)

    maternp_model.predict(xi, zi, xt)
    matern_model.predict(xi, zi, xt)
    matern_model.run_diagnosis(xi, zi)

    param = matern_model[0].get_param()
    lognu = param.get_by_path(["covparam", "regularity"])
    nu = gnp.to_scalar(gnp.exp(gnp.asarray(lognu).reshape(-1)[0]))
    print(f"Fixed Matérn regularity: nu = {p + 0.5:g}")
    print(f"Selected Matérn regularity: nu = {nu:.6g}")


if __name__ == "__main__":
    main()