1D interpolation

This example builds a one-dimensional noise-free Gaussian process model and selects covariance parameters by restricted maximum likelihood (REML). It is the recommended first complete example because it shows the basic GPmp sequence: construct a model, select covparam, predict, and plot the result.

What this example does

The script creates observation points xi and observations zi from a known reference function. It defines a constant-mean GP with an anisotropic Matern covariance and calls gp.kernel.select_parameters_with_reml. The selected covariance parameters are stored in model.covparam and are then used by model.predict on a dense grid xt.

Mathematical object

The noise-free observation model is

\[Z_i = Z(x_i), \qquad Z \sim \mathcal{GP}(m, k_\theta).\]

The array zi stores the realized values \(z_i\) of \(Z_i\). After selecting \(\theta\), model.predict(xi, zi, xt) computes the conditional Gaussian distribution of \(Z_t=Z(x_t)\) given \(Z_i=z_i\). With the usual block notation, this distribution has mean and covariance

\[m_t + k_{ti} k_{ii}^{-1}(z_i - m_i), \qquad k_{tt} - k_{ti} k_{ii}^{-1} k_{it}.\]

The plotted uncertainty envelope is built from the diagonal of the conditional covariance.

Outputs

The displayed quantities are the reference function, the observations, the posterior mean, and the posterior uncertainty envelope. Because the observations are treated as noise-free, the posterior mean interpolates the observed values. The uncertainty is small near observations and larger away from them.

API points

  • Use gpmp.core.Model or gp.Model to assemble a mean function and a covariance function.

  • Use select_parameters_with_reml when covariance parameters should be selected by REML.

  • Use model.predict(xi, zi, xt) to compute posterior mean and variance at prediction points.

../_images/interpolation_1d_0_0.png

Script: examples/gpmp_example02_1d_interpolation.py

  1"""
  2Plot and optimize the restricted negative log-likelihood
  3
  4Author: Emmanuel Vazquez <emmanuel.vazquez@centralesupelec.fr>
  5Copyright (c) 2022-2026, CentraleSupelec
  6License: GPLv3 (see LICENSE)
  7"""
  8
  9import gpmp.num as gnp
 10import gpmp as gp
 11import gpmp.plot as gpplot
 12
 13
 14def generate_data():
 15    """
 16    Data generation.
 17
 18    Returns
 19    -------
 20    tuple
 21        (xt, zt): target data
 22        (xi, zi): input dataset
 23    """
 24    dim = 1
 25    nt = 200
 26    box = [[-1], [1]]
 27    xt = gp.misc.designs.regulargrid(dim, nt, box)
 28    zt = gp.misc.testfunctions.twobumps(xt)
 29
 30    ni = 6  # FIXME improve stability for ni > 50
 31    xi = gp.misc.designs.ldrandunif(dim, ni, box)
 32    zi = gp.misc.testfunctions.twobumps(xi)
 33
 34    return xt, zt, xi, zi
 35
 36
 37def constant_mean(x, param):
 38    return gnp.ones((x.shape[0], 1))
 39
 40
 41def kernel(x, y, covparam, pairwise=False):
 42    p = 3
 43    return gp.kernel.maternp_covariance(x, y, p, covparam, pairwise)
 44
 45
 46def visualize_results(xt, zt, xi, zi, zpm, zpv):
 47    """
 48    Visualize the results using gp.plot (a matplotlib wrapper).
 49
 50    Parameters
 51    ----------
 52    xt : numpy.ndarray
 53        Target x values
 54    zt : numpy.ndarray
 55        Target z values
 56    xi : numpy.ndarray
 57        Input x values
 58    zi : numpy.ndarray
 59        Input z values
 60    zpm : numpy.ndarray
 61        Posterior mean
 62    zpv : numpy.ndarray
 63        Posterior variance
 64    """
 65    fig = gpplot.Figure(isinteractive=True)
 66    fig.plot(xt, zt, "k", linewidth=1, linestyle=(0, (5, 5)))
 67    fig.plotdata(xi, zi)
 68    fig.plotgp(xt, zpm, zpv, colorscheme="simple")
 69    fig.xylabels("$x$", "$z$")
 70    fig.title("Posterior GP with parameters selected by ReML")
 71    fig.show(grid=True, xlim=[-1.0, 1.0], legend=True, legend_fontsize=9)
 72
 73
 74def main():
 75    xt, zt, xi, zi = generate_data()
 76
 77    model = gp.Model(constant_mean, kernel)
 78
 79    # Automatic selection of parameters using REML
 80    model, info = gp.kernel.select_parameters_with_reml(model, xi, zi, info=True)
 81    gp.modeldiagnosis.diag(model, info, xi, zi)
 82
 83    # Prediction
 84    zpm, zpv = model.predict(xi, zi, xt)
 85
 86    # Visualization
 87    print("\nVisualization")
 88    print("-------------")
 89    plot_likelihood_cross_sections = True
 90    plot_likelihood_2d_profile = False
 91    if plot_likelihood_cross_sections:
 92        gp.modeldiagnosis.plot_selection_criterion_crosssections(
 93            info=info, delta=0.8, param_names=["sigma^2 (log)", "rho (log)"]
 94        )
 95    if plot_likelihood_2d_profile:
 96        gp.modeldiagnosis.plot_selection_criterion_sigma_rho(model, info)
 97
 98    visualize_results(xt, zt, xi, zi, zpm, zpv)
 99
100
101if __name__ == "__main__":
102    main()