Package 'gaussianprocesses'

Description

Measures exact GP fit time and retained model size over increasing training sizes using seeded RBF simulations. This function is intended for local empirical scaling studies rather than CI performance gating.

Usage

benchmark_exact_gp_scaling(
  sizes = c(50L, 100L, 200L, 400L),
  repeats = 3L,
  seed = 1729L
)

Arguments

Value

A data frame with one row per requested size.

Stability

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

Examples

# A deliberately small run; larger sizes show the cubic cost.
benchmark_exact_gp_scaling(sizes = c(20, 40), repeats = 1)

Description

Runs both models with identical fixed kernel and noise parameters, then compares retained model size, elapsed fitting/prediction time, posterior mean RMSE, and latent-variance MAE. The exact GP is treated as the numerical reference, not as ground truth.

Usage

benchmark_sparse_gp(
  x,
  y,
  kernel,
  n_inducing,
  prediction_x = x,
  noise_variance = 1e-06,
  mean = zero_mean(),
  selection = c("farthest", "quantile", "variance"),
  observation_noise_variance = NULL
)

Arguments

Details

The benchmark reports retained R object size, not peak process memory. Elapsed timings are simple wall-clock measurements intended for local comparison, not publication-grade benchmarking.

Value

A list containing a comparison data frame and both prediction objects.

Stability

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

Examples

x <- seq(-3, 3, length.out = 200)
y <- sin(x) + 0.1 * cos(7 * x)

benchmark <- benchmark_sparse_gp(
  x,
  y,
  kernel = rbf_kernel(length_scale = 0.8),
  n_inducing = 15,
  noise_variance = 0.01
)
benchmark$comparison

Description

Switches from one covariance to another along one input column, smoothly around an estimable location: a structural break in a time series, or a regime change in space.

Usage

changepoint_kernel(before, after, location = 0, steepness = 1, column = 1L)

Arguments

Details

With the sigmoid $s(x) = 1 / (1 + e^{-a (x_t - c)})$ along column $t$,

$k(x, x') = (1 - s(x)) k_1(x, x') (1 - s(x')) + s(x) k_2(x, x') s(x'),$

where $k_1$ is before and $k_2$ is after. Each term has the form $g(x) k(x, x') g(x')$, so the kernel is positive semidefinite. The two regimes are independent processes: as $a \to \infty$, the covariance between points on opposite sides of $c$ tends to 0.

location is a real-valued parameter, optimized on the identity scale, and steepness a positive one, optimized on the log scale. Parameter paths are location and steepness for the change itself, and before. and after. for the two kernels, so nested changepoints model several changes. optimize_gp() bounds every changepoint location by the training range of its column and, because the likelihood is often multimodal in the location, starts it at evenly spaced quantiles of that column.

The kernel is as smooth as its two components, so it supports kernel_input_gradient(), predict_gradient_gp(), and derivative observations when they do.

Value

A composite kernel specification.

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

time_series_kernel(), whose changepoint argument adds a change to the time-series kernel.

Examples

kernel <- changepoint_kernel(
  before = rbf_kernel(variance = 1, length_scale = 2),
  after = rbf_kernel(variance = 0.2, length_scale = 0.5),
  location = 5,
  steepness = 4
)
kernel
names(kernel_parameters(kernel, flatten = TRUE))

# Points on opposite sides of a steep change are nearly uncorrelated.
evaluate_kernel(kernel, c(4, 6))

Description

coef() returns the coefficients $\beta$ of a model's mean function, and vcov() their covariance given the hyperparameters.

Usage

## S3 method for class 'gaussianprocesses_latent_model'
coef(object, ...)

## S3 method for class 'gaussianprocesses_model'
coef(object, ...)

## S3 method for class 'gaussianprocesses_model'
vcov(object, ...)

Arguments

Details

For estimated coefficients, coef() is the generalized least-squares estimate $\hat\beta$ and vcov() its covariance $(H^\top C^{-1} H)^{-1}$. For marginalized coefficients they are the posterior mean $\bar\beta$ and covariance $(B^{-1} + H^\top C^{-1} H)^{-1}$. Fixed coefficients have zero covariance. Both treat the kernel and noise parameters as known. The zero mean has no coefficients.

For latent models (fit_latent_gp()), coef() returns the fixed coefficients, the estimate that maximizes the approximate marginal likelihood, or, under a Gaussian prior $N(b, B)$, the approximate posterior mean $b + B H^\top a$ with $a = \nabla \log p(y \mid \hat f)$. They have no vcov() method.

Value

coef(): a named numeric vector. vcov(): a matrix with matching row and column names.

Stability

Stable: from version 1.0.0 this interface changes incompatibly only in a major release, after a deprecation period. Results and options that concern an experimental model class, kernel, or argument follow that interface's tier. See gaussianprocesses-package for the policy.

See Also

mean_coefficients, linear_mean().

Examples

x <- seq(0, 4, length.out = 25)
model <- fit_gp(x, 1 + 0.5 * x + sin(2 * x), rbf_kernel(length_scale = 0.5),
  noise_variance = 0.01, mean = linear_mean())

coef(model)
vcov(model)

Description

Compare two GP posterior predictions numerically

Usage

compare_gp_predictions(reference, candidate)

Arguments

Value

A named list of numerical accuracy differences.

Stability

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

Examples

x <- seq(-3, 3, length.out = 100)
y <- sin(x)
kernel <- rbf_kernel(length_scale = 0.8)

exact <- predict_gp(fit_gp(x, y, kernel, noise_variance = 0.01), x)
sparse <- predict_sparse_gp(
  fit_sparse_gp(x, y, kernel, noise_variance = 0.01, n_inducing = 10),
  x
)

compare_gp_predictions(exact, sparse)

Description

A coregionalization kernel acts on one input column that holds an output index $p \in \{1, \ldots, P\}$ and gives the covariance between outputs, $k(p, q) = B_{pq}$ with

$B = W W^\top + \operatorname{diag}(\kappa),$

which is positive semidefinite by construction. Multiplied by a kernel on the other inputs it gives the intrinsic coregionalization model (ICM), $\operatorname{cov}(f_p(x), f_q(x')) = B_{pq}\, k(x, x')$; a sum of such products is the linear model of coregionalization (LMC).

Usage

coregionalization_kernel(n_outputs, rank = 1L, W = NULL, kappa = NULL)

Arguments

Details

The kernel needs exactly one input column, so in a model with other inputs it is restricted to the output-index column with select_dimensions(). With inputs in column 1 and the output index in column 2, an ICM kernel is product_kernel(select_dimensions(rbf_kernel(), 1), select_dimensions(coregionalization_kernel(P), 2)). gp_stack_outputs() builds such inputs from one response column per output, and initialize_coregionalization() gives a starting $W$ and $\kappa$.

The optimizer works on $W$ directly and on $\log \kappa$. Two cautions apply:

• $W$ is identified only up to rotation and sign: $W Q$ gives the same $B$ for any orthogonal $Q$. Report $B$, with coregionalization_matrix(), not $W$.

• $W = 0$ is a stationary point of the likelihood, where the gradient with respect to $W$ is zero, so an optimization started there keeps $W = 0$; start from initialize_coregionalization() instead.

The overall scale of an ICM term is shared by $B$ and the variance of the input kernel; fixing the latter at 1 removes the redundancy only if that parameter is not optimized, so expect the two to trade off.

Value

A gaussianprocesses_kernel.

Stability

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

References

Bonilla, E. V., Chai, K. M. A., and Williams, C. K. I. (2008). Multi-task Gaussian process prediction. Advances in Neural Information Processing Systems, 20.

Alvarez, M. A., Rosasco, L., and Lawrence, N. D. (2012). Kernels for vector-valued functions: a review. Foundations and Trends in Machine Learning, 4(3), 195–266.

See Also

gp_stack_outputs(), initialize_coregionalization(), coregionalization_matrix().

Examples

kernel <- coregionalization_kernel(
  n_outputs = 2,
  rank = 1,
  W = c(1, 0.5),
  kappa = c(0.1, 0.2)
)
coregionalization_matrix(kernel)
evaluate_kernel(kernel, matrix(c(1, 2, 2)))

# An ICM kernel for inputs in column 1 and the output index in column 2.
icm <- product_kernel(
  select_dimensions(rbf_kernel(length_scale = 0.5), 1),
  select_dimensions(kernel, 2)
)

Description

Returns $B = W W^\top + \operatorname{diag}(\kappa)$ for every coregionalization kernel in a kernel, such as the fitted kernel of a multi-output model.

Usage

coregionalization_matrix(kernel)

Arguments

Details

$B$ is identified, while $W$ is not: $W Q$ gives the same $B$ for any orthogonal $Q$. In an ICM term $B \otimes k$, the scale of $B$ also trades off with the variance of the input kernel $k$; the correlations stats::cov2cor(B) do not.

Value

A named list with one $P \times P$ matrix per coregionalization kernel, named by its parameter path ("" for a coregionalization kernel itself).

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

coregionalization_kernel().

Examples

kernel <- sum_kernel(
  product_kernel(
    select_dimensions(rbf_kernel(), 1),
    select_dimensions(coregionalization_kernel(2, W = c(1, 1)), 2)
  ),
  select_dimensions(coregionalization_kernel(2, kappa = c(0.5, 2)), 2)
)
coregionalization_matrix(kernel)

Description

The package offers two ways to compute a covariance matrix. Kernel specifications, such as rbf_kernel(), are objects: they can be combined with sum_kernel(), product_kernel(), and scale_kernel(), their parameters can be read, updated, and optimized, and models store them. evaluate_kernel() computes their covariance. The covariance functions, such as kernel_rbf(), take the parameter values as arguments and return one matrix directly. Both use the same formulas.

Usage

evaluate_kernel(kernel, x, y = NULL)

Arguments

Value

A numeric covariance matrix.

Stability

Stable: from version 1.0.0 this interface changes incompatibly only in a major release, after a deprecation period. Results and options that concern an experimental model class, kernel, or argument follow that interface's tier. See gaussianprocesses-package for the policy.

Examples

kernel <- sum_kernel(rbf_kernel(length_scale = 0.5), white_noise_kernel(0.1))
x <- c(0, 0.5, 1)

evaluate_kernel(kernel, x)

# Cross-covariance between training and new inputs.
evaluate_kernel(kernel, x, c(0.25, 2))

Description

Computes $\log p(y_i \mid f_i)$ and its first three derivatives with respect to the latent values $f_i$, the quantities that Laplace and related approximations need.

Usage

evaluate_likelihood(likelihood, y, f, exposure = NULL)

Arguments

Value

A data frame with one row per response: y in the internal coding (0/1 for binary responses), f, log_density, and the derivatives first, second, and third.

Stability

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

Examples

f <- c(-40, -2, 0, 2, 40)
evaluate_likelihood(bernoulli_likelihood("probit"), rep(1, 5), f)
evaluate_likelihood(poisson_likelihood(), c(0, 1, 2, 3, 4), f = log(1:5),
  exposure = 2)

Description

Evaluate a Gaussian-process mean function

Usage

evaluate_mean(mean, x)

Arguments

Details

Means with fixed coefficients evaluate to $h(x)^\top \beta$. A mean whose coefficients have a Gaussian prior evaluates to its prior mean $h(x)^\top b$. Estimated coefficients, and those with a vague prior, have no value before a model is fitted; predict_gp() and coef() use the fitted values.

Value

A numeric vector with one mean value per observation.

Stability

Stable: from version 1.0.0 this interface changes incompatibly only in a major release, after a deprecation period. Results and options that concern an experimental model class, kernel, or argument follow that interface's tier. See gaussianprocesses-package for the policy.

Examples

evaluate_mean(constant_mean(1), matrix(0, nrow = 4, ncol = 2))
evaluate_mean(linear_mean(c(1, 2)), c(0, 0.5, 1))

Description

Fits a Gaussian-process regression model by conditioning a Gaussian prior on observed responses. The implementation uses a Cholesky factorization and triangular solves; it never forms an explicit matrix inverse.

Usage

fit_gp(
  x,
  y,
  kernel,
  noise_variance = 1e-06,
  mean = zero_mean(),
  initial_jitter = 0,
  fallback_jitter = 1e-10,
  jitter_multiplier = 10,
  max_attempts = 8L,
  symmetry_tolerance = sqrt(.Machine$double.eps),
  derivative = NULL
)

Arguments

Details

A scalar noise_variance gives the usual homoscedastic model. A vector adds a known diagonal observation-noise variance for each training point. Observation noise and numerical jitter are always stored separately.

Value

An object of class gaussianprocesses_model.

Derivative observations

Differentiation is linear, so derivatives of a Gaussian process are jointly Gaussian with its values, and observing them is exact conditioning. With values at $X$ and derivatives at $X'$,

$\begin{pmatrix} y \\ y' \end{pmatrix} \sim N\left( \begin{pmatrix} m(X) \\ \partial_d m(X') \end{pmatrix}, \begin{pmatrix} K(X, X) & \partial_{x'_d} K(X, X') \\ \partial_{x_d} K(X', X) & \partial_{x_d} \partial_{x'_e} K(X', X') \end{pmatrix} + \Sigma \right),$

where each derivative row uses its own direction $d$. The blocks come from the kernel's input derivatives (kernel_input_gradient()), so the kernel must be mean-square differentiable; otherwise an error of class gaussianprocesses_smoothness_error is raised. Parametric means contribute their input derivatives. noise_variance may differ between value and derivative observations.

Conditioning, prediction, log_marginal_likelihood(), its gradient, and optimize_gp() work unchanged with this covariance; predict_gp() and predict_gradient_gp() predict values and derivatives at new inputs. Derivative blocks scale with $1 / \ell^2$, so gp_numerical_diagnostics() reports the conditioning of each observation type separately.

Numerical jitter

If $C = K + \Sigma_\varepsilon$ does not factorize, the factorization is retried with $C + \delta I$, where $\delta = \epsilon s$ and $s$ is the covariance scale, the mean of the diagonal of $C$. The relative jitter $\epsilon$ starts at initial_jitter (if positive), then takes the values fallback_jitter, fallback_jitter * jitter_multiplier, and so on, for at most max_attempts attempts. Because $\delta$ is proportional to $C$, results do not depend on the units of y: multiplying y by $c$ and every variance parameter by $c^2$ multiplies the jitter by $c^2$. The fitted model reports the absolute jitter $\delta$ in jitter, and gp_numerical_diagnostics() also reports $\epsilon$ and $s$.

Stability

Stable: from version 1.0.0 this interface changes incompatibly only in a major release, after a deprecation period. Results and options that concern an experimental model class, kernel, or argument follow that interface's tier. See gaussianprocesses-package for the policy. The derivative argument (derivative observations) is experimental.

See Also

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

Examples

x <- seq(-2, 2, length.out = 15)
y <- sin(2 * x)

model <- fit_gp(x, y, kernel = rbf_kernel(length_scale = 0.6), noise_variance = 0.01)
model

# Known, observation-specific noise variances.
fit_gp(
  x,
  y,
  kernel = rbf_kernel(length_scale = 0.6),
  noise_variance = seq(0.01, 0.1, length.out = 15)
)

Description

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

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.

Value

An object of class gaussianprocesses_heteroscedastic_model.

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
)
model

# 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]
)

Description

Fits a Gaussian process $f \sim GP(m, k)$ observed through a likelihood $p(y_i \mid f(x_i))$, such as Bernoulli responses for classification or Poisson counts, with the Laplace approximation to the posterior of $f$. With gaussian_likelihood() the approximation is exact and the model equals fit_gp().

Usage

fit_latent_gp(
  x,
  y,
  kernel,
  likelihood,
  mean = zero_mean(),
  exposure = NULL,
  method = "laplace",
  max_iterations = 100L,
  convergence_tolerance = 1e-14,
  initial_jitter = 0,
  fallback_jitter = 1e-10,
  jitter_multiplier = 10,
  max_attempts = 8L,
  symmetry_tolerance = sqrt(.Machine$double.eps)
)

Arguments

Details

The posterior of $f = f(X)$ is proportional to $p(y \mid f) N(f \mid m, K)$. The Laplace approximation replaces it by $N(\hat f, (K^{-1} + W)^{-1})$, where $\hat f$ maximizes

$\Psi(f) = \log p(y \mid f) - \tfrac12 (f - m)^\top K^{-1} (f - m)$

and $W = -\nabla\nabla \log p(y \mid \hat f)$ is diagonal. For a log-concave likelihood $\Psi$ is strictly concave, so the mode is unique; other likelihoods raise an error of class gaussianprocesses_log_concavity_error.

Newton's method finds the mode (Rasmussen and Williams, 2006, Algorithm 3.1) in the parameterization $f = m + K a$. Each iteration factorizes

$B = I + W^{1/2} K W^{1/2},$

whose eigenvalues are at least 1, so it is well conditioned even when $K$ is nearly singular; the factorization follows the jitter policy of fit_gp(). A step-halving line search on $\Psi$ prevents overshooting. The iterations stop when the Newton decrement

$\lambda^2 = g^\top (K^{-1} + W)^{-1} g, \qquad g = \nabla \Psi(f),$

which is twice the increase of $\Psi$ that the next Newton step predicts, is at most convergence_tolerance. $\lambda$ is about the distance to the mode measured in posterior standard deviations, so the default of $10^{-14}$ places the mode within about $10^{-7}$ standard deviations. ⁠$laplace⁠ records the number of iterations and step halvings, the final decrement, the largest component of $g$, the last change of $\Psi$, and whether the iterations converged; if they did not, a warning of class gaussianprocesses_convergence_warning is raised.

The approximate log marginal likelihood is

$\log q(y \mid X, \theta) = \log p(y \mid \hat f) - \tfrac12 (\hat f - m)^\top K^{-1} (\hat f - m) - \tfrac12 \log |B|,$

reported by log_marginal_likelihood() and logLik(), with its gradient in log_marginal_likelihood_gradient().

A mean with fixed coefficients contributes $m(x)$. Coefficients with a Gaussian prior $N(b, B_\beta)$ are integrated out: the prior becomes $m(x) = h(x)^\top b$ with covariance $k(x, x') + h(x)^\top B_\beta h(x')$. Coefficients estimated by maximum likelihood (estimate_coefficients("ml"), the default of linear_mean()) maximize $\log q(y \mid X, \theta, \beta)$ at the given hyperparameters, by BFGS with the analytical gradient $H^\top a + [(I + K W)^{-1} H]^\top s$, where $s$ is the derivative of $-\frac12 \log |B|$ in $\hat f$, starting from the generalized linear model without the latent process. The estimate is plugged in, like the hyperparameters; coef() returns it and ⁠$mean_fit⁠ records it. A constant mean estimated this way is the usual mean level of a count or classification model. Restricted likelihoods and vague priors are not supported.

Value

An object of class gaussianprocesses_latent_model with fields x, y (responses in the internal coding), response_levels (the levels of a factor response), exposure, kernel, mean, likelihood, method, mode ($\hat f$), alpha ($a$), curvature (the diagonal of $W$), b_cholesky ($R$ with $R^\top R = B + \delta I$), jitter ($\delta$), log_marginal_likelihood, laplace, n_observations, n_features, schema_version, and package_version. Other fields are internal.

Accuracy of the Laplace approximation

For a few observations the exact marginal likelihood $\log p(y \mid X)$ can be computed by tensor-product Gauss-Hermite quadrature. With an RBF kernel (variance 2, length scale 1) and the first $n$ of the inputs $x = 0, 0.5, 2$ and of the responses below, the error $\log q - \log p$ of the Laplace approximation is:

For binary responses and small counts the error is a few hundredths and grows with $n$; the approximation can over- or underestimate. For counts in the hundreds, where the likelihood is nearly Gaussian in $f$, it is ten times smaller. The test suite reproduces these values.

Stability

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

References

Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press. Chapters 3 and 5.

See Also

predict_latent_gp(), optimize_latent_gp(), sample_gp_posterior().

Examples

x <- seq(-3, 3, length.out = 30)
y <- as.numeric(sin(x) + 0.3 * cos(5 * x) > 0)

model <- fit_latent_gp(x, y, rbf_kernel(variance = 4), bernoulli_likelihood())
model
predict_latent_gp(model, c(-2, 0, 2))$probability

Description

Fits a sparse approximation with fixed inducing locations: the Fully Independent Training Conditional (FITC) approximation or the collapsed variational free energy (VFE) of Titsias (2009). Exact GP regression remains the reference model; these provide a lower-cost path when the number of inducing points is much smaller than the number of observations.

Usage

fit_sparse_gp(
  x,
  y,
  kernel,
  noise_variance = 1e-06,
  mean = zero_mean(),
  n_inducing = NULL,
  inducing_points = NULL,
  selection = c("farthest", "quantile", "variance"),
  method = c("fitc", "vfe"),
  fitc_diagonal_floor = 1e-10,
  initial_jitter = 0,
  fallback_jitter = 1e-10,
  jitter_multiplier = 10,
  max_attempts = 8L,
  symmetry_tolerance = sqrt(.Machine$double.eps)
)

Arguments

Details

Let $m$ be the number of inducing points and $n$ the number of training observations. FITC replaces the full training covariance with a low-rank inducing covariance plus a diagonal correction. The dominant training cost is $O(n m^2 + m^3)$ rather than the exact GP's $O(n^3)$, with $m \ll n$. The implementation stores an $m \times n$ cross covariance rather than an $n \times n$ training covariance.

FITC preserves each training point's prior marginal variance through the diagonal correction, unlike simpler deterministic training conditional approximations. It is still an approximation: inducing locations and their number can affect posterior means, variances, and the approximate marginal likelihood.

If the inducing covariance needs numerical jitter, the jittered matrix is used consistently in the low-rank term, the posterior, and the approximate log marginal likelihood. The inducing posterior is factorized in whitened form, whose eigenvalues are bounded below by one.

Value

An object of class gaussianprocesses_sparse_model.

VFE

With $Q_{ff} = K_{fu} K_{uu}^{-1} K_{uf}$ and noise $\Sigma = \mathrm{diag}(\sigma_i^2)$, the collapsed variational bound of Titsias (2009) is

$\mathcal F = \log N(y \mid m(X), Q_{ff} + \Sigma) - \frac12 \sum_i \frac{k(x_i, x_i) - [Q_{ff}]_{ii}}{\sigma_i^2} \le \log p(y).$

It is a lower bound on the exact log marginal likelihood, reported as log_marginal_likelihood, so maximizing it over the hyperparameters cannot overfit the way FITC can: FITC's diagonal correction $\mathrm{diag}(K_{ff} - Q_{ff})$ can absorb noise, and FITC tends to underestimate the noise variance (Bauer, van der Wilk, and Rasmussen, 2016). VFE predicts with the deterministic training conditional, which has the same form as the FITC prediction with $\Lambda = \Sigma$.

Both models report the trace term $t = \sum_i (k(x_i, x_i) - [Q_{ff}]_{ii}) / \sigma_i^2$ in trace_term: the prior variance that the inducing points do not explain, in units of noise variance. A trace term much larger than 1 means that the inducing points do not cover the data well: VFE's bound is loose by at least $t/2$ and FITC leans on its diagonal correction. More or better placed inducing points reduce it.

Stability

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

References

Titsias, M. (2009). Variational learning of inducing variables in sparse Gaussian processes. Proceedings of the 12th International Conference on Artificial Intelligence and Statistics, 567–574.

Bauer, M., van der Wilk, M., and Rasmussen, C. E. (2016). Understanding probabilistic sparse Gaussian process approximations. Advances in Neural Information Processing Systems, 29.

See Also

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

Examples

x <- seq(-3, 3, length.out = 200)
y <- sin(x) + 0.1 * cos(7 * x)

model <- fit_sparse_gp(
  x,
  y,
  kernel = rbf_kernel(length_scale = 0.8),
  noise_variance = 0.01,
  n_inducing = 15
)
model

Description

Fit a GP to ordered time-series observations

Usage

fit_time_series_gp(
  time,
  y,
  kernel,
  unit = "auto",
  method = c("exact", "state_space"),
  ...
)

Arguments

Details

This function was called fit_gp_time() in version 0.1.0.

Value

A fitted model containing explicit time-index metadata used by forecast_gp(): a gaussianprocesses_model for method = "exact" and a gaussianprocesses_state_space_model for method = "state_space".

State-space inference

On one input, a Matérn process with $\nu = p + 1/2$ is the first component of the state $x = (f, f', \ldots, f^{(p)})$ of the linear stochastic differential equation $dx = F x\,dt + L\,dW$, with $F$ the companion matrix of $(s + \lambda)^{p+1}$, $\lambda = \sqrt{2\nu}/\ell$, and white noise $W$ of spectral density $q_c$ (Hartikainen and Särkkä, 2010). The stationary state covariance $P_\infty$ solves $F P_\infty + P_\infty F^\top + L q_c L^\top = 0$, and $k(\tau) = H e^{F\tau} P_\infty H^\top$ for $\tau \ge 0$ with $H = (1, 0, \ldots, 0)$. Sums of kernels stack their states, and scaled kernels scale $P_\infty$ and $q_c$. Between consecutive times the state moves by $A = e^{F\Delta t}$, which has a closed form because $F + \lambda I$ is nilpotent, with process noise $P_\infty - A P_\infty A^\top$.

The Kalman filter (with Joseph-form updates) gives the log marginal likelihood by the prediction-error decomposition, the Rauch-Tung-Striebel smoother gives the posterior at the training times, and forecast_gp() interleaves the requested times as unobserved steps. The results equal those of the exact GP up to rounding, at $O(n)$ cost instead of $O(n^3)$, and no numerical jitter is used. sample_gp_posterior() draws by forward filtering, backward sampling. Kernels without an exact finite-dimensional state (RBF, rational quadratic, periodic, linear, and products) raise an error of class gaussianprocesses_state_space_error.

Stability

Stable: from version 1.0.0 this interface changes incompatibly only in a major release, after a deprecation period. Results and options that concern an experimental model class, kernel, or argument follow that interface's tier. See gaussianprocesses-package for the policy. method = "state_space" is experimental.

References

Hartikainen, J., and Särkkä, S. (2010). Kalman filtering and smoothing solutions to temporal Gaussian process regression models. IEEE International Workshop on Machine Learning for Signal Processing, 379–384.

Särkkä, S., and Solin, A. (2019). Applied Stochastic Differential Equations. Cambridge University Press.

See Also

optimize_time_series_gp() to estimate the hyperparameters.

Examples

time <- as.Date("2026-01-01") + 0:23
y <- sin(2 * pi * (0:23) / 12)

model <- fit_time_series_gp(
  time,
  y,
  kernel = time_series_kernel(period = 12),
  noise_variance = 0.05
)
model$time[c("origin", "unit")]

# The same posterior by Kalman smoothing, for a Matern kernel.
kernel <- matern32_kernel(length_scale = 3)
exact <- fit_time_series_gp(time, y, kernel, noise_variance = 0.05)
kalman <- fit_time_series_gp(time, y, kernel, noise_variance = 0.05,
  method = "state_space")
c(logLik(exact), logLik(kalman))

Description

Converts future time values using the training model's original time origin and units, then evaluates the GP posterior: exactly, or by Kalman smoothing for state-space models. Each requested point is labelled as interpolation, backward extrapolation, or forward forecast.

Usage

forecast_gp(
  model,
  future_time,
  interval_level = 0.95,
  include_covariance = FALSE,
  variance_tolerance = sqrt(.Machine$double.eps),
  observation_noise_variance = NULL
)

Arguments

Value

An object of class gaussianprocesses_time_forecast.

Stability

Stable: from version 1.0.0 this interface changes incompatibly only in a major release, after a deprecation period. Results and options that concern an experimental model class, kernel, or argument follow that interface's tier. See gaussianprocesses-package for the policy.

Examples

time <- as.Date("2026-01-01") + 0:23
model <- fit_time_series_gp(
  time,
  sin(2 * pi * (0:23) / 12),
  kernel = time_series_kernel(period = 12),
  noise_variance = 0.05
)

forecast <- forecast_gp(model, as.Date("2026-01-01") + c(10.5, 25, 30))
data.frame(
  region = forecast$region,
  mean = forecast$mean,
  forecast$prediction_interval
)

Description

Summarizes a rolling-origin GP backtest with both point and probabilistic forecast metrics.

Usage

forecast_metrics(backtest)

Arguments

Value

A named list containing RMSE, MAE, standardized RMSE, empirical interval coverage, mean interval width, and mean log predictive density.

Stability

Stable: from version 1.0.0 this interface changes incompatibly only in a major release, after a deprecation period. Results and options that concern an experimental model class, kernel, or argument follow that interface's tier. See gaussianprocesses-package for the policy.

Examples

time <- 0:23
backtest <- rolling_origin_gp(
  time,
  sin(2 * pi * time / 12),
  kernel = time_series_kernel(period = 12),
  initial_window = 12,
  horizon = 3,
  step = 3,
  noise_variance = 0.05
)

forecast_metrics(backtest)

Description

The structure, versioning, and persistence of the models returned by fit_gp(), optimize_gp(), fit_time_series_gp(), fit_sparse_gp(), fit_heteroscedastic_gp(), fit_latent_gp(), optimize_latent_gp(), and optimize_sparse_gp().

Exact models

Objects of class gaussianprocesses_model, from fit_gp(), optimize_gp(), and fit_time_series_gp(), document these fields:

x, y

Training inputs as a matrix and responses.

derivative

For derivative observations, the type of each observation: 0 for a function value and $d$ for $\partial f / \partial x_d$. NULL when every observation is a function value.

kernel, mean

The kernel and mean specifications. A polynomial mean records the centring and scaling of its inputs.

mean_fit

For estimated or marginalized mean coefficients, a list with treatment, method ("ml", "reml", "gaussian", or "vague"), the coefficients $\bar\beta$, the basis matrix basis ($H$), and coefficient_factor, an upper-triangular $G$ with $G^\top G = H^\top C^{-1} H$ plus $B^{-1}$ under a Gaussian prior. NULL for other means. Use coef() and vcov() to read the coefficients.

noise_variance, training_noise_variance, noise_structure

The noise variance as supplied, one value per training observation, and "homoscedastic" or "known_heteroscedastic".

cholesky

Upper-triangular $R$ with $R^\top R = C + \delta I$, where $C = K + \Sigma_\varepsilon$ and $\delta$ is jitter.

alpha

$\alpha = (C + \delta I)^{-1}(y - m(X))$, with $m(X) = H \bar\beta$ for estimated or marginalized coefficients.

jitter

Diagonal jitter $\delta$ that the factorization needed; usually 0. It is a multiple of the mean diagonal of $C$ (see fit_gp()).

n_observations, n_features

Training size and input dimension.

optimization

Optimizer diagnostics, for optimize_gp() models.

time

Time metadata, for fit_time_series_gp() models.

schema_version, package_version

See below.

The kernel matrix $K$ and the observation covariance $C$ are not stored, because they can be recomputed exactly: evaluate_kernel(model$kernel, model$x) is $K$, and adding diag(model$training_noise_variance) gives $C$. Storing only $R$ keeps a model with $n$ observations at about $8 n^2$ bytes. Other fields, such as centered_y, prior_mean, cholesky_attempts, and numerical_policy, are internal and may change without notice.

Sparse and heteroscedastic models

Sparse models (gaussianprocesses_sparse_model) document x, y, kernel, mean, noise_variance, inducing_points, n_inducing, method ("fitc" or "vfe"), approximation ("FITC" or "VFE"), log_marginal_likelihood (the FITC objective or the VFE bound), trace_term, and, for optimize_sparse_gp() models, optimization; their other fields are internal. Heteroscedastic models (gaussianprocesses_heteroscedastic_model) document x, y, kernel, noise_kernel, training_noise_variance, converged, iterations, and history, and contain two exact models, mean_model and noise_model. Both record schema_version and package_version.

Latent models

Latent models (gaussianprocesses_latent_model) document the fields listed in fit_latent_gp(). Like exact models, they store the Cholesky factor of $B = I + W^{1/2} K W^{1/2}$ but not $K$, and record schema_version and package_version; they were added at schema 3.

State-space models

State-space models (gaussianprocesses_state_space_model), from fit_time_series_gp() with method = "state_space", store the training data, kernel, mean, noise, log marginal likelihood, and the smoothed fitted means and variances, but no covariance matrix or factor: every prediction filters and smooths again. They were added at schema 4.

Schema versions

schema_version identifies the structure of a model, and package_version is the version of gaussianprocesses that created it. The current schema is 4: schema 2 added mean_fit, schema 3 added derivative, and schema 4 added the method and trace_term of sparse models. Schema 4 is the format of version 1.0. Models saved by versions that predate schema versions, which stored $K$ and $C$, are schema 0, and models saved before parametric means were added are schema 1. Every function that accepts a model first converts older schemas to schema 4; the converted model gives the same results, and a converted sparse model is a FITC model. A converted schema-0 model has package_version NA, because the creating version is unknown.

A model with a newer or unknown schema, or without the fields its class requires, raises an error of class gaussianprocesses_model_format_error. The condition object has fields found_schema, supported_schemas, and package_version.

Persistence

Models can be saved with saveRDS() and restored with readRDS(); a restored model gives the same predictions, likelihoods, gradients, diagnostics, and seeded draws as the original.

Every 1.x release reads models saved by any earlier version, back to 0.1.0. A 1.x release may introduce a new schema only together with an upgrade from schema 4, so models saved by version 1.0 stay readable and give the same results. A model saved by a newer version than the one reading it raises gaussianprocesses_model_format_error rather than giving wrong results. The tests read models saved by the 0.1.0 release and schema-4 models written before 1.0, and check that they give the results they were saved with. Keeping the training data with the kernel, mean, and noise also reproduces a model exactly: refitting them with the same numerical settings gives the same model.

Stability

Stable: from version 1.0.0 this interface changes incompatibly only in a major release, after a deprecation period. Results and options that concern an experimental model class, kernel, or argument follow that interface's tier. See gaussianprocesses-package for the policy.

Examples

x <- seq(-2, 2, length.out = 10)
model <- fit_gp(x, sin(x), rbf_kernel(), noise_variance = 0.01)

model$schema_version

# C is recomputed from the kernel, and R'R reproduces it.
C <- evaluate_kernel(model$kernel, model$x) +
  diag(model$training_noise_variance)
all.equal(crossprod(model$cholesky), C + diag(model$jitter, nrow(C)))

# A saved model gives identical predictions.
path <- tempfile(fileext = ".rds")
saveRDS(model, path)
identical(predict_gp(readRDS(path), 0.5), predict_gp(model, 0.5))

Description

Summarizes exact LOO predictive diagnostics. The returned statistics are descriptive diagnostics rather than a pass/fail calibration test.

Usage

gp_calibration_diagnostics(model, interval_level = 0.95)

Arguments

Details

PIT values should be assessed as a distribution rather than individually. Under a correctly specified continuous predictive distribution they are uniform in repeated sampling.

Value

A named list of calibration summaries.

Stability

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

Examples

x <- seq(-2, 2, length.out = 20)
model <- fit_gp(
  x,
  sin(2 * x),
  kernel = rbf_kernel(length_scale = 0.6),
  noise_variance = 0.01
)

gp_calibration_diagnostics(model, interval_level = 0.9)

Description

Scores predicted class probabilities against binary responses with proper scoring rules (log loss and Brier score), discrimination (ROC AUC), and calibration (reliability diagram data and the expected calibration error). gp_holdout_scores() returns these scores for latent models with a Bernoulli likelihood.

Usage

gp_classification_scores(
  observed,
  probability,
  n_bins = 10L,
  log_density = NULL
)

Arguments

Details

log loss

$-\frac1n \sum_i \log p(y_i)$, the negative mean log predictive probability. Proper; lower is better; $\log 2 \approx 0.693$ for $p_i = 1/2$.

Brier score

$\frac1n \sum_i (p_i - y_i)^2$. Proper; lower is better; 0.25 for $p_i = 1/2$.

ROC AUC

The probability that a random event receives a higher forecast than a random non-event, with ties counted as one half: the Mann-Whitney statistic from mid-ranks. 0.5 for uninformative forecasts; NA without both classes. It measures discrimination only and ignores calibration.

reliability

The forecasts are grouped into n_bins equal-width bins of $[0, 1]$. For calibrated forecasts the observed frequency of each bin is close to its mean probability, within about standard_error $= \sqrt{\bar p (1 - \bar p) / n_b}$.

ECE

The expected calibration error $\sum_b \frac{n_b}{n} |\bar p_b - \bar y_b|$. It depends on the binning, and even calibrated forecasts have a positive ECE of the order of the bin standard errors.

Value

An object of class gaussianprocesses_classification_scores with

summary

Named numeric vector: n, log_loss, brier_score, auc, ece, base_rate (the mean response), and mean_probability.

reliability

Data frame with one row per non-empty bin: bin, lower, upper, n, mean_probability, observed_frequency, and standard_error.

pointwise

Data frame with observed (0/1), probability, log_density, and squared_error.

Stability

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

References

Gneiting, T. and Raftery, A. E. (2007). Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association, 102(477), 359–378.

See Also

gp_scores() for Gaussian predictions, predict_latent_gp().

Examples

set.seed(1)
x <- runif(80, 0, 6)
y <- rbinom(80, 1, plogis(2 * sin(x)))
training <- 1:50

model <- fit_latent_gp(x[training], y[training], rbf_kernel(4, 1),
  bernoulli_likelihood())
scores <- gp_holdout_scores(model, x[-training], y[-training], n_bins = 5)
scores
scores$reliability

Description

Scores Poisson-lognormal predictive distributions of counts, those of latent Gaussian-process models with a poisson_likelihood(), with proper scores, a dispersion check, interval coverage, and randomized probability integral transform (PIT) values. gp_holdout_scores() returns these scores for latent Poisson models.

Usage

gp_count_scores(
  observed,
  latent_mean,
  latent_variance,
  exposure = NULL,
  interval_levels = c(0.5, 0.8, 0.9, 0.95),
  seed = NULL,
  quadrature_order = 40L
)

Arguments

Details

Given $f \sim N(\mu, \sigma^2)$, a count $y \sim \mathrm{Poisson}(E e^f)$ has the Poisson-lognormal predictive distribution, with mean $E e^{\mu + \sigma^2 / 2}$ and probabilities, distribution function, and quantiles computed by Gauss-Hermite quadrature.

log score

The mean log predictive probability of the counts. Proper; higher is better.

Dawid-Sebastiani score

$(y - m)^2 / s^2 + 2 \log s$ with the predictive mean $m$ and standard deviation $s$; proper for the first two moments and recommended for counts by Czado, Gneiting, and Held (2009). Lower is better.

dispersion

The mean squared Pearson residual $(y - m)^2 / s^2$. It is about 1 when the predictive variance is right; clearly larger values mean that the counts are overdispersed relative to the model.

coverage

Of the equal-tailed prediction intervals at each level. Counts are discrete, so the intervals cover at least their level: their coverage is conservative, by up to the probability of their end points.

randomized PIT

$u = F(y - 1) + v\, p(y)$ with $v \sim U(0, 1)$: uniform on $[0, 1]$ when the forecasts are calibrated, unlike the plain PIT $F(y)$ of a discrete distribution. The summary reports their mean (1/2), variance (1/12), and the p-value of a Kolmogorov-Smirnov test of uniformity.

Value

An object of class gaussianprocesses_count_scores with

summary

Named numeric vector: n, log_score, dawid_sebastiani, dispersion, pit_mean, pit_variance, and pit_ks_p_value.

calibration

Data frame with one row per interval level: the level, empirical coverage, and mean width (upper minus lower count).

pointwise

Data frame with observed, predictive_mean, predictive_variance, log_score, dawid_sebastiani, pearson_residual, and pit.

Stability

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

References

Czado, C., Gneiting, T., and Held, L. (2009). Predictive model assessment for count data. Biometrics, 65(4), 1254–1261.

See Also

gp_scores(), gp_classification_scores(), predict_latent_gp().

Examples

set.seed(1)
x <- seq(0, 10, length.out = 60)
counts <- rpois(60, exp(1 + sin(x)))
training <- seq(1, 60, by = 2)

model <- fit_latent_gp(x[training], counts[training], rbf_kernel(0.5, 1.5),
  poisson_likelihood(), mean = constant_mean(1))
gp_holdout_scores(model, x[-training], counts[-training], seed = 1)

Description

Combines exact leave-one-out predictive diagnostics, uncertainty-calibration summaries, and numerical-conditioning diagnostics. The components are kept separate so predictive misfit is not automatically attributed to numerical instability, and numerical instability is not mistaken for model misfit.

Usage

gp_diagnostics(
  model,
  interval_level = 0.95,
  condition_tolerance = sqrt(.Machine$double.eps)
)

Arguments

Value

An object of class gaussianprocesses_diagnostics.

Stability

Stable: from version 1.0.0 this interface changes incompatibly only in a major release, after a deprecation period. Results and options that concern an experimental model class, kernel, or argument follow that interface's tier. See gaussianprocesses-package for the policy. Its numerical component comes from the experimental gp_numerical_diagnostics().

Examples

x <- seq(-2, 2, length.out = 20)
model <- fit_gp(
  x,
  sin(2 * x),
  kernel = rbf_kernel(length_scale = 0.6),
  noise_variance = 0.01
)

gp_diagnostics(model)

Description

Approximates the uncertainty of a model's hyperparameters from the curvature of the log marginal likelihood at its optimum, and flags hyperparameters that the data identify poorly.

Usage

gp_hyperparameter_uncertainty(
  model,
  step = 1e-04,
  interval_level = 0.95,
  condition_threshold = 1e+06,
  log_scale_threshold = log(2)
)

Arguments

Details

The hyperparameters are those optimize_gp() estimated (all kernel parameters and the noise for other models), on its coordinate scale $\eta$: logarithms of positive parameters and real parameters themselves. The observed information $I = -\nabla^2 \ell(\hat\eta)$ is computed by central differences of the analytical gradient (log_marginal_likelihood_gradient()), one refit per displaced coordinate. Its inverse is the covariance of the Laplace approximation with a flat prior on the coordinate scale:

$\widehat{\mathrm{Cov}}(\hat\eta) \approx I^{-1}.$

Standard errors and correlations are reported on the coordinate scale, and approximate intervals $\hat\eta \pm z\,\mathrm{se}$ are mapped back to the natural scale, so that intervals for positive parameters stay positive. Estimated mean coefficients are profiled out (the likelihood is the profile, restricted, or marginal likelihood of the mean), and the approximation conditions on them.

The approximation needs a maximum: gradient reports the gradient at the estimate, which should be close to zero. Some combinations of hyperparameters can be nearly unidentifiable: under infill asymptotics, only $\sigma^2 / \ell^{2\nu}$ is consistently estimable for a Matérn kernel of smoothness $\nu$ (Zhang, 2004). The likelihood is then nearly flat along a ridge: the information is nearly singular in absolute terms, though not always in relative ones, because the curvature across the ridge can be small too. poorly_identified is therefore TRUE when the information is not positive definite, when its condition number exceeds condition_threshold, or when the standard deviation along the least determined combination of the log-scale parameters, weakest_log_scale_sd, exceeds log_scale_threshold. weakest_direction is the eigenvector of the smallest eigenvalue, the direction along the ridge; for a Matérn kernel under infill it is close to $(2\nu, 1) / \sqrt{4\nu^2 + 1}$ in (log variance, log length scale). Real-valued parameters, such as changepoint locations, have no common scale with the others and enter only the first two criteria. gp_profile_likelihood() shows the likelihood along one parameter without the quadratic approximation.

Value

An object of class gaussianprocesses_hyperparameter_uncertainty with

parameters

Data frame: parameter, estimate, coordinate ("log" or "identity"), coordinate_estimate, standard_error on the coordinate scale, and lower and upper on the natural scale.

covariance, correlation

On the coordinate scale; NULL when the information is not positive definite.

information, eigenvalues, condition_number, weakest_log_scale_sd, poorly_identified, weakest_direction

The observed information and its diagnostics.

gradient

The gradient at the estimate.

Stability

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

References

Zhang, H. (2004). Inconsistent estimation and asymptotically equal interpolations in model-based geostatistics. Journal of the American Statistical Association, 99(465), 250–261.

See Also

gp_profile_likelihood(), and the article on hyperparameter uncertainty on the documentation site.

Examples

set.seed(1)
x <- seq(0, 10, length.out = 60)
y <- sin(x) + rnorm(60, sd = 0.2)
model <- optimize_gp(x, y, rbf_kernel(), noise_variance = 0.1, n_starts = 1)

uncertainty <- gp_hyperparameter_uncertainty(model)
uncertainty
uncertainty$correlation

Description

Observation models $p(y_i \mid f_i)$ that link the latent Gaussian process $f$ to the responses. Like kernel specifications, they are plain data with a class, a print method, and constrained parameters.

Usage

gaussian_likelihood(variance = 1)

bernoulli_likelihood(link = c("logit", "probit"))

poisson_likelihood(link = "log")

Arguments

Details

gaussian_likelihood()

$y_i \sim N(f_i, \sigma^2)$, the observation model of fit_gp() with noise_variance $\sigma^2$.

bernoulli_likelihood()

$\Pr(y_i = 1) = \sigma(f_i)$ with the logistic function, $\log p = y f - \log(1 + e^f)$, or $\Phi(f_i)$ with the probit link. Responses may be 0/1 numbers, logical values (TRUE is 1), or a two-level factor whose second level is 1, as in stats::glm(). Multi-class (softmax) likelihoods are not supported.

poisson_likelihood()

$y_i \sim \mathrm{Poisson}(E_i e^{f_i})$ with exposure $E_i$, so $\log p = y (f + \log E) - E e^f - \log y!$. The exposure is data, supplied where the likelihood is used, not part of the specification. Counts more variable than a Poisson process with a smooth rate allows need negative-binomial or zero-inflated likelihoods, which are not yet supported; gp_count_scores() checks for such overdispersion.

Every likelihood provides, vectorized over observations, the log density and its first three derivatives in $f_i$ (evaluate_likelihood()), derivatives in its parameters, the predictive distribution of a new response when $f_*$ is Gaussian (likelihood_predictive()), and whether the log density is concave in $f$. All three are log-concave, which makes Laplace approximations well behaved. Gaussian expectations without a closed form use Gauss-Hermite quadrature: nodes from the Golub-Welsch eigenvalue method, weights from Christoffel sums, cached by order.

The probit derivatives use the inverse Mills ratio $r = \varphi(z) / \Phi(z)$ computed in log space, with $z = (2y - 1) f$. For $z < -5$, where $r$ grows like $-z$ and $z + r$ cancels, they use the continued fraction of the Mills ratio instead. The logistic derivatives compute $\pi$ and $1 - \pi$ separately. All derivatives stay finite and accurate for large $|f|$.

Latent models with any of these likelihoods are fitted with fit_latent_gp() and optimize_latent_gp(); fit_gp() remains the exact Gaussian model.

Value

An object of class gaussianprocesses_likelihood.

Stability

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

References

Golub, G. H. and Welsch, J. H. (1969). Calculation of Gauss quadrature rules. Mathematics of Computation, 23(106), 221–230.

Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press. Chapter 3.

See Also

evaluate_likelihood(), likelihood_predictive(), fit_latent_gp().

Examples

gaussian_likelihood(variance = 0.1)
bernoulli_likelihood("probit")
poisson_likelihood()

evaluate_likelihood(bernoulli_likelihood(), c(TRUE, FALSE), f = c(1, 1))

Description

Mean specifications separate the deterministic part of the model from the covariance kernel. zero_mean() is zero everywhere and constant_mean() is a constant, fixed or with an estimated or marginalized value. For regression trends see linear_mean(), polynomial_mean(), and basis_mean().

Usage

zero_mean()

constant_mean(value = 0)

Arguments

Value

An object of class gaussianprocesses_mean.

Stability

Stable: from version 1.0.0 this interface changes incompatibly only in a major release, after a deprecation period. Results and options that concern an experimental model class, kernel, or argument follow that interface's tier. See gaussianprocesses-package for the policy.

Examples

zero_mean()
constant_mean(2.5)

evaluate_mean(constant_mean(2.5), 1:3)

# Estimate the constant instead.
x <- seq(0, 3, length.out = 15)
model <- fit_gp(x, 4 + sin(2 * x), rbf_kernel(), noise_variance = 0.01,
  mean = constant_mean(estimate_coefficients()))
coef(model)

Description

Methods for stats::predict(), stats::logLik(), stats::nobs(), stats::fitted(), stats::residuals(), and summary() on exact (gaussianprocesses_model), sparse (gaussianprocesses_sparse_model), heteroscedastic (gaussianprocesses_heteroscedastic_model), latent (gaussianprocesses_latent_model), and state-space (gaussianprocesses_state_space_model) models.

Usage

## S3 method for class 'gaussianprocesses_model'
predict(object, newdata, ...)

## S3 method for class 'gaussianprocesses_sparse_model'
predict(object, newdata, ...)

## S3 method for class 'gaussianprocesses_state_space_model'
predict(object, newdata, ...)

## S3 method for class 'gaussianprocesses_heteroscedastic_model'
predict(object, newdata, ...)

## S3 method for class 'gaussianprocesses_latent_model'
predict(object, newdata, ...)

## S3 method for class 'gaussianprocesses_model'
logLik(object, ...)

## S3 method for class 'gaussianprocesses_sparse_model'
logLik(object, ...)

## S3 method for class 'gaussianprocesses_state_space_model'
logLik(object, ...)

## S3 method for class 'gaussianprocesses_heteroscedastic_model'
logLik(object, ...)

## S3 method for class 'gaussianprocesses_latent_model'
logLik(object, ...)

## S3 method for class 'gaussianprocesses_model'
nobs(object, ...)

## S3 method for class 'gaussianprocesses_sparse_model'
nobs(object, ...)

## S3 method for class 'gaussianprocesses_state_space_model'
nobs(object, ...)

## S3 method for class 'gaussianprocesses_heteroscedastic_model'
nobs(object, ...)

## S3 method for class 'gaussianprocesses_latent_model'
nobs(object, ...)

## S3 method for class 'gaussianprocesses_model'
fitted(object, ...)

## S3 method for class 'gaussianprocesses_sparse_model'
fitted(object, ...)

## S3 method for class 'gaussianprocesses_state_space_model'
fitted(object, ...)

## S3 method for class 'gaussianprocesses_heteroscedastic_model'
fitted(object, ...)

## S3 method for class 'gaussianprocesses_latent_model'
fitted(object, ...)

## S3 method for class 'gaussianprocesses_model'
residuals(object, ...)

## S3 method for class 'gaussianprocesses_sparse_model'
residuals(object, ...)

## S3 method for class 'gaussianprocesses_state_space_model'
residuals(object, ...)

## S3 method for class 'gaussianprocesses_heteroscedastic_model'
residuals(object, ...)

## S3 method for class 'gaussianprocesses_model'
summary(object, ...)

## S3 method for class 'gaussianprocesses_sparse_model'
summary(object, ...)

## S3 method for class 'gaussianprocesses_state_space_model'
summary(object, ...)

## S3 method for class 'gaussianprocesses_heteroscedastic_model'
summary(object, ...)

## S3 method for class 'gaussianprocesses_latent_model'
summary(object, ...)

## S3 method for class 'summary.gaussianprocesses_model'
print(x, ...)

Arguments

Details

predict() returns exactly what the package's prediction function for the model returns. For models from fit_time_series_gp(), newdata are time values and predict() calls forecast_gp(); this includes every state-space model with time metadata.

logLik() returns the log marginal likelihood $\log p(y \mid X, \hat\theta)$: log_marginal_likelihood() for exact and state-space models and, for sparse models, the FITC approximate marginal likelihood or the VFE lower bound. This is a type-II likelihood, with the latent function integrated out and the hyperparameters $\hat\theta$ plugged in. Its df attribute counts the hyperparameters estimated from the data: the parameters optimized by optimize_gp() or optimize_sparse_gp(), and 0 for models whose hyperparameters were supplied, plus the mean coefficients when they are estimated. For restricted likelihoods (REML estimates and vague coefficient priors) the nobs attribute is $n - p$, as for restricted likelihoods elsewhere in R, and such likelihoods are comparable only between models with the same mean. stats::AIC() and stats::BIC() therefore work, but comparing kernels by them is a heuristic: it is not an ordinary maximum-likelihood comparison, and the penalty does not account for the flexibility of the kernel. Heteroscedastic models have no logLik() method, because their noise function is a plug-in estimate rather than a likelihood maximizer. For latent models it is the Laplace approximation to the marginal likelihood, and df counts the parameters optimized by optimize_latent_gp() or the mean coefficients estimated by fit_latent_gp(). For state-space models df counts the parameters optimized by optimize_time_series_gp().

fitted() returns the posterior mean of the latent function at the training inputs, and residuals() returns the observed responses minus those fitted values. For exact models the fitted values are $y - (\Sigma_\varepsilon + \delta I)\alpha$, which equals $m(X) + K\alpha$ without forming $K$; $\delta$ is the numerical jitter. For state-space models they are the Kalman-smoother means. For latent models fitted() returns the posterior mode $\hat f$ of the latent function, and there is no residuals() method, because the responses are not on the scale of $f$.

summary() collects the training size and dimension, kernel, mean, noise, log marginal likelihood, optimization status, and numerical status, and prints them compactly.

Value

predict(): a prediction object. logLik(): an object of class logLik. nobs(): the number of training observations. fitted() and residuals(): numeric vectors with one value per training observation. summary(): an object of class summary.gaussianprocesses_model.

Stability

Stable: from version 1.0.0 this interface changes incompatibly only in a major release, after a deprecation period. Results and options that concern an experimental model class, kernel, or argument follow that interface's tier. See gaussianprocesses-package for the policy.

Examples

x <- seq(-2, 2, length.out = 20)
y <- sin(2 * x) + 0.1 * cos(7 * x)
model <- optimize_gp(
  x,
  y,
  kernel = rbf_kernel(),
  noise_variance = 0.05,
  n_starts = 1
)

summary(model)
predict(model, c(-1, 0.5))$mean
logLik(model)
AIC(model)
head(residuals(model))

Description

Reports numerical quantities separately from statistical fit diagnostics. In particular, observation-noise variance and numerical jitter remain distinct.

Usage

gp_numerical_diagnostics(
  model,
  condition_tolerance = sqrt(.Machine$double.eps)
)

Arguments

Value

A named list of numerical diagnostics, including:

observation_noise_variance

The statistical noise variance.

numerical_jitter

The diagonal jitter $\delta$ that the factorization of the training covariance $C$ needed.

covariance_scale

The mean of the diagonal of $C$.

relative_jitter

numerical_jitter / covariance_scale.

cholesky_attempts

Factorization attempts, 1 if no jitter was needed.

reciprocal_condition_number, condition_number

The 1-norm condition estimate of $C + \delta I$ and its reciprocal.

potentially_ill_conditioned

Whether the reciprocal condition number is below condition_tolerance.

observation_types

A data frame with one row per observation type ("value", or "derivative x1", ... for derivative observations): the number of observations, the mean diagonal of their block of $C$, and the reciprocal condition number of that block with the jitter. Derivative blocks scale with $1 / \ell^2$, so their conditioning can differ from that of the values.

Stability

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

Examples

# Duplicated inputs without observation noise need numerical jitter, which is
# reported separately from the (zero) noise variance.
model <- fit_gp(c(0, 0, 1), c(1, 1, 0), kernel = rbf_kernel(), noise_variance = 0)

diagnostics <- gp_numerical_diagnostics(model)
diagnostics[c(
  "observation_noise_variance",
  "numerical_jitter",
  "relative_jitter",
  "cholesky_attempts"
)]