Fit an approximate heteroscedastic Gaussian process

Fits two coupled Gaussian processes: one for the latent regression function and one for the logarithm of input-dependent observation-noise variance. Noise targets are constructed from exact leave-one-out residuals of the current mean GP, standardized by their leave-one-out variance so that the latent function's uncertainty is not counted as noise, and corrected for the known mean of a log chi-squared variable with one degree of freedom.

Usage

fit_heteroscedastic_gp(
  x,
  y,
  kernel,
  noise_kernel,
  mean = zero_mean(),
  initial_noise_variance = NULL,
  noise_floor = 1e-06,
  noise_ceiling = Inf,
  damping = 0.5,
  max_iterations = 50L,
  convergence_tolerance = 0.001,
  initial_jitter = 0,
  fallback_jitter = 1e-10,
  jitter_multiplier = 10,
  max_attempts = 8L,
  symmetry_tolerance = sqrt(.Machine$double.eps)
)

Arguments

x

Numeric vector or matrix of training inputs.

y

Numeric response vector.

kernel

Kernel specification for the latent regression function.

noise_kernel

Kernel specification for log observation-noise variance.

mean

Mean specification for the latent regression function, with fixed coefficients.

initial_noise_variance

Optional positive scalar used for the initial homoscedastic mean fit. If NULL, a small fraction of response variance is used, bounded below by noise_floor.

noise_floor

Strictly positive lower bound for estimated noise variances.

noise_ceiling

Optional upper bound for estimated noise variances.

damping

Value in \((0, 1]\) controlling how strongly each noise update replaces the previous estimate.

max_iterations

Maximum number of mean/noise refinement iterations. In a simulation study every fit that converged did so within 50 iterations, and most within 15.

convergence_tolerance

Convergence tolerance on the maximum absolute log change in training noise variance.

initial_jitter

Non-negative jitter tried first, relative to the covariance scale (the mean of its diagonal).

fallback_jitter

Positive relative jitter tried after the first failed factorization.

jitter_multiplier

Multiplicative jitter escalation factor.

max_attempts

Maximum Cholesky attempts.

symmetry_tolerance

Relative tolerance for checking that covariance matrices are symmetric.

Value

An object of class gaussianprocesses_heteroscedastic_model.

Details

This is a deterministic plug-in approximation inspired by most-likely heteroscedastic GP methods. It is not full joint Bayesian inference over the latent function and noise process.

For a Gaussian residual \(e \sim N(0, \sigma^2)\), \(\log(e^2)\) has additive log-chi-squared noise with known mean \(\psi(1/2) + \log 2\) and variance \(\psi_1(1/2)\). The implementation uses these known moments when constructing and fitting the log-noise GP.

A leave-one-out residual \(e_i\) has variance \(v_i + \sigma_i^2\), where \(v_i\) is the leave-one-out variance of the latent function, so \(\log(e_i^2)\) alone overestimates the noise. The target is therefore \(\log(e_i^2) - \psi(1/2) - \log 2 - \log((v_i + s_i^2) / s_i^2)\), with \(s_i^2\) the current noise estimate, which is unbiased at a fixed point of the iteration. The research note system.file("notes", "heteroscedastic-gp.md", package = "gaussianprocesses") derives the correction and reports a simulation study of its bias and convergence.

The approach alternates plug-in estimates rather than integrating jointly over the latent function and noise process. If the estimate has not converged after max_iterations iterations, the function warns with class gaussianprocesses_convergence_warning and returns the last iterate with converged = FALSE.

Stability

Experimental: this interface may change in a minor release, and every change is listed in NEWS. See gaussianprocesses-package for the policy.

See also

gaussianprocesses_model for the structure, versioning, and persistence of fitted models.

Examples

simulation <- gp_simulation_scenario("heteroscedastic_1d", n = 60)

model <- fit_heteroscedastic_gp(
  simulation$x,
  simulation$observed,
  kernel = matern52_kernel(length_scale = 0.2),
  noise_kernel = rbf_kernel(length_scale = 0.4),
  max_iterations = 30
)
#> Warning: The heteroscedastic noise estimate did not converge in 30 iterations: the last largest change in log noise variance was 0.0326, above convergence_tolerance = 0.001. Increase max_iterations.
model
#> Approximate heteroscedastic Gaussian-process model
#>   observations: 60
#>   iterations: 30
#>   converged: no
#>   training noise variance range: [0.0168854, 0.134745]
#>   approximation: iterative_log_residual_plugin

# The estimated noise variance rises with x, as in the simulation.
rows <- c(1, 20, 40, 60)
cbind(
  x = simulation$x[rows],
  true = simulation$noise_variance[rows],
  estimated = model$training_noise_variance[rows]
)
#>              x       true  estimated
#> [1,] 0.0000000 0.01500000 0.01688545
#> [2,] 0.3220339 0.04092646 0.03880570
#> [3,] 0.6610169 0.12423585 0.10815773
#> [4,] 1.0000000 0.26500000 0.13068690