Skip to content

Posterior Sampling

Condition GP samples on observed data.

Basic Posterior

Observe f at sparse points, predict at dense points:

import numpy as np
from gp4c import sample_posterior, Observations, SamplingSpec

# Training data
x_train = np.array([0.0, 1.0, 2.0, 3.0, 4.0])
y_train = np.sin(x_train)
obs = Observations(x_f=x_train, y_f=y_train, noise_f=0.01)

# Prediction points
x_test = np.linspace(0, 5, 100)
spec = SamplingSpec(x_f=x_test)

# Sample posterior
result = sample_posterior(obs, spec, ell=1.0, n_samples=10)

print(result.f_mean.shape)  # (100,) - posterior mean
print(result.f.shape)       # (10, 100) - posterior samples

The posterior mean result.f_mean interpolates through the training points.

Observation Noise

Specify noise level for each observation type:

obs = Observations(
    x_f=x_train,
    y_f=y_train,
    noise_f=0.1  # 10% noise on function observations
)

Higher noise allows more flexibility around observations. Lower noise forces posterior to pass closer to data.

Derivative Prediction

Train on f, predict the derivative h:

# Observe function
obs = Observations(x_f=x_train, y_f=y_train, noise_f=0.01)

# Predict derivative
spec = SamplingSpec(x_h=x_test)
result = sample_posterior(obs, spec, ell=1.0, n_samples=10)

print(result.h_mean.shape)  # (100,) - predicted derivative mean
print(result.h.shape)       # (10, 100) - derivative samples

This is useful for sensitivity analysis and gradient estimation.

Inverse Problems

Recover f from Derivatives

Observe derivatives, infer the function:

# Observe derivative data
x_deriv = np.linspace(0, 5, 20)
y_deriv = np.cos(x_deriv)  # derivative of sin
obs = Observations(x_h=x_deriv, y_h=y_deriv, noise_h=0.01)

# Infer function
x_dense = np.linspace(0, 5, 100)
spec = SamplingSpec(x_f=x_dense)
result = sample_posterior(obs, spec, ell=1.0, n_samples=10)

The recovered function captures the integral structure consistent with observed derivatives.

Recover f from Integrals

Know the integral, infer the function:

# Observe integral values
x_integral = np.array([1.0, 2.0, 3.0, 4.0])
y_integral = np.array([0.5, 1.2, 2.0, 2.5])
obs = Observations(x_g=x_integral, y_g=y_integral, noise_g=0.01)

# Infer function
spec = SamplingSpec(x_f=x_dense)
result = sample_posterior(obs, spec, ell=1.0, n_samples=10)

Gradient-Enhanced GP

Use both function and derivative observations:

obs = Observations(
    x_f=x_func, y_f=y_func, noise_f=0.01,
    x_h=x_deriv, y_h=y_deriv, noise_h=0.02
)
spec = SamplingSpec(x_f=x_test)
result = sample_posterior(obs, spec, ell=1.0, n_samples=10)

Combining both observation types improves interpolation quality, especially in regions between training points.

Mixed Observations

Observe multiple quantities simultaneously:

obs = Observations(
    x_f=x_func, y_f=y_func, noise_f=0.01,
    x_g=x_integral, y_g=y_integral, noise_g=0.05,
    x_h=x_deriv, y_h=y_deriv, noise_h=0.02
)

# Predict all three
spec = SamplingSpec(x_f=x_test, x_g=x_test, x_h=x_test)
result = sample_posterior(obs, spec, ell=1.0, n_samples=10)

Warning

At least one observation type must be provided. Each observation needs both x and y arrays.

Posterior Uncertainty

Posterior samples quantify predictive uncertainty:

result = sample_posterior(obs, spec, ell=1.0, n_samples=100)

# Compute credible intervals
f_std = np.std(result.f, axis=0)
f_lower = result.f_mean - 2 * f_std
f_upper = result.f_mean + 2 * f_std

Uncertainty is small near observations and increases away from data.

Next Steps