Posterior parameter sampling

This example demonstrates REMAP-based parameter selection followed by posterior sampling of GP covariance parameters. It is intended for cases where point estimates of covparam are not enough and uncertainty over covariance parameters should be explored explicitly. Here “posterior” means the unnormalized density obtained by exponentiating the selected criterion. With a REMAP criterion, this is the restricted-likelihood posterior induced by the chosen prior.

What this example does

The rendered preview performs REMAP parameter selection, plots the posterior prediction, runs a short adaptive Metropolis-Hastings sampler, and displays the sampled distribution of the covariance parameters. The full script also shows selection-criterion cross sections and can be switched between Metropolis-Hastings and NUTS. The available public wrappers cover adaptive Metropolis-Hastings, NUTS, and tempered Sequential Monte Carlo (SMC).

Mathematical target

The samplers work on the covariance-parameter vector \(\theta=\mathrm{covparam}\). With the REMAP criterion used in the script, the target is defined by an unnormalized density on \(\theta\). GPmp uses the selection criterion as a negative log density:

\[J(\theta) = J_{\mathrm{REMAP}}(\theta) = -\log \widetilde\pi(\theta\mid z_i),\]

where \(\widetilde\pi\) denotes the target density up to a normalizing constant. Equivalently,

\[\log \widetilde\pi(\theta\mid z_i) = -J_{\mathrm{REMAP}}(\theta).\]

The normalizing constant is not needed by MH, NUTS, or SMC. When a sampling_box is supplied, GPmp truncates the target support: outside the box the log density is \(-\infty\). The init_box controls initialization only. For SMC it is mandatory because the initial particle cloud must be drawn from a user-specified region.

The sampled object is \(\theta\), not a conditional sample path of \(Z\). For each sampled \(\theta\), one may compute a GP conditional distribution for \(Z_t=Z(x_t)\). This separates uncertainty on covariance parameters from conditional GP uncertainty at fixed covariance parameters.

Sampling methods

Metropolis-Hastings

The public MH function uses the lower-level gpmp.mcmc.MetropolisHastings sampler. It runs several chains in parallel. For chain \(c\), the default proposal is a Gaussian random-walk proposal

\[\theta' = \theta_{c,t} + \eta_{c,t}, \qquad \eta_{c,t}\sim\mathcal N(0,\Sigma_c).\]

With the default symmetric proposal, the log acceptance ratio is

\[\log a = \log \widetilde\pi(\theta'\mid z_i) - \log \widetilde\pi(\theta_{c,t}\mid z_i).\]

If symmetric=False is used in the lower-level sampler, GPmp adds the usual reverse-proposal correction from the Metropolis-Hastings ratio [5, 7].

The wrapper only needs criterion values. It builds

\[\log \widetilde\pi_T(\theta\mid z_i) = -J(\theta) / T,\]

where temperature=T. Values \(T>1\) flatten the target and may help exploration. Samples at \(T\ne 1\) are samples from the tempered target, not from the nominal posterior.

The default adaptation method in the public wrapper is Haario adaptation [4]. Every adaptation_interval iterations (50 by default in the wrapper), GPmp computes a block acceptance rate \(r_c\) for each chain. It also computes empirical covariance matrices from the recent chain states. The n_pool argument groups chains before estimating these covariances. If \(\widehat C_g\) is the empirical covariance for the group containing chain \(c\), the proposal covariance is updated as

\[s_c \leftarrow s_c \exp\{\gamma(r_c-r_\star)\}, \qquad \Sigma_c \leftarrow s_c\,\widehat C_g + \varepsilon I.\]

Here \(r_\star\) is the target acceptance rate (0.3 in the public wrapper), \(\gamma\) is the phase-dependent adaptation factor, and \(\varepsilon I\) is a small diagonal stabilization term. The public wrapper currently sets freeze_adaptation=False: adaptation is used during burn-in and continues during the sampling phase, with a smaller sampling-phase adaptation factor. The lower-level MHOptions object can be used directly when a frozen post-burn-in proposal is preferred. The lower-level covariance update also supports the common \(2.38^2/d\) scaling, motivated by optimal-scaling results for random-walk Metropolis algorithms on regular high-dimensional targets [12]. This scaling is a heuristic for practical GP parameter posteriors, not a guarantee.

The lower-level sampler also implements a Robbins-Monro scale update [11], selected with adaptation_method="RM". In that mode, scalar or diagonal proposal scales are multiplied by \(\exp\{\gamma_t(r_c-r_\star)\}\).

The MH object stores the full chain states in mh.x, acceptance indicators in mh.accept, and cached log-target values in mh.log_target_values. The function gpmp.mcmc.get_log_target_values() extracts these values, optionally after burn-in. GPmp evaluates the log target inside each MH step and reuses the cached value of the current state at the next step, because a criterion evaluation can be costly.

NUTS

NUTS is an adaptive Hamiltonian Monte Carlo method [6, 8]. It uses gradients of \(J(\theta)\) to propose distant moves while reducing random-walk behavior. It can be effective when the posterior is smooth and moderately well scaled, but each transition may require many criterion and gradient evaluations.

Sequential Monte Carlo

SMC evolves a population of particles through a sequence of tempered targets [2]. GPmp uses targets proportional to \(\exp[-J(\theta)/T]\), starting from a high temperature and ending at the requested final temperature. Resampling and MH rejuvenation steps are used to keep particles in high-density regions. SMC is useful for comparing separated posterior regions, but it is sensitive to the initialization region and the tempering schedule.

How to interpret the output

REMAP provides a regularized point estimate and a prior definition. Sampling then explores nearby and competing covariance-parameter values according to the same criterion interpreted as a negative log density. Chains or particles concentrated near a small region indicate locally well-identified covariance parameters. Broad, skewed, or separated clouds indicate weak identification or competing covariance explanations.

For MH and NUTS, inspect chain traces, acceptance rates, and log-probability traces before interpreting summaries. For SMC, inspect the final particle cloud, effective sample size behavior, and whether resampling has collapsed particles into a small number of distinct states.

The parameter-distribution figure below uses corner on the transformed coordinates \(\log_{10}(\sigma)\) and \(\log_{10}(\rho)\). The diagonal panels show marginal distributions. The off-diagonal panel shows the two-dimensional sample density. The dark-gold marker and lines show the REMAP estimate used to initialize the posterior exploration.

Functions used

  • select_parameters_sigma2_rho_with_remap_logsigma2_logrho_prior selects a regularized covariance vector and defines criterion callables in info.

  • Posterior samplers in gpmp.mcmc use info.selection_criterion_nograd as a negative log density when possible. NUTS requires the differentiable info.selection_criterion.

  • sampling_box defines hard bounds for the sampling target. init_box is only an initialization region, except that SMC requires it to draw the initial particles.

  • Use sampler diagnostics before interpreting posterior parameter samples.

../_images/posterior_sampling_0_0.png

The corner plot is produced by the following code.

import corner
import numpy as np
from matplotlib.lines import Line2D

log10_sigma = samples[..., 0].reshape(-1) / (2.0 * np.log(10.0))
log10_rho = -samples[..., 1].reshape(-1) / np.log(10.0)
samples_corner = np.column_stack((log10_sigma, log10_rho))
remap = (
    float(info.covparam[0] / (2.0 * np.log(10.0))),
    float(-info.covparam[1] / np.log(10.0)),
)
plot_range = []
for j in range(samples_corner.shape[1]):
    lo = min(np.min(samples_corner[:, j]), remap[j])
    hi = max(np.max(samples_corner[:, j]), remap[j])
    margin = max(0.20 * (hi - lo), 0.05)
    plot_range.append((lo - margin, hi + margin))

fig = corner.corner(
    samples_corner,
    labels=[r"$\log_{10}(\sigma)$", r"$\log_{10}(\rho)$"],
    range=plot_range,
    bins=28,
    smooth=1.0,
    smooth1d=1.0,
    quantiles=[0.16, 0.5, 0.84],
    show_titles=True,
    title_fmt=".2f",
    plot_datapoints=True,
    plot_density=True,
    plot_contours=True,
    fill_contours=True,
    contour_kwargs={"colors": "#003d45", "linewidths": 0.9},
    pcolor_kwargs={"cmap": "cividis"},
    hist_kwargs={"color": "#027bab", "linewidth": 1.2},
)
axes = np.asarray(fig.axes).reshape((2, 2))
remap_color = "#b8860b"
axes[0, 0].axvline(remap[0], color=remap_color, linewidth=1.5)
axes[1, 0].axvline(remap[0], color=remap_color, linewidth=1.5)
axes[1, 0].axhline(remap[1], color=remap_color, linewidth=1.5)
axes[1, 0].plot(
    remap[0],
    remap[1],
    marker="o",
    color=remap_color,
    markersize=5,
)
axes[1, 1].axvline(remap[1], color=remap_color, linewidth=1.5)
axes[0, 1].axis("off")
axes[0, 1].legend(
    [Line2D([0], [0], color=remap_color, marker="o", linewidth=1.5)],
    ["REMAP estimate"],
    loc="center",
    frameon=False,
)
fig.suptitle("Sampled covariance-parameter distribution", y=1.02)
fig.tight_layout()
../_images/posterior_sampling_1_0.png

Script: examples/gpmp_example23_1d_interpolation_posterior_sampling.py

  1"""
  2Demonstrates ReMAP-based GP parameter selection, posterior sampling,
  3and visualization of a 1D Gaussian process model.
  4
  5Author: Emmanuel Vazquez <emmanuel.vazquez@centralesupelec.fr>
  6Copyright (c) 2022-2026, CentraleSupelec
  7License: GPLv3 (see LICENSE)
  8"""
  9
 10import gpmp.num as gnp
 11import gpmp as gp
 12from gpmp.mcmc.param_posterior import (
 13    sample_from_selection_criterion_mh,
 14    sample_from_selection_criterion_nuts,
 15)
 16import matplotlib.pyplot as plt
 17from matplotlib import interactive
 18
 19
 20def generate_data():
 21    """
 22    Data generation.
 23
 24    Returns
 25    -------
 26    tuple
 27        (xt, zt): target data
 28        (xi, zi): input dataset
 29    """
 30    dim = 1
 31    nt = 200
 32    box = [[-1], [1]]
 33    xt = gp.misc.designs.regulargrid(dim, nt, box)
 34    zt = gp.misc.testfunctions.twobumps(xt)
 35
 36    ni = 8
 37    xi = gp.misc.designs.ldrandunif(dim, ni, box)
 38    zi = gp.misc.testfunctions.twobumps(xi)
 39
 40    return xt, zt, xi, zi
 41
 42
 43def constant_mean(x, param):
 44    return gnp.ones((x.shape[0], 1))
 45
 46
 47def kernel(x, y, covparam, pairwise=False):
 48    p = 3
 49    return gp.kernel.maternp_covariance(x, y, p, covparam, pairwise)
 50
 51
 52def visualize_results(xt, zt, xi, zi, zpm, zpv):
 53    """
 54    Visualize the results using gp.plot.plotutils (a matplotlib wrapper).
 55
 56    Parameters
 57    ----------
 58    xt : numpy.ndarray
 59        Target x values
 60    zt : numpy.ndarray
 61        Target z values
 62    xi : numpy.ndarray
 63        Input x values
 64    zi : numpy.ndarray
 65        Input z values
 66    zpm : numpy.ndarray
 67        Posterior mean
 68    zpv : numpy.ndarray
 69        Posterior variance
 70    """
 71    fig = gp.plot.Figure(isinteractive=True)
 72    fig.plot(xt, zt, "k", linewidth=1, linestyle=(0, (5, 5)))
 73    fig.plotdata(xi, zi)
 74    fig.plotgp(xt, zpm, zpv, colorscheme="simple")
 75    fig.xylabels("$x$", "$z$")
 76    fig.title("Posterior GP with parameters selected by ReMAP")
 77    fig.show(grid=True, xlim=[-1.0, 1.0], legend=True, legend_fontsize=9)
 78
 79
 80def main():
 81    xt, zt, xi, zi = generate_data()
 82
 83    model = gp.core.Model(constant_mean, kernel)
 84
 85    # Automatic selection of parameters using ReMAP
 86    model, info = (
 87        gp.kernel.select_parameters_sigma2_rho_with_remap_logsigma2_logrho_prior(
 88            model, xi, zi, info=True
 89        )
 90    )
 91    gp.modeldiagnosis.diag(
 92        model, "linear_mean_maternp_anisotropic", info, xi, zi
 93    )
 94
 95    # Prediction
 96    zpm, zpv = model.predict(xi, zi, xt)
 97
 98    sampler = "nuts"  # "mh" or "nuts"
 99    n_chains = 4
100    nuts_init_box = None
101    if hasattr(info, "bounds") and info.bounds is not None:
102        nuts_init_box = [
103            [b[0] for b in info.bounds],
104            [b[1] for b in info.bounds],
105        ]
106
107    if sampler == "mh":
108        samples, _sampler_state = sample_from_selection_criterion_mh(
109            info,
110            n_steps_total=10_000,
111            burnin_period=5_000,
112            n_chains=n_chains,
113            show_progress=True,
114        )
115    elif sampler == "nuts":
116        samples, _sampler_state = sample_from_selection_criterion_nuts(
117            info,
118            num_samples=500,
119            num_warmup=1_000,
120            n_chains=n_chains,
121            init_box=nuts_init_box,
122            progress=True,
123        )
124    else:
125        raise ValueError("Unknown sampler. Use 'mh' or 'nuts'.")
126
127    # Visualization
128    print("\nVisualization")
129    print("-------------")
130    interactive(True)
131    plot_likelihood_cross_sections = True
132    plot_likelihood_2d_profile = True
133    if plot_likelihood_cross_sections:
134        gp.modeldiagnosis.plot_selection_criterion_crosssections(
135            info=info, delta=0.6, param_names=["log(sigma^2)", "log(1/rho)"]
136        )
137    if plot_likelihood_2d_profile:
138        gp.modeldiagnosis.plot_selection_criterion_sigma_rho(
139            model, info, criterion_name="log posterior"
140        )
141
142    plt.scatter(
143        gnp.log10(gnp.exp(samples[0, :, 0] / 2)),
144        gnp.log10(gnp.exp(-samples[0, :, 1])),
145        alpha=0.2,
146    )
147
148    visualize_results(xt, zt, xi, zi, zpm, zpv)
149
150
151if __name__ == "__main__":
152    main()