Skip to content

Python Examples

Code examples from the examples/ directory.

Basic Usage

File: examples/basic_usage.py

Sample f and its integral g:

import numpy as np
import matplotlib.pyplot as plt
from gp4c import sample_prior, SamplingSpec

# Define grid
x = np.linspace(0, 5, 100)

# Sample joint GP
spec = SamplingSpec(x_f=x, x_g=x)
result = sample_prior(
    spec,
    sigma2=1.0,
    ell=0.5,
    n_samples=5,
    seed=42
)

# Plot results
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))

ax1.plot(x, result.f.T, alpha=0.7)
ax1.set_title('f(x) samples')
ax1.grid(True)

ax2.plot(x, result.g.T, alpha=0.7)
ax2.set_title('g(x) = integral of f')
ax2.grid(True)

plt.tight_layout()
plt.show()

Run: 

python examples/basic_usage.py

Verification

File: examples/verification.py

Verify that g matches numerical integration of f:

import numpy as np
from gp4c import sample_prior, SamplingSpec
from scipy.integrate import cumulative_trapezoid

x = np.linspace(0, 5, 200)
spec = SamplingSpec(x_f=x, x_g=x)
result = sample_prior(spec, ell=0.5, n_samples=1)

# Numerically integrate f
f = result.f[0]
g_numerical = cumulative_trapezoid(f, x, initial=0)
g_analytical = result.g[0]

# Compare
error = np.abs(g_analytical - g_numerical)
print(f"Max absolute error: {error.max():.6f}")
print(f"Mean absolute error: {error.mean():.6f}")

# Plot comparison
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 6))
plt.subplot(2, 1, 1)
plt.plot(x, f, label='f(x)')
plt.legend()
plt.grid(True)

plt.subplot(2, 1, 2)
plt.plot(x, g_analytical, label='g(x) analytical')
plt.plot(x, g_numerical, '--', label='g(x) numerical')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

Run: 

python examples/verification.py

Parameter Study

File: examples/parameter_study.py

Explore different length scale parameters:

import numpy as np
import matplotlib.pyplot as plt
from gp4c import sample_prior, SamplingSpec

x = np.linspace(0, 5, 100)
length_scales = [0.2, 0.5, 1.0, 2.0]

fig, axes = plt.subplots(2, 2, figsize=(12, 10))
axes = axes.flatten()

for i, ell in enumerate(length_scales):
    spec = SamplingSpec(x_f=x)
    result = sample_prior(spec, ell=ell, n_samples=5)

    axes[i].plot(x, result.f.T, alpha=0.7)
    axes[i].set_title(f'Length scale ell = {ell}')
    axes[i].grid(True)

plt.tight_layout()
plt.show()

Run: 

python examples/parameter_study.py

Derivative Sampling

Sample f and its derivative h:

import numpy as np
import matplotlib.pyplot as plt
from gp4c import sample_prior, SamplingSpec

x = np.linspace(0, 5, 100)
spec = SamplingSpec(x_f=x, x_h=x)
result = sample_prior(spec, ell=0.5, n_samples=3)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))

ax1.plot(x, result.f.T, alpha=0.7)
ax1.set_title('f(x)')
ax1.grid(True)

ax2.plot(x, result.h.T, alpha=0.7)
ax2.set_title("h(x) = f'(x)")
ax2.grid(True)

plt.tight_layout()
plt.show()

Posterior Regression

Condition on sparse observations:

import numpy as np
import matplotlib.pyplot as plt
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)

# Test points
x_test = np.linspace(0, 5, 100)

# Posterior sampling
obs = Observations(x_f=x_train, y_f=y_train, noise_f=0.01)
spec = SamplingSpec(x_f=x_test)
result = sample_posterior(obs, spec, ell=1.0, n_samples=20)

# Plot
plt.figure(figsize=(10, 6))
plt.plot(x_test, result.f.T, color='lightblue', alpha=0.3)
plt.plot(x_test, result.f_mean, 'b-', linewidth=2, label='Posterior mean')
plt.plot(x_train, y_train, 'ro', markersize=8, label='Observations')
plt.legend()
plt.grid(True)
plt.show()

Hyperparameter Optimization

Use fit() to automatically find the kernel parameters that best explain the data, then feed the result into sample_posterior for a fully data-driven model.

import numpy as np
import gp4c

rng = np.random.default_rng(42)

# Noisy observations of a sine wave
x_train = np.linspace(0, 10, 30)
y_train = np.sin(x_train) + rng.normal(0, 0.1, len(x_train))

obs = gp4c.Observations(x_f=x_train, y_f=y_train, noise_f=0.01)

# Find optimal hyperparameters
result = gp4c.fit(obs, kernel='rbf')
print(f"sigma2 = {result.sigma2:.3f}")
print(f"ell    = {result.ell:.3f}")
print(f"log p(y|X,θ) = {result.log_likelihood:.3f}")

# Posterior at fine grid using optimized parameters
x_test = np.linspace(0, 12, 300)
spec = gp4c.SamplingSpec(x_f=x_test)

posterior = gp4c.sample_posterior(
    obs, spec,
    sigma2=result.sigma2,
    ell=result.ell,
    kernel='rbf',
    n_samples=50,
)

# Plot
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 5))
plt.fill_between(
    x_test,
    posterior.f_mean - 2 * posterior.f_std,
    posterior.f_mean + 2 * posterior.f_std,
    alpha=0.3, label='95% credible band',
)
plt.plot(x_test, posterior.f_mean, label='Posterior mean')
plt.scatter(x_train, y_train, color='red', zorder=5, label='Observations')
plt.legend()
plt.grid(True)
plt.title(f'RBF GP (sigma2={result.sigma2:.2f}, ell={result.ell:.2f})')
plt.show()

For periodic data, constrain the period search range:

# True period is 2.0
x_p = np.linspace(0, 20, 80)
y_p = np.sin(np.pi * x_p) + rng.normal(0, 0.1, len(x_p))  # period = 2.0

obs_p = gp4c.Observations(x_f=x_p, y_f=y_p, noise_f=0.01)

result_p = gp4c.fit(
    obs_p,
    kernel='periodic',
    bounds={'period': (1.0, 4.0)},
    method='differential_evolution',
    seed=42,
)
print(f"Recovered period = {result_p.period:.3f}")  # Should be close to 2.0

Interactive Notebooks

See notebooks/demo_derivative_sampling.ipynb for interactive demonstrations.

Run: 

jupyter notebook notebooks/demo_derivative_sampling.ipynb

Next Steps