API Reference: Functions¶

Prior Sampling¶

sample_prior¶

sample_prior(spec, sigma2=1.0, ell=1.0, kernel='rbf', n_samples=1, seed=42, period=None, ell_p=None)

Sample from the prior distribution of f, g=integral(f), h=f', and/or u=f''.

Parameters:

Parameter	Type	Default	Description
`spec`	`SamplingSpec`	—	Which functions to sample and at which points
`sigma2`	`float`	`1.0`	Kernel variance
`ell`	`float`	`1.0`	Length scale (SE decay scale for locally-periodic)
`kernel`	`str`	`'rbf'`	Kernel type (see below)
`n_samples`	`int`	`1`	Number of joint samples
`seed`	`int`	`42`	Random seed
`period`	`float \\| None`	`None`	Period for periodic/locally-periodic kernels
`ell_p`	`float \\| None`	`None`	Periodic length scale; defaults to `ell` if not set

Supported kernels:

`kernel`	Description	Supports h?	Supports g?	Supports u?
`'rbf'`	RBF / Squared Exponential	Yes	Yes	Yes
`'matern52'`	Matérn 5/2 (twice differentiable)	Yes	Yes	Yes
`'matern32'`	Matérn 3/2 (once mean-square differentiable)	Yes	Yes	No
`'periodic'`	Periodic (requires `period`)	Yes	No	Yes
`'locally_periodic'`	Locally-periodic (requires `period`)	Yes	No	Yes

Returns: GPSamples — named tuple with fields f, g, h, u (each ndarray or None)

f: shape (n_samples, n_f) if spec.x_f is set, else None
g: shape (n_samples, n_g) if spec.x_g is set, else None
h: shape (n_samples, n_h) if spec.x_h is set, else None
u: shape (n_samples, n_u) if spec.x_u is set, else None

Example:

RBFPeriodicLocally-Periodic

import numpy as np
from gp4c import sample_prior, SamplingSpec

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

print(result.f.shape)  # (10, 100)
print(result.g.shape)  # (10, 100)
print(result.h.shape)  # (10, 100)

import numpy as np
from gp4c import sample_prior, SamplingSpec

x = np.linspace(0, 10, 200)
spec = SamplingSpec(x_f=x)
result = sample_prior(spec, kernel='periodic', period=2.0, sigma2=1.0, ell=1.0, n_samples=5)

print(result.f.shape)  # (5, 200)

import numpy as np
from gp4c import sample_prior, SamplingSpec

x = np.linspace(0, 10, 200)
spec = SamplingSpec(x_f=x)
# ell controls the SE envelope decay; ell_p controls the periodic wiggle scale
result = sample_prior(spec, kernel='locally_periodic', period=2.0, ell=5.0, ell_p=1.0, n_samples=5)

Posterior Sampling¶

sample_posterior¶

sample_posterior(obs, spec, sigma2=1.0, ell=1.0, kernel='rbf', n_samples=1, seed=42, period=None, ell_p=None)

Sample from the posterior distribution conditioned on observations.

Parameters:

Parameter	Type	Default	Description
`obs`	`Observations`	—	Observed data; at least one of (x_f,y_f), (x_g,y_g), (x_h,y_h), (x_u,y_u) required
`spec`	`SamplingSpec`	—	Prediction points for f, g, h, and/or u
`sigma2`	`float`	`1.0`	Kernel variance
`ell`	`float`	`1.0`	Length scale
`kernel`	`str`	`'rbf'`	Kernel type (same options as `sample_prior`)
`n_samples`	`int`	`1`	Number of posterior samples
`seed`	`int`	`42`	Random seed
`period`	`float \\| None`	`None`	Period for periodic/locally-periodic kernels
`ell_p`	`float \\| None`	`None`	Periodic length scale; defaults to `ell` if not set

Returns: PosteriorSamples — named tuple with fields:

Field	Shape	Description
`f`, `g`, `h`, `u`	`(n_samples, n_x)` or `None`	Posterior samples
`f_mean`, `g_mean`, `h_mean`, `u_mean`	`(n_x,)` or `None`	Analytical posterior means
`f_std`, `g_std`, `h_std`, `u_std`	`(n_x,)` or `None`	Analytical posterior standard deviations

Example:

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])
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, x_h=x_test)

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

print(result.f.shape)       # (100, 100) — samples
print(result.f_mean.shape)  # (100,) — posterior mean
print(result.f_std.shape)   # (100,) — posterior std dev
print(result.h.shape)       # (100, 100) — derivative samples

Hyperparameter Optimization¶

fit¶

fit(obs, kernel='rbf', bounds=None, method='differential_evolution', x0=None, **optimizer_kwargs)

Optimize GP hyperparameters by maximizing the log marginal likelihood.

Finds the hyperparameters (sigma2, ell, and optionally period) that best explain the observed data. Supports both global optimization strategies (recommended) and local methods (faster when a good initial guess is available).

Parameters:

Parameter	Type	Default	Description
`obs`	`Observations`	—	Training data; f and/or h observations. g observations are not supported.
`kernel`	`str`	`'rbf'`	Kernel type (same options as `sample_prior`)
`bounds`	`dict \\| None`	`None`	Parameter bounds as `{'sigma2': (lo, hi), 'ell': (lo, hi), 'period': (lo, hi)}`. Merged with defaults if partial.
`method`	`str`	`'differential_evolution'`	Optimization method (see below)
`x0`	`dict \\| None`	`None`	Initial parameter values; required for local methods
`**optimizer_kwargs`	—	—	Forwarded to the underlying scipy optimizer

Default bounds:

Parameter	Default range
`sigma2`	`(1e-4, 100.0)`
`ell`	`(1e-4, 10.0)`
`period`	`(1e-2, 100.0)` (periodic kernels only)

Supported optimization methods:

Global (recommended)Local

Method	Notes
`'differential_evolution'`	Genetic algorithm; robust, no initial guess needed
`'dual_annealing'`	Simulated annealing variant
`'shgo'`	Simplicial homology global optimization
`'basinhopping'`	Basin-hopping; uses `x0` if provided, else midpoint of bounds

Method	Notes
`'L-BFGS-B'`	Quasi-Newton; fast, requires `x0`
`'Nelder-Mead'`	Derivative-free simplex; requires `x0`
`'Powell'`	Derivative-free; requires `x0`
`'SLSQP'`	Sequential Least Squares; requires `x0`
Any `scipy.optimize.minimize` method	Requires `x0`

Returns: OptimizationResult — dataclass with:

Field	Type	Description
`sigma2`	`float`	Optimized kernel variance
`ell`	`float`	Optimized length scale
`period`	`float \\| None`	Optimized period (None for non-periodic kernels)
`log_likelihood`	`float`	Log marginal likelihood at optimum
`success`	`bool`	Whether the optimizer converged
`scipy_result`	`Any`	Raw scipy optimization result

Warning

Integral observations (obs.x_g) are not supported by fit(). Use only f and/or h observations.

Note

The likelihood surface can be non-convex. Global methods are recommended for initial hyperparameter search. Once a good region is found, a local method can refine the result.

Examples:

Basic RBFL-BFGS-BPeriodic kernelDerivative observations

import numpy as np
import gp4c

x = np.linspace(0, 10, 50)
y = np.sin(x) + np.random.default_rng(0).normal(0, 0.1, len(x))
obs = gp4c.Observations(x_f=x, y_f=y, noise_f=0.01)

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

result = gp4c.fit(
    obs,
    kernel='rbf',
    method='L-BFGS-B',
    x0={'sigma2': 1.0, 'ell': 1.0},
    bounds={'sigma2': (0.1, 5.0), 'ell': (0.1, 3.0)},
)

result = gp4c.fit(
    obs,
    kernel='periodic',
    bounds={'period': (0.5, 4.0)},
    method='differential_evolution',
    maxiter=300,
    seed=42,
)
print(f"period={result.period:.3f}")

x_train = np.linspace(0, 5, 20)
obs = gp4c.Observations(
    x_f=x_train, y_f=np.sin(x_train),   noise_f=0.01,
    x_h=x_train, y_h=np.cos(x_train),   noise_h=0.01,
)
result = gp4c.fit(obs, kernel='matern52')

Log Marginal Likelihood¶

log_marginal_likelihood / LML¶

log_marginal_likelihood(obs, sigma2, ell, ell_p, period, kernel) -> float

Compute the log marginal likelihood: \(\log p(\mathbf{y} \mid X, \theta)\).

\[\log p(\mathbf{y} \mid X, \theta) = -\tfrac{1}{2}\mathbf{y}^\top K^{-1}\mathbf{y} - \tfrac{1}{2}\log|K| - \tfrac{n}{2}\log 2\pi\]

This is a thin Python wrapper around the C implementation. LML is an alias for log_marginal_likelihood.

Parameters:

Parameter	Type	Description
`obs`	`Observations`	Training data
`sigma2`	`float`	Kernel variance
`ell`	`float`	SE decay length scale
`ell_p`	`float \\| None`	Periodic length scale; pass `None` to use `ell`
`period`	`float`	Period (for periodic kernels; pass any value for non-periodic)
`kernel`	`str`	Kernel type

Returns: float — log marginal likelihood (more positive is better).

Example:

import numpy as np
import gp4c

obs = gp4c.Observations(x_f=np.linspace(0, 5, 20), y_f=np.sin(np.linspace(0, 5, 20)), noise_f=0.01)

lml = gp4c.log_marginal_likelihood(obs, sigma2=1.0, ell=0.5, ell_p=None, period=1.0, kernel='rbf')
print(f"log p(y|X,θ) = {lml:.3f}")

# Alias
lml2 = gp4c.LML(obs, sigma2=1.0, ell=0.5, ell_p=None, period=1.0, kernel='rbf')

Kernel Parameters¶

sigma2 (variance)¶

Controls the amplitude of GP samples.

Small values (e.g., 0.1): Low-amplitude fluctuations
Large values (e.g., 10.0): High-amplitude fluctuations
Default: 1.0

ell (length scale)¶

Controls the smoothness and correlation distance.

Small values (e.g., 0.1): Rapidly varying, rough samples
Large values (e.g., 5.0): Slowly varying, smooth samples
Default: 1.0

Warning

Very small length scales with derivatives can cause numerical instability. Use ell >= 0.2 for derivative sampling.

period (periodic kernels only)¶

Controls the repetition distance for 'periodic' and 'locally_periodic' kernels.

Required when using either periodic kernel; no default is provided.
For 'locally_periodic', the ell parameter additionally controls an SE envelope that damps the periodicity over long distances.

Next Steps¶

Types Reference — Data structures and type definitions
Hyperparameter Optimization — Full optimization API
Tutorials — Usage examples