2D interpolation

This example builds an anisotropic Matérn GP on a two-dimensional test function. It illustrates the same sequence as the 1D interpolation example, but with a two-dimensional input space and contour plots for spatial diagnostics.

What this example does

The script selects a two-dimensional test function, draws observation points, defines a GP, and selects covariance parameters with a REMAP criterion. The fixed-regularity Matérn covariance uses the sigma2_rho layout [log(sigma2), -log(rho_0), -log(rho_1)]. The selected predictor is evaluated on a regular grid to compare the reference function and the GP prediction.

Mathematical description

The GP uses the same noise-free conditional formulation as in the 1D interpolation example, but with two-dimensional points \(x=(x_0,x_1)\). The anisotropic Matérn kernel uses the scaled distance

\[h(x,x') = \left[ \left(\frac{x_0-x'_0}{\rho_0}\right)^2 + \left(\frac{x_1-x'_1}{\rho_1}\right)^2 \right]^{1/2}.\]

The two lengthscales control smoothing along the two coordinate axes. Small \(\rho_j\) allows faster variation along coordinate \(j\). Large \(\rho_j\) makes the posterior vary more slowly along that coordinate.

Outputs

The four panels contain the reference function, GP posterior mean, absolute prediction error, and posterior standard deviation. Red dots mark observation points. The error and posterior standard deviation panels should be read together: a well-calibrated model should assign larger uncertainty in regions with sparse observations or difficult extrapolation.

Functions used

  • gp.misc.designs.ldrandunif creates a space-filling observation design.

  • gp.kernel.select_parameters_sigma2_rho_with_remap_logsigma2_logrho_prior selects covariance parameters with the default REMAP criterion.

  • model.predict returns posterior means and variances at the grid points.

  • Convert backend arrays with gnp.to_np before sending them to Matplotlib functions that expect NumPy arrays.

../_images/interpolation_2d_0_0.png

Script: examples/gpmp_example03_2d.py

  1"""
  2Gaussian processes in 2D
  3
  4An anisotropic Matern covariance function is used for the Gaussian
  5Process (GP) prior. The parameters of this covariance function
  6(variance and ranges) are estimated using the Restricted Maximum
  7A Posteriori (ReMAP) method.
  8
  9The mean function of the GP prior is assumed to be constant and
 10unknown.
 11
 12The function is sampled on a space-filling Latin Hypercube design, and
 13the data is assumed to be noiseless.
 14
 15----
 16Author: Emmanuel Vazquez <emmanuel.vazquez@centralesupelec.fr>
 17Copyright (c) 2022-2026, CentraleSupelec
 18License: GPLv3 (see LICENSE)
 19----
 20This example is based on the file stk_example_kb03.m from the STK at
 21https://github.com/stk-kriging/stk/ 
 22by Julien Bect and Emmanuel Vazquez, released under the GPLv3 license.
 23
 24Original copyright notice:
 25
 26   Copyright (c) 2015, 2016, 2018 CentraleSupelec
 27   Copyright (c) 2011-2014 SUPELEC
 28----
 29"""
 30
 31import numpy as np
 32import gpmp.num as gnp
 33import gpmp as gp
 34import matplotlib.pyplot as plt
 35
 36
 37# Test function selection
 38def select_test_function(case_num):
 39    if case_num == 1:
 40        f = gp.misc.testfunctions.braninhoo
 41        dim = 2
 42        box = [[-5, 0], [10, 15]]
 43        ni = 20
 44    elif case_num == 2:
 45        f = gp.misc.testfunctions.wave
 46        dim = 2
 47        box = [[-1, -1], [1, 1]]
 48        ni = 40
 49    return f, dim, box, ni
 50
 51
 52def create_model():
 53    def constant_mean(x, param):
 54        return gnp.ones((x.shape[0], 1))
 55
 56    def kernel(x, y, covparam, pairwise=False):
 57        p = 6
 58        return gp.kernel.maternp_covariance(x, y, p, covparam, pairwise)
 59
 60    return gp.Model(constant_mean, kernel)
 61
 62
 63def main():
 64    case_num = 1
 65    f, dim, box, ni = select_test_function(case_num)
 66
 67    # Compute the function on a 80 x 80 regular grid
 68    nt = [80, 80]
 69    xt = gp.misc.designs.regulargrid(dim, nt, box)
 70    zt = f(xt)
 71
 72    design_type = "ld"
 73    if design_type == "lhs":
 74        xi = gp.misc.designs.maximinlhs(dim, ni, box)
 75    elif design_type == "ld":
 76        xi = gp.misc.designs.ldrandunif(dim, ni, box)
 77    zi = f(xi)
 78
 79    model = create_model()
 80
 81    # Parameter selection
 82    model, info = gp.kernel.select_parameters_sigma2_rho_with_remap_logsigma2_logrho_prior(
 83        model, xi, zi, info=True
 84    )
 85    gp.modeldiagnosis.diag(
 86        model, "linear_mean_maternp_anisotropic", info, xi, zi
 87    )
 88
 89    # Prediction
 90    (zpm, zpv) = model.predict(xi, zi, xt)
 91
 92    # Visualization
 93    contour_lines = 30
 94    xt_np = gnp.to_np(xt)
 95    xi_np = gnp.to_np(xi)
 96    zt_np = gnp.to_np(zt).reshape(nt)
 97    zpm_np = gnp.to_np(zpm).reshape(nt)
 98    zsd_np = np.sqrt(np.maximum(gnp.to_np(zpv).reshape(nt), 0.0))
 99
100    fig, axes = plt.subplots(nrows=2, ncols=2)
101    data = [zt_np, zpm_np, np.abs(zpm_np - zt_np), zsd_np]
102    titles = [
103        "function to be approximated",
104        f"approximation from {ni} points",
105        "true approx error",
106        "posterior std",
107    ]
108    cmaps = ["PiYG", "PiYG", "magma_r", "viridis"]
109
110    for ax, z, title, cmap in zip(axes.flat, data, titles, cmaps):
111        cs = ax.contourf(
112            xt_np[:, 0].reshape(nt),
113            xt_np[:, 1].reshape(nt),
114            z,
115            levels=contour_lines,
116            cmap=cmap,
117        )
118        ax.plot(xi_np[:, 0], xi_np[:, 1], "ro", label="data")
119        ax.set_title(title)
120        ax.set_xlabel("$x_1$")
121        ax.set_ylabel("$x_2$")
122        ax.legend()
123        fig.colorbar(cs, ax=ax, shrink=0.9)
124
125    plt.show()
126
127    # Predictions vs truth
128    plt.figure()
129    plt.plot(zt, zpm, "ko")
130    (xmin, xmax), (ymin, ymax) = plt.xlim(), plt.ylim()
131    xmin = min(xmin, ymin)
132    xmax = max(xmax, ymax)
133    plt.plot([xmin, xmax], [xmin, xmax], "--")
134    plt.xlabel("true values")
135    plt.ylabel("predictions")
136    plt.show()
137
138    # LOO predictions
139    zloom, zloov, eloo = model.loo(xi, zi)
140    gp.plot.plot_loo(zi, zloom, zloov)
141
142    gp.plot.crosssections(model, xi, zi, box, ind_i=[0, 10], ind_dim=[0, 1])
143
144
145if __name__ == "__main__":
146    main()