This is the first release on CRAN and the first since 0.1.0.
Breaking changes. The names renamed after 0.1.0 are removed, without the release of deprecation warnings that NEWS promised: version 1.0.0 is the first release after 0.1.0, and the stable API starts without deprecated names (#40).
fit_gp_time() is removed; use
fit_time_series_gp().tolerance argument of fit_gp(),
optimize_gp(), fit_sparse_gp(),
sample_gp_prior(), sample_gp_posterior(), and
simulate_gp_data() is removed; use
symmetry_tolerance. In sample_gp_posterior(),
which passes ... to its methods, an old
tolerance argument is now ignored like any other unknown
argument.fit_heteroscedastic_gp(), tolerance and
numerical_tolerance are removed; use
convergence_tolerance and
symmetry_tolerance.predict_heteroscedastic_gp() no longer returns
noise_variance; use
observation_noise_variance.CRAN preparation (#40).
citation("gaussianprocesses") gives a citation from
inst/CITATION.DESCRIPTION has no BugReports field:
GitLab serves the issue list only to logged-in users, and R’s incoming
check flags the public work items page. The README links the
tracker.cran-check job that
installs every suggested package and checks the built tarball with its
tests, and a manual text-checks job for spelling and URLs.
RELEASING.md describes CRAN submission and resubmission.Stability tiers and the versioning policy for 1.0.
?gaussianprocesses sets out the
policy. From 1.0.0, stable interfaces change incompatibly only in a
major release, after a deprecation period, and experimental interfaces
may change in a minor release, with a NEWS entry.optimize_time_series_gp() and
fit_time_series_gp(method = "state_space"));fit_gp(derivative = )) and
predict_gradient_gp();Validation against independent references.
inst/validation/references.R defines 35 references, at
least one for every inference method: other implementations
(DiceKriging, kernlab, nlme, gplite, and LAPACK’s pivoted Cholesky
decomposition), brute-force computations (dense linear algebra, refits,
numerical integration and differentiation, optim()), and
closed forms. Each has a tolerance and, for other packages, the version
it was set with. A test runs every reference, skipping those whose
package is not installed, and checks that every exported function
belongs to a validated method or is listed as performing no
inference.inst/benchmarks and the gplite, nlme, and DiceKriging
comparison tests are replaced by references. CI now installs DiceKriging
and kernlab, so every reference runs there.datasets::co2, a
classifier for versicolor and virginica irises against logistic
regression, and a Poisson model of datasets::discoveries
against a quadratic Poisson regression.State-space (Kalman) inference for time series.
fit_time_series_gp(method = "state_space") computes the
exact GP posterior and log marginal likelihood in O(n) time by Kalman
filtering and Rauch-Tung-Striebel smoothing, for Matérn-1/2, 3/2, and
5/2 kernels, their sums, and scaled versions of them. The kernel
registry supplies their stochastic differential equations, and the
transition matrices are in closed form. Results agree with exact
inference to a relative 1e-8 on irregular and tied inputs. Other kernels
raise an error of class
gaussianprocesses_state_space_error.forecast_gp(), predict(),
logLik(), fitted(), residuals(),
summary(), log_marginal_likelihood(), and
rolling_origin_gp(method = "state_space") accept the new
models, and sample_gp_posterior() draws by forward
filtering, backward sampling.optimize_time_series_gp() estimates the
hyperparameters of a time-series GP with either method; the state-space
likelihood is optimized with finite-difference gradients.inst/benchmarks/state-space-scaling.R times the method
from 1000 to 100 000 observations, and the time-series vignette compares
state-space and exact forecasts.Inducing-point selection and optimization for sparse Gaussian processes.
select_inducing_points(method = "variance", kernel = ...)
performs greedy conditional-variance selection, the pivot order of a
pivoted Cholesky decomposition of K(X, X), in O(n m^2) time without
forming K(X, X). The sparse fitters accept
selection = "variance".optimize_sparse_gp(optimize_inducing = TRUE) optimizes
the inducing inputs jointly with the hyperparameters, with analytical
gradients from the kernel input derivatives. It is recommended for VFE;
with FITC it overfits.$inducing_collapse, and warn
(class gaussianprocesses_inducing_collapse_warning).Variational sparse Gaussian processes, gradients, and sparse sampling.
fit_sparse_gp() gains method = "vfe", the
collapsed variational bound of Titsias (2009), which lies below the
exact log marginal likelihood and equals it when every input is an
inducing point. VFE predicts with the deterministic training
conditional. FITC stays the default.method; models now use schema 4, and older sparse
models are read as FITC.optimize_sparse_gp() optimizes kernel parameters
and noise for VFE (the default) or FITC with fixed inducing points,
using the optimizer of optimize_gp() and analytical
gradients computed in O(n m^2) time without n x n matrices.sample_gp_posterior() draws from FITC and VFE
posteriors, and predict_sparse_gp() predicts marginal
variances in blocks (block_size), with memory linear in the
number of prediction inputs.Poisson count models with latent Gaussian processes.
estimate_coefficients("ml"), the default of
linear_mean()): fit_latent_gp() maximizes the
approximate marginal likelihood over them with its analytical gradient,
and optimize_latent_gp() estimates them jointly with the
hyperparameters as mean.<name>. A constant mean
estimated this way is the log level of a count model.
coef() returns the coefficients of latent models.predict_latent_gp() returns, for Poisson likelihoods,
rate_interval, the exact interval of the lognormal rate,
and prediction_interval, equal-tailed intervals for new
counts from the Poisson-lognormal predictive distribution, which are
conservative because counts are discrete; for Gaussian likelihoods it
returns the prediction_interval of the exact model.
sample_gp_posterior() returns rate draws.gp_count_scores() scores count forecasts with the
log score, Dawid-Sebastiani score, Pearson dispersion (an overdispersion
check), coverage of the count intervals, and randomized PIT values;
gp_holdout_scores() returns them for latent Poisson
models.datasets::discoveries and
reports a simulation: rate intervals cover the true rate at their
nominal level with known hyperparameters, count intervals are
conservative, randomized PIT values are uniform, and the dispersion
check detects overdispersion.gplite.Binary classification with latent Gaussian processes.
gp_classification_scores() scores probability
forecasts with the log loss, Brier score, ROC AUC, reliability bins, and
expected calibration error. gp_holdout_scores() returns
these scores for latent models with a Bernoulli likelihood, with the log
loss from predictive log probabilities, and the scores of the Gaussian
predictive for a Gaussian likelihood.predict_latent_gp() returns
probability_interval, the interval of the class probability
itself, and sample_gp_posterior() returns
probability draws for Bernoulli likelihoods.optimize_latent_gp() warns (class
gaussianprocesses_bound_warning) when a parameter ends at a
bound of the search and records it in
$optimization$at_bound. With separable classes the signal
variance grows to a bound or to a very large value.gplite package (now
suggested) to 1e-6 in the mode, approximate marginal likelihood, and
latent predictions.Latent Gaussian-process models with the Laplace approximation, for responses observed through a likelihood such as binary labels or counts.
fit_latent_gp(x, y, kernel, likelihood, ...) finds
the posterior mode by Newton’s method on the well-conditioned matrix B =
I + W^(1/2) K W^(1/2) (Rasmussen and Williams, Algorithm 3.1), with the
package’s jitter policy, a step-halving line search, and a convergence
test on the Newton decrement. $laplace reports the
iterations, step halvings, decrement, and gradient norm, and a warning
is raised if the iterations do not converge. Likelihoods that are not
log-concave are refused. The Newton update is written so that it stays
accurate for counts in the thousands, where the usual form loses all
digits.log_marginal_likelihood() and
log_marginal_likelihood_gradient() return the Laplace
approximation and its analytical gradient, including the implicit
dependence of the mode on the hyperparameters, for kernel and likelihood
parameters.optimize_latent_gp() optimizes the hyperparameters
with the optimizer of optimize_gp(), which now shares its
setup, multi-start, and diagnostics code, and warm-starts each Newton
iteration from the previous mode.predict_latent_gp() gives the latent posterior
(Algorithm 3.2) and the predictive distribution of new responses through
the likelihood, including class probabilities.
sample_gp_posterior() draws latent values from the Gaussian
approximation and observations from the likelihood.
predict(), logLik(), nobs(),
fitted(), and summary() have methods for
latent models.gaussian_likelihood() the approximation is exact,
and the model equals fit_gp() to 1e-10.
?fit_latent_gp reports the error of the approximation
against tensor-product quadrature for up to three observations.Likelihoods for non-Gaussian responses, the foundation for latent Gaussian-process models.
gaussian_likelihood(),
bernoulli_likelihood() (logit and probit links), and
poisson_likelihood() (log link, with exposure) build
likelihood specifications: plain data with a print method, like kernel
specifications. Binary responses may be 0/1, logical, or a two-level
factor whose second level is the event, as in glm().evaluate_likelihood() returns the log density and
its first three derivatives in the latent values. They stay finite and
accurate for latent values up to 40 in magnitude: the probit derivatives
use the inverse Mills ratio in log space, and its continued fraction in
the far tail.likelihood_predictive() gives the predictive mean,
variance, class probability, and log density of new responses when the
latent value is Gaussian, in closed form where one exists and otherwise
by adaptive Gauss-Hermite quadrature centred at the mode of the
integrand.fit_gp() remains the exact Gaussian
model (see the API audit).Proper scoring rules, hyperparameter uncertainty, and multi-start diagnostics.
gp_scores() scores Gaussian predictions with the
log score, the CRPS in closed form, interval scores and coverage across
a grid of levels, PIT values, SMSE, and MSLL.
gp_holdout_scores() scores any model on held-out data, and
gp_loo_scores() on its exact leave-one-out
predictions.gp_hyperparameter_uncertainty() computes the
observed information at the estimate from differences of the analytical
gradient on the optimizer scale, with standard errors, correlations,
intervals, and a flag for poorly identified hyperparameters. It detects
the variance-length-scale ridge of a Matérn process observed within one
length scale, and stays silent when the data span twenty length scales.
New gp_profile_likelihood() profiles one
hyperparameter.optimize_gp() records the distinct optima its starts
reached in $optimization$optima, and summary()
reports them.New spectral_mixture_kernel() (Wilson and Adams,
2013): a sum of components w exp(-2 pi^2 tau^2 v) cos(2 pi tau mu), one
weight, frequency, and spectral variance per component and input
dimension. A component with frequency 0 is an RBF kernel.
initialize_spectral_mixture() places the components
at the peaks of the data’s spectrum: the periodogram for regular inputs
and a Lomb-Scargle periodogram for irregular ones, with a mixture of
Gaussians and a uniform background, for the noise floor, fitted by EM.
On noisy sums of two sinusoids it recovers both frequencies within a
fraction of a frequency-resolution bin, for regular and irregular
sampling.optimize_gp() bounds every frequency by 0 and the
Nyquist frequency of its input column, and its later starts displace
frequencies by one frequency-resolution bin rather than one unit.datasets::co2, by held-out log predictive
density.New
changepoint_kernel(before, after, location, steepness, column)
switches between two kernels around an estimable location along one
input column: (1 - s(x)) k1 (1 - s(x’)) + s(x) k2 s(x’) with a sigmoid
s.
location,
steepness, before., and after.,
and nesting gives several changes.optimize_gp() bounds every changepoint location by the
training range of its column and starts it at evenly spaced quantiles of
that column, because the likelihood is often multimodal in the
location.time_series_kernel() gains changepoint,
which wraps the trend-periodic-local kernel in a changepoint with
independent copies before and after the break.Derivative observations. fit_gp() and
optimize_gp() take derivative, the type of
each observation: 0 for a function value and d for the derivative with
respect to input d.
predict_gp()) and
derivatives (predict_gradient_gp()), posterior draws,
leave-one-out diagnostics, and the likelihood and its gradient use
it.input_parameter_gradient, for every differentiable kernel,
combined for sums, products, scaled kernels, and selections of input
columns.gaussianprocesses_smoothness_error.gp_numerical_diagnostics() reports conditioning per
observation type in observation_types.derivative; older
models are read as before.Parametric mean functions — universal kriging.
linear_mean(), polynomial_mean(), and
basis_mean() define trends m(x) = h(x)’β, and
constant_mean() accepts an unknown constant. Polynomial
bases standardize their inputs to [-1, 1] for conditioning.estimate_coefficients(), with
the profile likelihood or REML), or marginalized under a Gaussian prior
or its vague limit (coefficient_prior()). Marginalized
coefficients add their uncertainty to predictions, posterior draws, and
leave-one-out diagnostics, and a Gaussian prior is included in prior
draws.log_marginal_likelihood(), its gradient, and
optimize_gp() use the profile, restricted, or marginal
likelihood that matches the mean. Coefficients are profiled out in
closed form, and $optimization$likelihood and
$optimization$mean_coefficients record the result.
coef() and vcov() return the coefficients and
their covariance. Rank-deficient bases raise
gaussianprocesses_rank_deficiency_error.nlme::gls() (ML and REML), and
DiceKriging::km() (estimates, profile likelihood, and
universal kriging), and the Gaussian-prior model against the Gaussian
process with the augmented kernel.mean_fit; older
models are read as before.inst/benchmarks/external-dicekriging.R used
DiceKriging length scales √2 times too long. DiceKriging’s Gaussian
covariance has the same form as the RBF kernel, so the length scales are
equal.
New select_dimensions() restricts a kernel to some
input columns, k(x_S, x’_S). Sums of restricted kernels give additive
models such as f(x) = f1(x1) + f2(x2, x3).
kernel2.length_scale.New kernel_input_gradient() and
predict_gradient_gp().
kernel_input_gradient() differentiates a kernel with
respect to its inputs: the covariance between the gradient of the
process and the process.predict_gradient_gp() returns the posterior mean,
variance, and credible intervals of each partial derivative of the
latent function of an exact model, computed in blocks like
predict_gp().gaussianprocesses_smoothness_error.smoothness and their input derivatives (see
CONTRIBUTING.md).gp_diagnostics() computes leave-one-out once instead
of twice, which removes about a third of its run time (1.11 s to 0.78 s
at n = 1500).
Test coverage is 96% of lines, up from 92%, and every file in
R/ is above 89%. A coverage CI job, which runs
in scheduled pipelines or when started by hand, measures it with covr
and fails if a file falls below 85%. The built-in kernels are now
registered on first use, so coverage tools see their code; their results
are unchanged.
Behaviour change:
fit_heteroscedastic_gp() corrects its noise targets for the
leave-one-out uncertainty of the latent function. Squared leave-one-out
residuals estimate the noise variance plus that uncertainty, so the
noise was overestimated. In a 50-replicate simulation study the median
log noise error was +0.24, +0.14, and +0.08 for n = 40, 100, and 200. It
is now -0.06, +0.02, and +0.02, each within one Monte Carlo standard
error of zero, with a lower root-mean-square error at every
n. inst/benchmarks/heteroscedastic-noise.R holds the
study, and the heteroscedastic research note derives the
correction.
Behaviour change:
fit_heteroscedastic_gp() defaults to
max_iterations = 50 instead of 6. In the study no fit
converged within 6 iterations, and every fit that converged did so
within 50. It now warns, with class
gaussianprocesses_convergence_warning, when it stops before
converging; about 2 to 10% of the simulated fits cycle without
converging.
Observation noise is handled consistently.
optimize_gp() accepts one noise variance per
observation. With optimize_noise = FALSE the vector is
fixed. With optimize_noise = TRUE it is a known relative
noise pattern whose overall scale is estimated with an analytical
gradient, and $optimization$noise_parameterization records
which case applied.rolling_origin_gp() accepts one noise variance per
observation and gives each training window and each forecast the values
of its rows.benchmark_sparse_gp() gives both models the same scalar
or per-observation noise and takes prediction-time noise through the new
observation_noise_variance argument.API audit (inst/notes/api-audit.md records the
decision for every exported function). Renamed functions and arguments
keep working for at least one minor release, with a deprecation
warning:
fit_gp_time() is now
fit_time_series_gp().tolerance, which meant different things in different
functions, is now symmetry_tolerance in
fit_gp(), optimize_gp(),
fit_sparse_gp(), sample_gp_prior(),
sample_gp_posterior(), and
simulate_gp_data().fit_heteroscedastic_gp(), tolerance is
now convergence_tolerance and
numerical_tolerance is now
symmetry_tolerance.Predictions share one set of fields with the same meaning,
documented in ?predict_gp.
predict_heteroscedastic_gp() now returns
observation_noise_variance and, with
include_covariance = TRUE,
observation_covariance; its noise_variance
field is a deprecated copy of
observation_noise_variance.
Exact, sparse, and heteroscedastic models have methods for
predict(), logLik(), nobs(),
fitted(), residuals(), and
summary(), so AIC() and BIC()
work too (see ?gp_model_methods). logLik() is
the type-II marginal likelihood; comparing kernels by AIC or BIC is a
heuristic. Heteroscedastic models have no logLik(), because
their noise function is a plug-in estimate. The package now declares
stats in Imports.
The documentation explains the two kernel families:
specifications such as rbf_kernel() for models, and
functions such as kernel_rbf() that return one covariance
matrix.
The project is public: the README no longer asks for an access
token, and BugReports points to the issue tracker page that
GitLab serves to visitors who are not signed in.
Behaviour change: numerical jitter is now
relative to the scale of the covariance being factorized. When a
Cholesky factorization fails, the diagonal jitter is the jitter ladder
(initial_jitter, fallback_jitter,
jitter_multiplier) times the mean of the covariance’s
diagonal, instead of the ladder itself. Results therefore no longer
depend on the units of y: rescaling y by 1e-4
or 1e4, which previously changed the log marginal likelihood, posterior
mean, and function draws by up to 170%, now changes them by at most
about 1e-6 (relative, after undoing the scaling). Fits that need no
jitter, which is most of them, are unchanged; fits that need jitter on
covariances whose mean diagonal is not 1 now use a different jitter.
This applies to every factorization: exact and sparse fits,
optimization, sampling, and simulation.
inst/benchmarks/jitter-scaling.R compares the two
policies.
Jitter is reported consistently.
gp_numerical_diagnostics() returns
relative_jitter and covariance_scale (the mean
diagonal of the training covariance) next to
numerical_jitter; they replace
jitter_relative_to_covariance_scale, which divided by the
larger of 1 and the largest diagonal entry. Sample objects record
relative_jitter and covariance_scale, and
models, diagnostics, and samples print the jitter with its size relative
to the covariance scale.
The log marginal likelihood gradient keeps forming the inverse
covariance from its Cholesky factor with chol2inv().
Against an exact reference it was as accurate as two triangular solves
per parameter for condition numbers up to 1e12, and 4 to 18 times
faster. inst/benchmarks/gradient-numerics.R holds the
study, and vignette("v04-numerical-stability") summarizes
both studies.
Fitted exact models no longer store the kernel matrix and the
observation covariance, only the Cholesky factor of the covariance,
which is all that prediction, likelihoods, gradients, leave-one-out
diagnostics, and sampling need. A model with 2000 observations shrinks
from 91.6 MB to 30.6 MB. The matrices are recomputed exactly when
needed. Code that read model$kernel_covariance or
model$observation_covariance should use
evaluate_kernel(model$kernel, model$x), plus
diag(model$training_noise_variance) for the observation
covariance.
Every fitted model now records schema_version and
package_version, and exact models record the numerical
settings used to fit them. Every function that accepts a model converts
models saved by earlier versions, which give the same results, and
refuses models with an unknown schema with a
gaussianprocesses_model_format_error. The new help page
?gaussianprocesses_model documents the fields of fitted
models, their schema, and how to save them.
predict_gp() no longer forms the covariance between
prediction inputs when only marginal predictions are requested (the
default, include_covariance = FALSE). It evaluates the
prior variances with kernel_diagonal() and the
training-by-prediction cross-covariance in blocks, so memory no longer
grows with the square of the number of prediction inputs. With 500
training inputs, 8000 prediction inputs previously needed about 2 GB and
now need about 65 MB; 50000 now need about 100 MB. Results agree with
the full-covariance path to about 1e-14, and
include_covariance = TRUE is unchanged. The new
block_size argument sets the block size; by default it is
chosen from a memory budget. predict_heteroscedastic_gp()
and forecast_gp() use the same path.
inst/benchmarks/prediction-memory.R measures time and
memory.
Internal: each leaf kernel is now defined once, in a kernel
registry that supplies its parameters, constraints, covariance,
diagonal, and derivatives. Constructors, evaluate_kernel(),
kernel_diagonal(), kernel_gradient(),
parameter paths and updates, printing, sparse FITC models, and
optimize_gp() all read the registry instead of listing
kernel types separately. Results are unchanged: golden tests generated
before the change reproduce every covariance, diagonal, gradient,
parameter path, and printout exactly, and optimizer results and
evaluation counts are identical. optimize_gp() now records
each parameter’s optimizer coordinate in
$optimization$coordinates; for the built-in kernels every
coordinate is the log. CONTRIBUTING.md describes how to add a
kernel.
New sample_gp_prior() and
sample_gp_posterior() draw reproducible functions from a GP
prior and from the latent posterior, and noisy observations from the
posterior predictive distribution. Draws use the Cholesky factor of the
covariance with the package’s jitter policy, never an inverse. Models
with known observation-specific noise take prediction-time noise
variances explicitly, and heteroscedastic models are sampled
conditionally on their plug-in noise estimate. The regression vignette
now shows function draws.
optimize_gp() diagnostics now report
function_evaluations, gradient_evaluations,
and model_fits, per start and summed over starts, and the
printout shows them. Unlike optim()’s own
counts, the function evaluations include the objective
calls made for finite-difference gradients.
inst/benchmarks/optimizer-gradients.R compares analytical
and numerical gradients with these counts.
Every exported function’s help page now has a runnable example.
R CMD check and the documentation site run them
all.
optimize_gp() now uses the exact gradient of the log
marginal likelihood instead of finite differences, which needs one GP
fit per optimizer step rather than two per parameter.
gradient = "numerical" restores the previous behaviour, and
the optimization diagnostics and printout record which gradient was
used.
optimize_gp() now allows 500 L-BFGS-B iterations per
start by default instead of 100, so composite kernels with periodic
components, which can need a few hundred iterations, converge.
New log_marginal_likelihood_gradient() returns that
gradient with respect to log kernel parameters and, for homoscedastic
models, the log noise variance.
New kernel_gradient() returns analytical derivatives
of a covariance matrix with respect to the log of every kernel
parameter, named like
kernel_parameters(kernel, flatten = TRUE). It supports all
kernels, ARD length scales, and sum, product, and scaled kernels. A new
vignette derives each formula.
The documentation site now includes the research notes, a worked time-series example, and installation instructions.
First release.